*   MOVE3.FOR
 
*
      SUBROUTINE MOVE3
*
*--the NM molecules are moved over the time interval DTM
*
      include 'parameter.inc'
*
      DOUBLE PRECISION COL(MNSP,MNSP),MOVT,NCOL,SELT,SEPT,CS(7,MNC,MNSP)
      DOUBLE PRECISION CSS(11,MNSE,MNSP)
      integer mcell(MNM), ncell(MNM)
      integer cellout(MNC), cellx(MNC), wgtcell(MNC)
      double precision commx(MNC), avgcomm, variance, stdev
      integer nremain, nmove, nxmove, mvfar
      save nremain, nmove, nxmove, mvfar
*
      IFT=-1
*--a negative IFT indicates that molecules have not entered at this step
      N=0
100   N=N+1
      IF (N.LE.NM) THEN
        IF (IFT.LT.0) AT=DTM
        IF (IFT.GT.0) AT=RND()*DTM
*--the time step is a random fraction of DTM for entering molecules
150     MOVT=MOVT+1
        MSC=IPL(N)
        MC=ISC(MSC)
	mcell(N) = MC
*--MC is the initial cell number
        XI=PP(1,N)
        IF (XI.LT.CB(1).OR.XI.GT.CB(2)) THEN
          WRITE (*,*) ' MOL ',N,' X COORD OUTSIDE FLOW ',XI
          CALL REMOVE(N)
          GO TO 100
        END IF
        YI=PP(2,N)
        IF (YI.LT.CB(3).OR.YI.GT.CB(4)) THEN
          WRITE (*,*) ' MOL ',N,' Y COORD OUTSIDE FLOW ',YI
          CALL REMOVE(N)
          GO TO 100
        END IF
        ZI=PP(3,N)
        IF (ZI.LT.CB(5).OR.ZI.GT.CB(6)) THEN
          WRITE (*,*) ' MOL ',N,' Z COORD OUTSIDE FLOW ',ZI
          CALL REMOVE(N)
          GO TO 100
        END IF
        DX=PV(1,N)*AT
        DY=PV(2,N)*AT
        DZ=PV(3,N)*AT
        X=XI+DX
        Y=YI+DY
        Z=ZI+DZ
        DO 200 KS=1,2
*--check the surfaces
          IF (ISURF(KS).GT.0) THEN
            IF (ISURF(KS).EQ.1.OR.ISURF(KS).EQ.2) THEN
              L1=LIMS(KS,1)
              YS=(L1-1)*CH
              IF ((ISURF(KS).EQ.1.AND.(YI.GT.YS.AND.Y.LT.YS)).OR.
     &            (ISURF(KS).EQ.2.AND.(YI.LT.YS.AND.Y.GT.YS))) THEN
                XC=XI+(YS-YI)*DX/DY
                ZC=ZI+(YS-YI)*DZ/DY
                IF (XC.LE.CB(1).AND.IB(1).EQ.2) THEN
                  XC=2.*CB(1)-XC
                  PV(1,N)=-PV(1,N)
                END IF
                IF (XC.GE.CB(2).AND.IB(2).EQ.2) THEN
                  XC=2.*CB(2)-XC
                  PV(1,N)=-PV(1,N)
                END IF
                IF (ZC.LE.CB(5).AND.IB(5).EQ.2) THEN
                  ZC=2.*CB(5)-ZC
                  PV(3,N)=-PV(3,N)
                END IF
                IF (ZC.GE.CB(6).AND.IB(6).EQ.2) THEN
                  ZC=2.*CB(6)-ZC
                  PV(3,N)=-PV(3,N)
                END IF
                L2=LIMS(KS,2)
                L3=LIMS(KS,3)
                XSU=(L2-1)*CW
                XSD=L3*CW
                L4=LIMS(KS,4)
                L5=LIMS(KS,5)
                ZSU=(L4-1)*CD
                ZSD=L5*CD
                IF ((XC.GT.XSU.AND.XC.LT.XSD).AND.
     &              (ZC.GT.ZSU.AND.ZC.LT.ZSD)) THEN
*--molecule collides with surface at XC,ZC
                  MCX=(XC-CB(1))/CW+0.99999
                  IF (MCX.LT.1) MCX=1
                  IF (MCX.GT.NCX) MCX=NCX
                  MCZ=(ZC-CB(5))/CD+0.99999
                  IF (MCZ.LT.1) MCZ=1
                  IF (MCZ.GT.NCZ) MCZ=NCZ
                  MCSX=MCX-(L2-1)
                  MCSZ=MCZ-(L4-1)
                  MCS=(MCSZ-1)*(L3-L2+1)+MCSX
                  IF (KS.EQ.2) MCS=MCS+(LIMS(1,3)-LIMS(1,2)+1)
     &                             *(LIMS(1,5)-LIMS(1,4)+1)
