*   INIT3.FOR
*
      SUBROUTINE INIT3
*
      include 'parameter.inc'
*
      DOUBLE PRECISION COL(MNSP,MNSP),MOVT,NCOL,SELT,SEPT,CS(7,MNC,MNSP)
      DOUBLE PRECISION CSR(MNC,MNSP)
 
*
*--set constants
*
      PI=3.141592654
      SPI=SQRT(PI)
      BOLTZ=1.380622E-23
*
*--set data variables to default values that they retain if the data
*----does not reset them to specific values
      FND=0.
      FTMP=273.
      VFX=0.
      VFY=0.
      VFZ=0.
      DO 100 N=1,6
        IB(N)=3
        DO 50 L=1,MNSP
          ISP(L)=1
          FSP(L)=0.
          BME(N,L)=0.
          BMR(N,L)=0.
50      CONTINUE
100   CONTINUE
      IJET=0
      DO 200 L=1,MNSP
        BMEJ(L)=0.
        BMRJ(L)=0.
200   CONTINUE

      CALL DATA3
*
*--set additional data on the gas
*
      IF (MNSP.EQ.1) ISPD=0
      DO 300 N=1,MNSP
        DO 250 M=1,MNSP
          IF ((ISPR(3,N).EQ.0).AND.(M.NE.N)) THEN
            SPR(1,N,M)=SPR(1,N,N)
            SPR(2,N,M)=SPR(2,N,N)
            SPR(3,N,M)=SPR(3,N,N)
          END IF
          IF ((ISPD.EQ.0).OR.(N.EQ.M)) THEN
            SPM(1,N,M)=0.25*PI*(SP(1,N)+SP(1,M))**2
*--the collision cross section is assumed to be given by eqn (1.35)
            SPM(2,N,M)=0.5*(SP(2,N)+SP(2,M))
            SPM(3,N,M)=0.5*(SP(3,N)+SP(3,M))
            SPM(4,N,M)=0.5*(SP(4,N)+SP(4,M))
*--mean values are used for ISPD=0
          ELSE
            SPM(1,N,M)=PI*SPM(1,N,M)**2
*--the cross-collision diameter is converted to the cross-section
          END IF
          SPM(5,N,M)=(SP(5,N)/(SP(5,N)+SP(5,M)))*SP(5,M)
*--the reduced mass is defined in eqn (2.7)
          SPM(6,N,M)=GAM(2.5-SPM(3,N,M))
250     CONTINUE
300   CONTINUE
*
*--initialise variables
*
      TIME=0.
      NM=0
      NPR=0
      NCOL=0
      MOVT=0.
      SELT=0.
      SEPT=0.
*
      DO 400 M=1,MNSP
        DO 350 N=1,MNSP
          COL(M,N)=0.
350     CONTINUE
400   CONTINUE
*
      FW=CB(2)-CB(1)
      FH=CB(4)-CB(3)
      FD=CB(6)-CB(5)
      CG(1,1)=CB(1)
      CW=FW/NCX
*--CW is the uniform cell width
      CG(3,1)=CB(3)
      CH=FH/NCY
*--CH is the uniform cell height
      CG(5,1)=CB(5)
      CD=FD/NCZ
*--CD is the uniform cell depth (in the z direction)
      CC=CW*CH*CD
*--CC is the uniform cell volume


        DO MC = 1, MNC
          CT(MC) = FTMP
        END DO

        DO mc = 1, MNC
          mx = ncordx(mc)
          if (mx .eq. 1) CG(1,mc) = CG(1,1)
        if (mx .gt. 1) CG(1,mc) = CG(2,mc-1)
        CG(2,mc) = CG(1,mc) + CW

        my = ncordy(mc)
        if (my .eq. 1) CG(3,mc) = CG(3,1)
        if (my.gt.1 .and. mx.eq.1) CG(3,mc)=CG(4,mc-1)
        if (my.gt.1 .and. mx.gt.1) CG(3,mc)=CG(3,mc-1)
        CG(4,mc) = CG(3,mc) + CH

        mz = ncordz(mc)
        if (mz .eq. 1) CG(5,mc)=CG(5,1)
        if (mz.gt.1 .and. mx.eq.1 .and. my.eq.1) CG(5,mc)=CG(6,mc-1)
        if (mz.gt.1 .and. (mx.gt.1 .or. my.gt.1)) CG(5,mc)=CG(5,mc-1)
        CG(6,mc)=CG(5,mc)+CD
      enddo

*--CCG is distributed
      do mc = 1, MNC
        do l = 1, MNSG
          do k = 1, MNSG
                CCG(2,mc,l,k)=0.
                CCG(1,mc,l,k)=SPM(1,1,1)*300.*SQRT(FTMP/300.)
          enddo
        enddo
      enddo


*
*--set sub-cells
*
      DO 600 N=1,MNC
        DO 550 J=1,NSCZ
          DO 520 M=1,NSCY
            DO 510 K=1,NSCX
              L=(N-1)*NSCX*NSCY*NSCZ+(J-1)*NSCX*NSCY+(M-1)*NSCX+K
              ISC(L)=N
510         CONTINUE
520       CONTINUE
550     CONTINUE
600   CONTINUE
*
      IF (IIS.GT.0.AND.ISG.GT.0) THEN
*--if IIS=1 generate initial gas with temperature FTMP
*
        DO 650 L=1,MNSP
          REM=0
          IF (IIS.EQ.1) VMP=SQRT(2.*BOLTZ*FTMP/SP(5,L))
*--VMP is the most probable speed in species L, see eqns (4.1) and (4.7)
          DO 620 N=1,MNC
