*   ENTER3.FOR
*
      SUBROUTINE ENTER3
*
*--new molecules enter at boundaries
*
      include 'parameter.inc'
*
      DOUBLE PRECISION COL(MNSP,MNSP),MOVT,NCOL,SELT,SEPT,CS(7,MNC,MNSP)
*
      DO 100 N=1,6
*--consider each boundary in turn
        IF (IB(N).EQ.1) THEN
          IF (N.LT.3) THEN
            NCS1=NCY
            NCS2=NCZ
          END IF
          IF (N.GT.2.AND.N.LT.5) THEN
            NCS1=NCX
            NCS2=NCZ
          END IF
          IF (N.GT.4) THEN
            NCS1=NCX
            NCS2=NCY
          END IF
          DO 20 NC1=1,NCS1
            DO 10 NC2=1,NCS2
              DO 5 L=1,MNSP
*--consider each species in turn
                VMP=SQRT(2.*BOLTZ*FTMP/SP(5,L))
                IF (N.LT.3) A=BME(N,L)*CH*CD/(FH*FD)+BMR(N,L)
                IF (N.GT.2.AND.N.LT.5) A=BME(N,L)*CW*CD/(FW*FD)+BMR(N,L)
                IF (N.GT.4) A=BME(N,L)*CW*CH/(FW*FH)+BMR(N,L)
                M=A
                BMR(N,L)=A-M
*--M molecules enter, remainder has been reset
                IF (M.GT.0) THEN
                  IF (N.EQ.1.OR.N.EQ.2) THEN
                    IF (ABS(VFX).GT.1.E-6) THEN
                      IF (N.EQ.1) SC=VFX/VMP
                      IF (N.EQ.2) SC=-VFX/VMP
                    END IF
                  END IF
                  IF (N.EQ.3.OR.N.EQ.4) THEN
                    IF (ABS(VFY).GT.1.E-6) THEN
                      IF (N.EQ.3) SC=VFY/VMP
                      IF (N.EQ.4) SC=-VFY/VMP
                    END IF
                  END IF
                  FS1=SC+SQRT(SC*SC+2.)
                  FS2=0.5*(1.+SC*(2.*SC-FS1))
* the above constants are required for the entering distn. of eqn (12.?)
                  DO 4 K=1,M
                    IF (NM.LT.MNM) THEN
                      NM=NM+1
*--NM is now the number of the new molecule
                      IF ((N.LT.3.AND.ABS(VFX).GT.1.E-6).OR.
     &                    (N.GT.2.AND.N.LT.5.AND.ABS(VFY).GT.1.E-6).OR.
     &                    (N.GT.4.AND.ABS(VFZ).GT.1.E-6)) THEN
                        QA=3.
                        IF (SC.LT.-3.) QA=ABS(SC)+1.
2                       U=-QA+2.*QA*RND()
*--U is a potential normalised thermal velocity component
                        UN=U+SC
*--UN is a potential inward velocity component
                        IF (UN.LT.0.) GO TO 2
                        A=(2.*UN/FS1)*EXP(FS2-U*U)
                        IF (A.LT.RND()) GO TO 2
*--the inward normalised vel. component has been selected (eqn (12.?))
                        IF (N.EQ.1) PV(1,NM)=UN*VMP
                        IF (N.EQ.2) PV(1,NM)=-UN*VMP
                        IF (N.EQ.3) PV(2,NM)=UN*VMP
                        IF (N.EQ.4) PV(2,NM)=-UN*VMP
                        IF (N.EQ.5) PV(3,NM)=UN*VMP
                        IF (N.EQ.6) PV(3,NM)=-UN*VMP
                      ELSE
                        IF (N.EQ.1) PV(1,NM)=SQRT(-LOG(RND()))*VMP
                        IF (N.EQ.2) PV(1,NM)=-SQRT(-LOG(RND()))*VMP
                        IF (N.EQ.3) PV(2,NM)=SQRT(-LOG(RND()))*VMP
                        IF (N.EQ.4) PV(2,NM)=-SQRT(-LOG(RND()))*VMP
                        IF (N.EQ.5) PV(3,NM)=SQRT(-LOG(RND()))*VMP
                        IF (N.EQ.6) PV(3,NM)=-SQRT(-LOG(RND()))*VMP
*--for a stationary external gas, use eqn (12.3)
                      END IF
                      IF (N.LT.3) THEN
                        CALL RVELC(PV(2,NM),PV(3,NM),VMP)
                        PV(2,NM)=PV(2,NM)+VFY
                        PV(3,NM)=PV(3,NM)+VFZ
                      END IF
                      IF (N.GT.2.AND.N.LT.5) THEN
                        CALL RVELC(PV(1,NM),PV(3,NM),VMP)
                        PV(1,NM)=PV(1,NM)+VFX
                        PV(3,NM)=PV(3,NM)+VFZ
                      END IF
                      IF (N.GT.4) THEN
                        CALL RVELC(PV(1,NM),PV(2,NM),VMP)
                        PV(1,NM)=PV(1,NM)+VFX
                        PV(2,NM)=PV(2,NM)+VFY
                      END IF
*--a single call of RVELC generates the two normal velocity components
                      IF (ISPR(1,L).GT.0)
     &                     CALL SROT(PR(NM),FTMP,ISPR(1,L))
                      IF (N.EQ.1) PP(1,NM)=CB(1)+0.001*CW
                      IF (N.EQ.2) PP(1,NM)=CB(2)-0.001*CW
                      IF (N.EQ.3) PP(2,NM)=CB(3)+0.001*CH
                      IF (N.EQ.4) PP(2,NM)=CB(4)-0.001*CH
                      IF (N.EQ.5) PP(3,NM)=CB(5)+0.001*CD
                      IF (N.EQ.6) PP(3,NM)=CB(6)-0.001*CD
