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```