*   OUT3.FOR
*
#if (ALPHA | RS6K)
#define IFIX    IDINT
#endif
*
      SUBROUTINE OUT3
*
*--output a progressive set of results to file DSMC3.OUT.
*
      include 'parameter.inc'
*
      DOUBLE PRECISION COL(MNSP,MNSP),MOVT,NCOL,SELT,SEPT,CS(7,MNC,MNSP)
      DOUBLE PRECISION CSR(MNC,MNSP)
      DOUBLE PRECISION CSS(11,MNSE,MNSP)
*
      DOUBLE PRECISION VEL(3),SMU(3),SVEL(3,MNC),SN,SM,SMCC,SRDF,SRE,TT,
     &                 TROT,DBOLTZ,SS(11)
      DBOLTZ=BOLTZ
*
      OPEN (4,FILE='DSMC3.OUT',FORM='FORMATTED')
*
      WRITE (4,*) ' FLOW SAMPLED FROM TIME ',TIMI,' TO TIME ',TIME
      WRITE (4,*) ' COLLISIONS:-'
      WRITE (4,99001) ((IFIX(COL(M,L)),M=1,MNSP),L=1,MNSP)
99001 FORMAT (5I12)
      WRITE (4,*) ' TOTAL NUMBER OF SAMPLES ',NSMP
      WRITE (4,*) NM,' MOLECULES'
      WRITE (4,*) MOVT,' TOTAL MOLECULAR MOVES'
      IF (NCOL.GT.0) THEN
        WRITE (4,*) INT(SELT),' SELECTIONS ',INT(NCOL),
     &              ' COLLISION EVENTS, RATIO  ',REAL(NCOL/SELT)
        IF (NCOL.GT.0) WRITE (4,*) ' MEAN COLLISION SEPARATION ',
     &                             REAL(SEPT/NCOL)
      END IF
*
      WRITE (4,*)
      DO 100 KS=1,2
        IF (ISURF(KS).GT.0) THEN
          WRITE (4,*) ' SURFACE ',KS
          WRITE (4,*)
          NEL=(LIMS(KS,3)-LIMS(KS,2)+1)*(LIMS(KS,5)-LIMS(KS,4)+1)
          IF (KS.EQ.1) THEN
            NEL1=1
            NEL2=NEL
          ELSE
            NEL1=(LIMS(1,3)-LIMS(1,2)+1)*(LIMS(1,5)-LIMS(1,4)+1)+1
            NEL2=NEL1+NEL-1
          END IF
          A=FNUM/(TIME-TIMI)
          WRITE (4,*) 
     &'    X COORD     Y COORD      Z COORD     SAMPLE     FRACTIONSPECI
     &ES 1,    SPECIES2....'
          DO 20 K=NEL1,NEL2
            NC=ICS(K)
            IF (ICS(K).GT.0) THEN
              IF (ISURF(KS).LT.3) THEN
                Y=(LIMS(KS,1)-1)*CH
                X=0.5*(CG(1,NC)+CG(2,NC))
                Z=0.5*(CG(5,NC)+CG(6,NC))
              END IF
              IF (ISURF(KS).GT.2.AND.ISURF(KS).LT.5) THEN
                X=(LIMS(KS,1)-1)*CW
                Y=0.5*(CG(3,NC)+CG(4,NC))
                Z=0.5*(CG(5,NC)+CG(6,NC))
              END IF
              IF (ISURF(KS).GT.4) THEN
                Z=(LIMS(KS,1)-1)*CD
                Y=0.5*(CG(3,NC)+CG(4,NC))
                X=0.5*(CG(1,NC)+CG(2,NC))
              END IF
              SS(1)=0.
              DO 5 L=1,MNSP
                SS(1)=SS(1)+CSS(1,K,L)
5             CONTINUE
              WRITE (4,99002) X,Y,Z,SS(1),(CSS(1,K,L)/SS(1),L=1,MNSP)
99002         FORMAT (3F12.5,F12.1,7F13.6)
            END IF