*--the molecule is moved just off the boundary
                      IF (MNSP.GT.1) IPS(NM)=L
                      IF (N.LT.3) THEN
                        IF (N.EQ.1) MC=(NC2-1)*NCX*NCY+(NC1-1)*NCX+1
                        IF (N.EQ.2) MC=(NC2-1)*NCX*NCY+NC1*NCX
                        PP(2,NM)=CG(3,MC)+RND()*CH
                        PP(3,NM)=CG(5,MC)+RND()*CD
                      END IF
                      IF (N.GT.2.AND.N.LT.5) THEN
                        IF (N.EQ.3) MC=(NC2-1)*NCX*NCY+NC1
                        IF (N.EQ.4) MC=(NC2-1)*NCX*NCY+(NCY-1)*NCX+NC1
                        PP(1,NM)=CG(1,MC)+RND()*CW
                        PP(3,NM)=CG(5,MC)+RND()*CD
                      END IF
                      IF (N.GT.4) THEN
                        IF (N.EQ.5) MC=(NC2-1)*NCX+NC1
                        IF (N.EQ.6) MC=(NCZ-1)*NCX*NCY+(NC2-1)*NCX+NC1
                        PP(1,NM)=CG(1,MC)+RND()*CW
                        PP(2,NM)=CG(3,MC)+RND()*CH
                      END IF
                      MSCX=((PP(1,NM)-CG(1,MC))/CW)*(NSCX-.001)+1
                      MSCY=((PP(2,NM)-CG(3,MC))/CH)*(NSCY-.001)+1
                      MSCZ=((PP(3,NM)-CG(5,MC))/CD)*(NSCZ-.001)+1
                      MSC=(MSCZ-1)*NSCX*NSCY+(MSCY-1)*NSCX+MSCX+
     &                    NSCX*NSCY*(MC-1)
*--MSC is the new sub-cell number
                      IF (MSC.LT.1) MSC=1
                      IF (MSC.GT.MNSC) MSC=MNSC
                      IPL(NM)=MSC
                    ELSE
                      WRITE (*,*) 
     &' WARNING: EXCESS MOLECULE LIMIT - RESTART WITH AN INCREASED FNUM'
                    END IF
4                 CONTINUE
                END IF
5             CONTINUE
10          CONTINUE
20        CONTINUE
        END IF
100   CONTINUE
*--now the jet molecules
      IF (IJET.GT.0) THEN
        DO 150 L=1,MNSP
*--consider each species in turn
          VMP=SQRT(2.*BOLTZ*TMPJ/SP(5,L))
          A=BMEJ(L)+BMRJ(L)
          M=A
          BMRJ(L)=A-M
*--M molecules enter, remainder has been reset
          IF (M.GT.0) THEN
            IF (ABS(FVJ).GT.1.E-6) SC=FVJ/VMP
            FS1=SC+SQRT(SC*SC+2.)
            FS2=0.5*(1.+SC*(2.*SC-FS1))
* the above constants are required for the entering distn. of eqn (12.?)
            DO 110 K=1,M
              IF (NM.LT.MNM) THEN
                NM=NM+1
*--NM is now the number of the new molecule
                IF (ABS(FVJ).GT.1.E-6) THEN
                  QA=3.
                  IF (SC.LT.-3.) QA=ABS(SC)+1.
102               U=-QA+2.*QA*RND()
*--U is a potential normalised thermal velocity component
                  UN=U+SC
*--UN is a potential inward velocity component
                  IF (UN.LT.0.) GO TO 102
                  A=(2.*UN/FS1)*EXP(FS2-U*U)
                  IF (A.LT.RND()) GO TO 102
*--the inward normalised vel. component has been selected (eqn (12.?))
                  PV(2,NM)=UN*VMP
                ELSE
                  PV(2,NM)=SQRT(-LOG(RND()))*VMP
*--for a stationary external gas, use eqn (12.3)
                END IF
                CALL RVELC(PV(1,NM),PV(3,NM),VMP)
*--a single call of RVELC generates the two normal velocity components
                IF (ISPR(1,L).GT.0) CALL SROT(PR(NM),TMPJ,ISPR(1,L))
*
                PP(2,NM)=0.001*CH
*--the molecule is moved just off the boundary
                IPS(NM)=L
                R=SQRT((K-RND())/FLOAT(M))*RJ
                A=2.*PI*(K-RND())/FLOAT(M)
                IF (IJET.EQ.2) A=0.5*A
*--R and A are the random radius and azimuth relative to the jet centre
                PP(1,NM)=XJ+R*COS(A)
                PP(3,NM)=ZJ+R*SIN(A)
*
                MCX=(PP(1,NM)-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)
                MCZ=(PP(3,NM)-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+MCX
                MSCX=((PP(1,NM)-CG(1,MC))/CW)*(NSCX-.001)+1
                MSCZ=((PP(3,NM)-CG(5,MC))/CD)*(NSCZ-.001)+1
                MSC=(MSCZ-1)*NSCX*NSCY+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(NM)=MSC
              ELSE
                WRITE (*,*) 
     &' WARNING: EXCESS MOLECULE LIMIT - RESTART WITH AN INCREASED FNUM'
              END IF
110         CONTINUE
          END IF
150     CONTINUE
      END IF
 
      RETURN
      END
