cTrkden.for c c

c
c
c
      SUBROUTINE TRKDEN(VF,LXU2D,LYV2D,IFACE)
C
C This subroutine follows a string of particles and
C determines the area delimited by the string and the
C edges of the computational cells.
C
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/tracmn'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
c...npcdim is the maximum number of particles allowed in
c   any segment of any string traversing a cell.
c
      PARAMETER (NPCDIM=100)
      DIMENSION VF(35,45),AREA2(35,45) ! what is the significance of the 
                 ! 35,45 ? They are the NY & NX of case g722 !!!
                 ! VF is volume fraction; area2 is area occupied by
                 ! phase ??
      DIMENSION DY(35),DX(45),XU(45),YV(35),VVOL(35,45)
      DIMENSION XPC(NPCDIM),YPC(NPCDIM),XP(1000),YP(1000)
      INTEGER OFACE
      LOGICAL NEZ
      SAVE
C
      DATA RMIN/1.E-06/
C
C
c..initialise array to contain the computed areas.
C      IF(ISTEP.EQ.1) THEN    ! first step only
c        LDEN1=L0F(DEN1)
c        DO 551 IX=1,NX
c          DO 551 IY=1,NY
c            I=IY+(IX-1)*NY
c            VF(IY,IX)=1.0
c            IF(F(LDEN1+I).LT.RHO2) VF(IY,IX)=0.0
c            VF(1,1)=1.0  ! take outside DO loop
c 551    CONTINUE
c      ENDIF
      DO 5 IY=1,NY           ! all time steps
        DO 5 IX=1,NX
        AREA2(IY,IX)=0.0
  5   CONTINUE
      DO 552 IX=1,NX         ! why do this at every time step ?
        DO 552 IY=1,NY
          I=IY+(IX-1)*NY
          XU(IX)=F(LXU2D+I)
          YV(IY)=F(LYV2D+I)
 552  CONTINUE
c..compute geometrical quantities only for the first
c  time step.
      IF(ISTEP.NE.1) GO TO 6
      DY(1)=YV(1)
      DX(1)=XU(1)
      DO 1 IY=2,NY
        DY(IY)=YV(IY)-YV(IY-1)
        DO 2 IX=2,NX
          DX(IX)=XU(IX)-XU(IX-1)
          VVOL(IY,IX)=DY(IY)*DX(IX)
  2     CONTINUE
  1   CONTINUE
      DO 3 IX=1,NX
  3   VVOL(1,IX)=DY(1)*DX(IX)
      DO 4 IY=1,NY
  4   VVOL(IY,1)=DY(IY)*DX(1)
  6   CONTINUE
C### JZW from GENTRA
      DO 7 IPR=1,NTRACK
        ILOC0=IPFSTR(IPR)
        IXYZP=ILOC0+IPRSTN
        XP(IPR)=F(IXYZP+1)
        YP(IPR)=F(IXYZP+2)
  7   CONTINUE
c
C
c..determine the indices of the cell containing the starting
c  point.
      DO 10 IYIN=1,NY
        IF(YP(1).LT.YV(IYIN)) GO TO 20
 10   CONTINUE
 20   DO 11 IXIN=1,NX
        IF(XP(1).LT.XU(IXIN)) GO TO 21
 11   CONTINUE
 21   CONTINUE
c..initialise string
      IPOUT=1
      XPC(1)=XP(1)
      YPC(1)=YP(1)
c..determine coordinates of the edges of the cell
      XB=XU(IXIN)
      YB=YV(IYIN)
      XL=XB-DX(IXIN)
      YL=YB-DY(IYIN)
C
c..move along the interface looking to see when it leaves.
c  when it does so, use ypc and xpc to store the coordinates
c  of the point at which the string left the cell.
 25   CONTINUE
      OFACE=0
      DO 30 NPTS=1,NPCDIM
c       write(14,*) 'npts= ',npts,'xpc= ',xpc(npts),'ypc= ',ypc(npts)
        XDIF=XPC(NPTS)-XP(IPOUT)
        IF(ABS(XDIF).GT.1.E-8) THEN
          UUM=(YPC(NPTS)-YP(IPOUT))/XDIF
        ELSE
          UUM=0.0
        ENDIF
        UUC=YPC(NPTS)-UUM*XPC(NPTS)
c..check if the track leaves through the west face.
        XOUT=XL
        YOUT=UUM*XOUT+UUC
        IF(YOUT.LE.YB.AND.YOUT.GE.YL.AND.XP(IPOUT).LT.XL) THEN
          OFACE=1
          GO TO 31
        ENDIF