20        CONTINUE
*
          WRITE (4,*) 
     &'     X COORD      Y COORD     Z COORD     NUM FLUX    INC PRESSRE
     &FL PRESS   INC SH ST1    REFL SH ST1   INC SH ST2   REFL SH ST2   
     &INC TR EN   REF TR EN   INC ROT EN   REFL ROTEN  NET HEAT FLUX'
          DO 40 K=NEL1,NEL2
            NC=ICS(K)
            IF (ICS(K).GT.0) THEN
              IF (ISURF(KS).LT.3) THEN
                Y=(LIMS(KS,1)-1)*CH
                X=0.5*(CG(1,NC)+CG(2,NC))
                Z=0.5*(CG(5,NC)+CG(6,NC))
                AR=CW*CD
              END IF
              IF (ISURF(KS).GT.2.AND.ISURF(KS).LT.5) THEN
                X=(LIMS(KS,1)-1)*CW
                Y=0.5*(CG(3,NC)+CG(4,NC))
                Z=0.5*(CG(5,NC)+CG(6,NC))
                AR=CH*CD
              END IF
              IF (ISURF(KS).GT.4) THEN
                Z=(LIMS(KS,1)-1)*CD
                Y=0.5*(CG(3,NC)+CG(4,NC))
                X=0.5*(CG(1,NC)+CG(2,NC))
                AR=CW*CH
              END IF
              DO 25 N=1,11
                SS(N)=0.
                DO 22 L=1,MNSP
                  SS(N)=SS(N)+CSS(N,K,L)
22              CONTINUE
25            CONTINUE
              DO 30 N=1,11
                SS(N)=SS(N)*A/AR
30            CONTINUE
              WRITE (4,99003) X,Y,Z,(SS(N),N=1,4),SS(9),SS(10),SS(11),
     &                        SS(5),SS(6),SS(7),SS(8),SS(5)+SS(6)+SS(7)
     &                        +SS(8)
99003         FORMAT (15E13.5)
            END IF
40        CONTINUE
        END IF
100   CONTINUE
*
      WRITE (4,*) ' FLOWFIELD PROPERTIES '
      WRITE (4,*) 'SAMPLES'
      WRITE (4,*) ' CELL     N SP 1    N SP 2     ETC '
      DO 200 N=1,MNC
        WRITE (4,99004) N,(IFIX(CS(1,N,L)),L=1,MNSP)
200   CONTINUE
99004 FORMAT (' ',6I9)
*
      WRITE (4,*) ' FLOWFIELD PROPERTIES'
      WRITE (4,*) 
     &'  CELL   X COORD   Y COORD   Z COORD   DENSITY   TR TEMP   ROTL T
     &EMP    OV. TEMP      U         V         W'
*--first the mixture properties
      DO 400 N=1,MNC
        A=FNUM/(CC*NSMP)
        SN=0.
        SM=0.
        DO 250 K=1,3
          SMU(K)=0.
250     CONTINUE
        SMCC=0.
        SRE=0.
        SRDF=0.
        DO 300 L=1,MNSP
          SN=SN+CS(1,N,L)
*--SN is the number sum
          SM=SM+SP(5,L)*CS(1,N,L)
*--SM is the sum of molecular masses
          DO 260 K=1,3
            SMU(K)=SMU(K)+SP(5,L)*CS(K+1,N,L)
*--SMU(1 to 3) are the sum of mu, mv, mw
260       CONTINUE
          SMCC=SMCC+(CS(5,N,L)+CS(6,N,L)+CS(7,N,L))*SP(5,L)
*--SMCC is the sum of m(u**2+v**2+w**2)
          SRE=SRE+CSR(N,L)
*--SRE is the sum of rotational energy
          SRDF=SRDF+ISPR(1,L)*CS(1,N,L)
*--SRDF is the sum of the rotational degrees of freedom
          SUU=SUU+SP(5,L)*CS(5,N,L)
*--SUU is the sum of m*u*u
300     CONTINUE
        DENN=SN*A
*--DENN is the number density, see eqn (1.34)
        DEN=DENN*SM/SN
*--DEN is the density, see eqn (1.42)
        DO 350 K=1,3
          VEL(K)=SMU(K)/SM
          SVEL(K,N)=VEL(K)