*--MCS is the code number of the surface element
                  MC=(MCZ-1)*NCX*NCY+MCX
                  IF (ISURF(KS).EQ.1) MC=MC+(L1-1)*NCX
                  IF (ISURF(KS).EQ.2) MC=MC+(L1-2)*NCX
*--MC is the cell number for the reflected molecule
                  AT=AT*(YS-YI)/DY
                  CALL REFLECT3(N,KS,MCS,XC,YS,ZC,MC)
                  GO TO 150
                END IF
              END IF
            END IF
            IF (ISURF(KS).EQ.3.OR.ISURF(KS).EQ.4) THEN
              L1=LIMS(KS,1)
              XS=(L1-1)*CW
              IF ((ISURF(KS).EQ.3.AND.(XI.GT.XS.AND.X.LT.XS)).OR.
     &            (ISURF(KS).EQ.4.AND.(XI.LT.XS.AND.X.GT.XS))) THEN
                YC=YI+(XS-XI)*DY/DX
                ZC=ZI+(XS-XI)*DZ/DX
                IF (YC.LE.CB(3).AND.IB(3).EQ.2) THEN
                  YC=2.*CB(3)-YC
                  PV(2,N)=-PV(2,N)
                END IF
                IF (YC.GE.CB(4).AND.IB(4).EQ.2) THEN
                  YC=2.*CB(4)-YC
                  PV(2,N)=-PV(2,N)
                END IF
                IF (ZC.LE.CB(5).AND.IB(5).EQ.2) THEN
                  ZC=2.*CB(5)-ZC
                  PV(3,N)=-PV(3,N)
                END IF
                IF (ZC.GE.CB(6).AND.IB(6).EQ.2) THEN
                  ZC=2.*CB(6)-ZC
                  PV(3,N)=-PV(3,N)
                END IF
                L2=LIMS(KS,2)
                L3=LIMS(KS,3)
                YSU=(L2-1)*CH
                YSD=L3*CH
                L4=LIMS(KS,4)
                L5=LIMS(KS,5)
                ZSU=(L4-1)*CD
                ZSD=L5*CD
                IF ((YC.GT.YSU.AND.YC.LT.YSD).AND.
     &              (ZC.GT.ZSU.AND.ZC.LT.ZSD)) THEN
*--molecule collides with surface at YC, ZC
                  MCY=(YC-CB(3))/CH+0.99999
                  IF (MCY.LT.1) MCY=1
                  IF (MCY.GT.NCY) MCY=NCY
                  MCZ=(ZC-CB(5))/CD+0.99999
                  IF (MCZ.LT.1) MCZ=1
                  IF (MCZ.GT.NCZ) MCZ=NCZ
                  MCSY=MCY-(L2-1)
                  MCSZ=MCZ-(L4-1)
                  MCS=(MCSZ-1)*(L3-L2+1)+MSCY
                  IF (KS.EQ.2) MCS=MCS+(LIMS(1,3)-LIMS(1,2)+1)
     &                             *(LIMS(1,5)-LIMS(1,4)+1)
*--MCS is the code number of the surface element
                  MC=(MCZ-1)*NCX*NCY+MCY
                  IF (ISURF(KS).EQ.3) MC=(MC-1)*NCX+L1
                  IF (ISURF(KS).EQ.4) MC=(MC-1)*NCX+L1-1
*--MC is the cell number for the reflected molecule
                  AT=AT*(XS-XI)/DX
                  CALL REFLECT3(N,KS,MCS,XS,YC,ZC,MC)
                  GO TO 150
                END IF
              END IF
            END IF
            IF (ISURF(KS).EQ.5.OR.ISURF(KS).EQ.6) THEN
              L1=LIMS(KS,1)
              IF (L1.LE.NCZ) THEN
                ZS=(L1-1)*CD
              ELSE
                ZS=CB(6)
                L1=L1-1
              END IF
              IF ((ISURF(KS).EQ.5.AND.(ZI.GT.ZS.AND.Z.LT.ZS)).OR.
     &            (ISURF(KS).EQ.6.AND.(ZI.LT.ZS.AND.Z.GT.ZS))) THEN
                YC=YI+(ZS-ZI)*DY/DZ
                XC=XI+(ZS-ZI)*DX/DZ
                IF (XC.LE.CB(1).AND.IB(1).EQ.2) THEN
                  XC=2.*CB(1)-XC
                  PV(1,N)=-PV(1,N)
                END IF
                IF (XC.GE.CB(2).AND.IB(2).EQ.2) THEN
                  XC=2.*CB(2)-XC
                  PV(1,N)=-PV(1,N)
                END IF
                IF (YC.LE.CB(3).AND.IB(3).EQ.2) THEN
                  YC=2.*CB(3)-YC
                  PV(2,N)=-PV(2,N)
                END IF
                IF (YC.GE.CB(4).AND.IB(4).EQ.2) THEN
                  YC=2.*CB(4)-YC
                  PV(2,N)=-PV(2,N)
                END IF
                L2=LIMS(KS,2)
                L3=LIMS(KS,3)
                XSU=(L2-1)*CW
                XSD=L3*CW
                L4=LIMS(KS,4)
                L5=LIMS(KS,5)
                YSU=(L4-1)*CH
                YSD=L5*CH
                IF ((XC.GT.XSU.AND.XC.LT.XSD).AND.
     &              (YC.GT.YSU.AND.YC.LT.YSD)) THEN