*           A=FND*CC*FSP(L)/FNUM+REM
            A=FND*CC*FSP(L)/FNUM
*--A is the number of simulated molecules of species L in cell N to
*--simulate the required concentrations at a total number density of FND
            IF (N.LT.MNC) THEN
              MM=A
              REM=(A-MM)
*--the remainder REM is carried forward to the next cell
            ELSE
              MM=NINT(A)
            END IF
            IF (MM.GT.0) THEN
              DO 605 M=1,MM
                IF (NM.LT.MNM) THEN
*--round-off error could have taken NM to MNM+1
                  NM=NM+1
                  IF (MNSP.GT.1) IPS(NM)=L
                  PP(1,NM)=CG(1,N)+RND()*(CG(2,N)-CG(1,N))
                  NCOLM=(PP(1,NM)-CG(1,N))*(NSCX-.001)/CW+1
                  PP(2,NM)=CG(3,N)+RND()*(CG(4,N)-CG(3,N))
                  NROW=(PP(2,NM)-CG(3,N))*(NSCY-.001)/CH+1
                  PP(3,NM)=CG(5,N)+RND()*(CG(6,N)-CG(5,N))
                  NSLC=(PP(3,NM)-CG(5,N))*(NSCZ-.001)/CD+1
                  IPL(NM)=(N-1)*NSCX*NSCY*NSCZ+(NSLC-1)*NSCX*NSCY+
     &                    (NROW-1)*NSCX+NCOLM
*--species, position, and sub-cell number have been set
                  DO 602 K=1,3
                    CALL RVELC(PV(K,NM),A,VMP)
602               CONTINUE
                  PV(1,NM)=PV(1,NM)+VFX
                  PV(2,NM)=PV(2,NM)+VFY
                  PV(3,NM)=PV(3,NM)+VFZ
*--velocity components have been set
*--set the rotational energy
                  IF (ISPR(1,L).GT.0) CALL SROT(PR(NM),FTMP,ISPR(1,L))
                END IF
605           CONTINUE
            END IF
620       CONTINUE
650     CONTINUE
*
        WRITE (*,99001) NM
99001   FORMAT (' ',I6,' MOLECULES')
      END IF
      IF (IIS.GT.0) THEN
*
*--calculate the number of molecules that enter at each time step
*--across the six sides of the simulated region
        DO 700 N=1,6
          IF (IB(N).EQ.1) THEN
            WRITE (*,*) 'side',N
*--molecules enter from an external stream
            DO 660 L=1,MNSP
              VMP=SQRT(2.*BOLTZ*FTMP/SP(5,L))
*--VMP is the most probable speed in species L, see eqns (4.1) and (4.7)
              IF (N.EQ.1) SC=VFX/VMP
              IF (N.EQ.2) SC=-VFX/VMP
              IF (N.EQ.3) SC=VFY/VMP
              IF (N.EQ.4) SC=-VFY/VMP
              IF (N.EQ.5) SC=VFZ/VMP
              IF (N.EQ.6) SC=-VFZ/VMP
*--SC is the inward directed speed ratio
              IF (ABS(SC).LT.10.1) A=(EXP(-SC*SC)+SPI*SC*(1.+ERF(SC)))
     &                               /(2.*SPI)
              IF (SC.GT.10.) A=SC
              IF (SC.LT.-10.) A=0.
*--A is the non-dimensional flux of eqn (4.22)
              IF (N.EQ.1.OR.N.EQ.2) BME(N,L)=FND*FSP(L)
     &            *A*VMP*DTM*FH*FD/FNUM
              IF (N.EQ.3.OR.N.EQ.4) BME(N,L)=FND*FSP(L)
     &            *A*VMP*DTM*FW*FD/FNUM
              IF (N.EQ.5.OR.N.EQ.6) BME(N,L)=FND*FSP(L)
     &            *A*VMP*DTM*FH*FW/FNUM
              WRITE (*,*) ' species ',L,' entering mols ',BME(N,L)
*--the entering number is given by eqn (12.?)
660         CONTINUE
          END IF
700     CONTINUE
      END IF
*--now calculate the number that enter in jet
      IF (IJET.GT.0) THEN
*--molecules enter from a jet
        DO 750 L=1,MNSP
          VMP=SQRT(2.*BOLTZ*TMPJ/SP(5,L))
*--VMP is the most probable speed in species L, see eqns (4.1) and (4.7)
          SC=FVJ/VMP
*--SC is the inward directed speed ratio
          IF (ABS(SC).LT.10.1) A=(EXP(-SC*SC)+SPI*SC*(1.+ERF(SC)))
     &                           /(2.*SPI)
          IF (SC.GT.10.) A=SC
          IF (SC.LT.-10.) A=0.
*--A is the non-dimensional flux of eqn (4.22)
          AJ=PI*RJ*RJ
          IF (IJET.EQ.2) AJ=0.5*AJ
          BMEJ(L)=FNDJ*FSPJ(L)*A*VMP*DTM*AJ/FNUM
          WRITE (*,*) ' entering mols in jet ',BMEJ(L)
*--the entering number is given by eqn (12.?)
750     CONTINUE
      END IF
*
      RETURN
      END

*----------------------------------------------------------------------------

      function ncordx(mc)
        include 'parameter.inc'
        ncordx = mod(mc-1,NCX) + 1
      return
      end

      function ncordy(mc)
        include 'parameter.inc'
        ncordy = mod((mc-1) / NCX, NCY) + 1
      return
      end

      function ncordz(mc)
        include 'parameter.inc'
        ncordz = (mc-1) / (NCX*NCY) + 1
      return
      end
*------------------