c..check if the track leaves through the east face.
        XOUT=XB
        YOUT=UUM*XOUT+UUC
        IF(YOUT.LE.YB.AND.YOUT.GE.YL.AND.XP(IPOUT).GE.XB) THEN
          OFACE=2
          GO TO 31
        ENDIF
c..check if the track leaves through the south face.
        YOUT=YL
        IF(NEZ(UUM)) THEN
          XOUT=(YOUT-UUC)/UUM
        ELSE
          XOUT=XPC(NPTS)
        ENDIF
        IF(XOUT.LE.XB .AND. XOUT.GE.XL .AND. YP(IPOUT).LT.YL) THEN
          OFACE=3
          GO TO 31
        ENDIF
c..check if the track leaves through the north face.
        YOUT=YB
        IF(NEZ(UUM)) THEN
          XOUT=(YOUT-UUC)/UUM
        ELSE
          XOUT=XPC(NPTS)
        ENDIF
        IF(XOUT.LE.XB.AND.XOUT.GE.XL.AND.YP(IPOUT).GE.YB) THEN
          OFACE=4
          GO TO 31
        ENDIF
        XPC(NPTS+1)=XP(IPOUT)
        YPC(NPTS+1)=YP(IPOUT)
        IPOUT=IPOUT+1
 30   CONTINUE
 31   CONTINUE
c      write(14,*) 'iface= ',iface,'oface= ',oface,ixin,iyin
      NPTS=NPTS+1
      YPC(NPTS)=YOUT
      XPC(NPTS)=XOUT
C
c..calculate the volume of the cell occupied by the fluid
c below the track.
c...  1-west, 2-east, 3-south, 4-north
c..if the track goes straight across the cell.
c
      IF(IFACE.EQ.2.AND. OFACE.EQ.1.OR.
     1  IFACE.EQ.1 .AND.OFACE.EQ.2) THEN
        IF(IFACE.EQ.2) THEN
          AREAA2=AREA(XPC,YPC,NPTS ,YB,-1.0)
        ELSE
          AREAA2=VVOL(IYIN,IXIN)-AREA(XPC,YPC,NPTS,YB,1.0)
        ENDIF
      ELSE
        IF(IFACE.EQ.3.AND.OFACE.EQ.4.OR.
     1IFACE.EQ.4.AND.OFACE.EQ.3 ) THEN
          IF(IFACE.EQ.3) THEN
       AREAA2=VVOL(IYIN,IXIN)-AREA(YPC,XPC,NPTS,XB,1.0)
          ELSE
            AREAA2=AREA(YPC,XPC,NPTS,XL,1.0)
          ENDIF
        ELSE
c..if the track cuts the cell across a corner.
c
          IF(IFACE.EQ.2) THEN
            IF(OFACE.EQ.3) THEN
              AREAA2=VVOL(IYIN,IXIN)-AREA(XPC,YPC,NPTS,YL,1.0)
            ELSE
              AREAA2=AREA(XPC,YPC,NPTS,YB,-1.0)
            ENDIF
          ELSE
            IF(IFACE.EQ.3) THEN
              IF(OFACE.EQ.1) THEN
              AREAA2=VVOL(IYIN,IXIN)-AREA(YPC,XPC,NPTS,XL,-1.0)
              ELSE
                AREAA2=AREA(YPC,XPC,NPTS,XB,1.0)
              ENDIF
            ELSE
              IF(IFACE.EQ.1) THEN
                IF(OFACE.EQ.4) THEN
                  AREAA2=VVOL(IYIN,IXIN)-
     1            AREA(XPC,YPC,NPTS,YB,1.0)
                ELSE
                  AREAA2=AREA(XPC,YPC,NPTS,YL,-1.0)
                ENDIF
              ELSE
                IF(IFACE.EQ.4) THEN
                  IF(OFACE.EQ.2) THEN
                    AREAA2=VVOL(IYIN,IXIN)-
     1              AREA(XPC,YPC,NPTS,YB,1.0)
                  ELSE
                    AREAA2=AREA(XPC,YPC,NPTS,YB,-1.0)
                  ENDIF
                ENDIF
              ENDIF
            ENDIF
          ENDIF
        ENDIF
      ENDIF