350     CONTINUE
*--VEL and SVEL are the stream velocity components, see eqn (1.43)
        UU=VEL(1)**2+VEL(2)**2+VEL(3)**2
        TT=(SMCC-SM*UU)/(3.D00*DBOLTZ*SN)
*--TT is the translational temperature, see eqn (1.51)
        IF (SRDF.GT.1.E-6) TROT=(2.D00/DBOLTZ)*SRE/SRDF
*--TROT is the rotational temperature, see eqn (11.11)
        TEMP=(3.D00*TT+(SRDF/SN)*TROT)/(3.+SRDF/SN)
*--TEMP is the overall temperature, see eqn (11.12)
        CT(N)=TEMP
        XC=0.5*(CG(1,N)+CG(2,N))
        YC=0.5*(CG(3,N)+CG(4,N))
        ZC=0.5*(CG(5,N)+CG(6,N))
*--XC,YC are the x,y coordinates of the midpoint of the cell
        WRITE (4,99005) N,XC,YC,ZC,DEN,TT,TROT,TEMP,VEL(1),VEL(2),VEL(3)
99005   FORMAT (' ',I5,3F11.5,1P,E12.4,0P,6F10.4,2E12.4)
400   CONTINUE
*
      WRITE (4,*)
      DO 500 L=1,MNSP
*--now the properties of the separate species
        WRITE (4,*) ' SPECIES ',L
        WRITE (4,*) 
     &' CELL   X COORD    Y COORD    Z COORD   N DENS    DENSITY    TTX 
     &    TTY       TTZ     TR TEMP    ROT TEMP     TEMP    U DIF VEL V 
     &DIF VEL   W DIF VEL '
        DO 450 N=1,MNC
          A=FNUM/(CC*NSMP)
          DENN=CS(1,N,L)*A
*--DENN is the partial number density
          DEN=SP(5,L)*DENN
*--DEN is the partial density, see eqn (1.13)
          DO 420 K=1,3
            VEL(K)=CS(K+1,N,L)/CS(1,N,L)
*--VEL defines the average velocity of the species L molecules
420       CONTINUE
          UU=VEL(1)**2+VEL(2)**2+VEL(3)**2
          TTX=(SP(5,L)/DBOLTZ)*(CS(5,N,L)/CS(1,N,L)-VEL(1)**2)
          TTY=(SP(5,L)/DBOLTZ)*(CS(6,N,L)/CS(1,N,L)-VEL(2)**2)
          TTZ=(SP(5,L)/DBOLTZ)*(CS(7,N,L)/CS(1,N,L)-VEL(3)**2)
*--the component temperatures are based on eqn (1.30)
          TT=(SP(5,L)/(3.D00*DBOLTZ))
     &       *((CS(5,N,L)+CS(6,N,L)+CS(7,N,L))/CS(1,N,L)-UU)
*--TT is the translational temperature, see eqn (1.29)
          IF (ISPR(1,L).GT.0) THEN
            TROT=2.D00*CSR(N,L)/(ISPR(1,L)*DBOLTZ*CS(1,N,L))
          ELSE
            TROT=0.
          END IF
*--TROT is the rotational temperature, see eqn (11.10)
          TEMP=(3.D00*TT+ISPR(1,L)*TROT)/(3.+ISPR(1,L))
          DO 440 K=1,3
            VEL(K)=VEL(K)-SVEL(K,N)
*--VEL now defines the diffusion velocity of species L, see eqn (1,45)
440       CONTINUE
          XC=0.5*(CG(1,N)+CG(2,N))
          YC=0.5*(CG(3,N)+CG(4,N))
          ZC=0.5*(CG(5,N)+CG(6,N))
          WRITE (4,99006) N,XC,YC,ZC,DENN,DEN,TTX,TTY,TTZ,TT,TROT,TEMP,
     &                    VEL(1),VEL(2),VEL(3)
99006     FORMAT (' ',I5,3F10.5,1P,2E12.4,0P,9F10.4)
450     CONTINUE
500   CONTINUE
*
*
      CLOSE (4)
*
      RETURN
      END