*--molecule collides with surface at YC, XC
                  MCY=(YC-CB(3))/CH+0.99999
                  IF (MCY.LT.1) MCY=1
                  IF (MCY.GT.NCY) MCY=NCY
                  MCX=(XC-CB(1))/CW+0.99999
                  IF (MCX.LT.1) MCX=1
                  IF (MCX.GT.NCX) MCX=NCX
                  MCSX=MCX-(L2-1)
                  MCSY=MCY-(L4-1)
                  MCS=(MCSY-1)*(L3-L2+1)+MCSX
                  IF (KS.EQ.2) MCS=MCS+(LIMS(1,3)-LIMS(1,2)+1)
     &                             *(LIMS(1,5)-LIMS(1,4)+1)
*--MCS is the code number of the surface element
                  MC=(MCY-1)*NCX+MCX
                  IF (ISURF(KS).EQ.5) MC=MC+(L1-1)*NCX*NCY
                  IF (ISURF(KS).EQ.6) MC=MC+(L1-2)*NCX*NCY
*--MC is the cell number for the reflected molecule
                  AT=AT*(ZS-ZI)/DZ
                  CALL REFLECT3(N,KS,MCS,XC,YC,ZS,MC)
                  GO TO 150
                END IF
              END IF
            END IF
          END IF
200     CONTINUE
        IF (X.LT.CB(1).OR.X.GT.CB(2)) THEN
          IF (X.LT.CB(1)) K=1
          IF (X.GT.CB(2)) K=2
*--intersection with boundary K
          IF (IB(K).EQ.2) THEN
*--specular reflection from the boundary (eqn (11.7))
            X=2.*CB(K)-X
            PV(1,N)=-PV(1,N)
          ELSE
*--molecule leaves flow
            CALL REMOVE(N)
            GO TO 100
          END IF
        END IF
        IF (Y.LT.CB(3).OR.Y.GT.CB(4)) THEN
          IF (Y.LT.CB(3)) K=3
          IF (Y.GT.CB(4)) K=4
*--intersection with boundary K
          IF (IB(K).EQ.2) THEN
*--specular reflection from the boundary (eqn (11.7))
            Y=2.*CB(K)-Y
            PV(2,N)=-PV(2,N)
          ELSE
*--molecule leaves flow
            CALL REMOVE(N)
            GO TO 100
          END IF
        END IF
        IF (Z.LT.CB(5).OR.Z.GT.CB(6)) THEN
          IF (Z.LT.CB(5)) K=5
          IF (Z.GT.CB(6)) K=6
*--intersection with boundary K
          IF (IB(K).EQ.2) THEN
*--specular reflection from the boundary (eqn (11.7))
            Z=2.*CB(K)-Z
            PV(3,N)=-PV(3,N)
          ELSE
*--molecule leaves flow
            CALL REMOVE(N)
            GO TO 100
          END IF
        END IF
*
        IF (X.LT.CG(1,MC).OR.X.GT.CG(2,MC).OR.Y.LT.CG(3,MC).OR.
     &      Y.GT.CG(4,MC).OR.Z.LT.CG(5,MC).OR.Z.GT.CG(6,MC)) THEN
*--the molecule has moved from the initial cell
          MCX=(X-CB(1))/CW+0.99999
          IF (MCX.LT.1) MCX=1
          IF (MCX.GT.NCX) MCX=NCX
*--MCX is the new cell column (note avoidance of round-off error)
          MCY=(Y-CB(3))/CH+0.99999
          IF (MCY.LT.1) MCY=1
          IF (MCY.GT.NCY) MCY=NCY
*--MCY is the new cell row (note avoidance of round-off error)
          MCZ=(Z-CB(5))/CD+0.99999
          IF (MCZ.LT.1) MCZ=1
          IF (MCZ.GT.NCZ) MCZ=NCZ