C
c..if the track enters and leaves the cell through the
c  same face.
      IF(IFACE.EQ.2 .AND.OFACE.EQ.2) THEN
        IF(YPC(1).GT.YPC(NPTS)) THEN
          AREAA2=VVOL(IYIN,IXIN)-AREA(YPC,XPC,NPTS,XB,-1.0)
        ELSE
          AREAA2=AREA(YPC,XPC,NPTS,XB,1.0)
        ENDIF
      ELSE
        IF(IFACE.EQ.3.AND.OFACE.EQ.3) THEN
          IF(XPC(1).GT.XPC(NPTS)) THEN
            AREAA2=VVOL(IYIN,IXIN)-AREA(XPC,YPC,NPTS,YL,1.0)
          ELSE
            AREAA2=AREA(XPC,YPC,NPTS,YL,-1.0)
          ENDIF
        ELSE
          IF(IFACE.EQ.1 .AND. OFACE.EQ.1) THEN
            IF(YPC(1).GT.YPC(NPTS)) THEN
              AREAA2=AREA(YPC,XPC,NPTS,XL,1.0)
            ELSE
              AREAA2=VVOL(IYIN,IXIN)-AREA(YPC,XPC,NPTS,XL,-1.0)
            ENDIF
          ELSE
            IF(IFACE.EQ.4.AND.OFACE.EQ.4) THEN
              IF(XPC(1).GT.XPC(NPTS)) THEN
                AREAA2=AREA(XPC,YPC,NPTS,YB,-1.0)
              ELSE
                AREAA2=VVOL(IYIN,IXIN)-AREA(XPC,YPC,NPTS,YB,1.0)
              ENDIF
            ENDIF
          ENDIF
        ENDIF
      ENDIF
C
c..sum areas so that they are overlayed to give the correct
c  value
c      write(14,*) 'AREAA2',AREAA2
      IF(AREAA2.GT.VVOL(IYIN,IXIN)-AREA2(IYIN,IXIN)) THEN
        AREA2(IYIN,IXIN)=AREA2(IYIN,IXIN)+AREAA2-VVOL(IYIN,IXIN)
      ELSE
        AREA2(IYIN,IXIN)=AREA2(IYIN,IXIN)+AREAA2
      ENDIF
C
c.. prepare parameters to enter the next cell.
        IF(OFACE.EQ.4) THEN
          IYIN=IYIN+1
          IFACE=3
          YL=YB
        YB=YB+DY(IYIN)
      ELSE
        IF(OFACE.EQ.3) THEN
          IYIN=IYIN-1
          IFACE=4
          YB=YL
          YL=YL-DY(IYIN)
        ELSE
          IF(OFACE.EQ.2) THEN
            IXIN=IXIN+1
            IFACE=1
            XL=XB
            XB=XB+DX(IXIN)
          ELSE
            IF(OFACE.EQ.1) THEN
              IXIN=IXIN-1
              IFACE=2
              XB=XL
              XL=XB-DX(IXIN)
            ENDIF
          ENDIF
        ENDIF
      ENDIF
      YPC(1)=YOUT
      XPC(1)=XOUT
C
c..check whether the last particle on the string has been
c reached.
      IF(IPOUT.GE.NTRACK) GO TO 100
C
c....if it has not, go and keep on moving along the track.
      GO TO 25
 100  CONTINUE
C
c..the last particle on the string has been reached.
c  see what has changed and set volume fraction
c  accordingly.
C
      DO 110 IY=1,NY
      DO 110 IX=1,NX
c..for those cells which were traversed by the string,
c  calculate the fraction of volume occuped by the fluid
c below the track.
      IF(AREA2(IY,IX).GT.VVOL(IY,IX)/100.0) THEN
        VF(IY,IX)= AREA2(IY,IX)/VVOL(IY,IX)
      ELSE
c..set volume fraction for those cells which do not contain
c  any particles.
      IF(VF(IY,IX).LT.0.5) THEN
        VF(IY,IX)=RMIN
      ELSE
        VF(IY,IX)=1.0-RMIN
      ENDIF
      ENDIF
c     WRITE(14,*) 'VF OLD=',VF(IY,IX)
 110  CONTINUE
C
      END
C
c
      FUNCTION AREA(XP,YP,NPTS,YY,SWAP) ! swap = 1.0 or -1.0
c   this function computes the area delimited by a
c   string of particles and the edges of the cell
c   which contains them.
c ???? y and x are not treated in the same way. Why not? 
c                                               And does it matter?
C
      DIMENSION XP(100),YP(100)
C
c.....initialise function to zero.
      AREA=0.0
C
c.....move along the segment calculating area.
      DO 10 IP=1,NPTS-1
        XOLD=XP(IP)
        XNEW=XP(IP+1)
        YOLD=YP(IP)
        YNEW=YP(IP+1)
        AREA=AREA+0.5*(2.0*YY-YNEW-YOLD)*(XNEW-XOLD)*SWAP
 10   CONTINUE
      END
C