*--MCZ is the new cell slice (note avoidance of round-off error)
          MC=(MCZ-1)*NCX*NCY+(MCY-1)*NCX+MCX
        END IF
        MSCX=((X-CG(1,MC))/CW)*(NSCX-.001)+1
        MSCY=((Y-CG(3,MC))/CH)*(NSCY-.001)+1
        MSCZ=((Z-CG(5,MC))/CD)*(NSCZ-.001)+1
        MSC=(MSCZ-1)*NSCX*NSCY+(MSCY-1)*NSCX+MSCX+NSCX*NSCY*NSCZ*(MC-1)
*--MSC is the new sub-cell number
        IF (MSC.LT.1) MSC=1
        IF (MSC.GT.MNSC) MSC=MNSC
        IPL(N)=MSC
        PP(1,N)=X
        PP(2,N)=Y
        PP(3,N)=Z
	ncell(N) = MSC
        GO TO 100
      ELSE IF (IFT.LT.0) THEN
        IFT=1
*--new molecules enter
        CALL ENTER3
        N=N-1
        GO TO 100
      END IF
c	do i = 1, NM
c	  write(*,*) 'mcell = ',mcell(i),'  ncell = ',ncell(i)
c	enddo

	nremain = 0
	nmove = 0
	nxmove = 0
	do i = 1, MNC
	  cellout(i) = 0
	  cellx(i) = 0
	  wgtcell(i) = 0
	enddo

	mvfar = 0
	do i = 1, NM
	  wgtcell(mcell(i)) = wgtcell(mcell(i)) + 1
	  if (mcell(i) .eq. ncell(i)) then
		nremain = nremain + 1
	  else
		nmove = nmove + 1
		cellout(mcell(i)) = cellout(mcell(i)) + 1
	  endif
	  if (mcell(i)+1 .eq. ncell(i)) then
		nxmove = nxmove + 1
		cellx(mcell(i)) = cellx(mcell(i)) + 1
	  endif
	  if (mcell(i) .ne. ncell(i) .and.
     &		is_neighbor(mcell(i),ncell(i),NCX,NCY,NCZ) .eq. 0) then
		mvfar = mvfar + 1
	  endif
	enddo

C	if (NPR .eq. NPS) then
C	  open(4,file='dsmc3.load',form='formatted')
C	  do i = 1, MNC
C		write(4,444) wgtcell(i),cellx(i)	! vertex and edge
C	  enddo
C	  close(4)
C	  stop
C	endif

	avgout = float(nmove) / MNC
	avgx = float(nxmove) / MNC

	maxout = cellout(1)
	maxx = cellx(1)
	do i = 1, MNC
	  if (cellout(i) .gt. maxout) maxout = cellout(i)
	  if (cellx(i) .gt. maxx) maxx = cellx(i)
	enddo

	nsteps = nsteps + 1
	NM_sum = NM_sum + NM
	nre_sum = nre_sum + nremain
	nmov_sum = nmov_sum + nmove
	nxmv_sum = nxmv_sum + nxmove
	Nmovefar = Nmovefar + mvfar

	avgcomm = 0.
	variance = 0.
	do i = 1, MNC
	  commx(i) = cellx(i)
	  avgcomm = avgcomm + commx(i)
	enddo
	avgcomm = avgcomm / MNC
	do i = 1, MNC
	  variance = variance + (commx(i)-avgcomm) * (commx(i)-avgcomm)
	enddo
	variance = variance / MNC
	stdev = sqrt(variance)

#ifdef MOVE_TRACE
	write(*,111) NM,nremain,nmove,nxmove,float(nxmove)/float(nmove)
	write(*,112) nmove,mvfar,float(nmove)/float(NM),
     &		float(mvfar)/float(nmove)
	write(*,122) maxout,avgout,float(maxout)/avgout,
     &		maxx,avgx,float(maxx)/avgx
	write(*,133) avgcomm,variance,stdev
111	format(' NM=',I8,' nremain=',I8,' nmove=',I7,' nxmove=',I7,
     &		' rate_X=',F5.3)
112	format(' nmove=',I8,' mvfar=',I8,' rate_move=',F5.3,
     &		' rate_FARmove=',F5.3)
122	format(' maxout=',I4,' avgout=',F5.1,'(R=',F7.3,') maxx=',I4,' avgx=',
     &		F5.1,'(R=',F5.3,')')
133	format(' COMM: avg=',F7.3,' variance=',F7.3,' std devi=',F7.3)
444	format(' ',I8,' ',I8)
#endif
      RETURN
      END
