*   OUTT.FOR
 
 
*
*
      SUBROUTINE OUTT
*
*--output a TECPLOT DATA FILE to file DSMC2.TEC
*
      include 'parameter.inc'
      PARAMETER (NVARS=8,NVAR=8,NCOLL=26,NROW=19,NSLC=19)
      DIMENSION GQS(2,NVARS,NCOLL,NROW),GTS(2,NVARS,25,18)
      DIMENSION GQ(NVAR,NSLC,NROW),GT(NVAR,18,18)
*
      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
*
*  FIRST THE SURFACE PROPERTIES
*
      OPEN (4,FILE='DSMC3S.TEC',FORM='FORMATTED')
      WRITE (4,*) 'TITLE = "DSMC3 OUTPUT"'
      WRITE (4,*) 'VARIABLES = X, Y, Z, P, S1, S2, HI, HR '
*
 
      DO 100 KS=1,2
        IF (ISURF(KS).GT.0) THEN
          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)
*
          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))
                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 5 N=1,11
                SS(N)=0.
                DO 2 L=1,MNSP
                  SS(N)=SS(N)+CSS(N,K,L)
2               CONTINUE
5             CONTINUE
              DO 10 N=1,11
                SS(N)=SS(N)*A/AR
10            CONTINUE
              NRW=(K-NEL1+1)/(NCOLL-1)+1
              NCL=(K-NEL1+1)-(NRW-1)*(NCOLL-1)
              GTS(KS,1,NCL,NRW)=X
              GTS(KS,2,NCL,NRW)=Y
              GTS(KS,3,NCL,NRW)=Z
              GTS(KS,4,NCL,NRW)=SS(2)+SS(3)
              GTS(KS,5,NCL,NRW)=SS(4)+SS(9)
              GTS(KS,6,NCL,NRW)=SS(10)+SS(11)
              GTS(KS,7,NCL,NRW)=SS(5)
              GTS(KS,8,NCL,NRW)=SS(6)
            END IF
20        CONTINUE
        END IF
100   CONTINUE
*
      DO 300 KS=1,2
        DO 150 NY=1,NROW
          DO 120 NX=1,NCOLL
            IF (KS.EQ.1) THEN
              IF (NX.EQ.1) GQS(KS,1,NX,NY)=CB(1)
              IF (NX.NE.1) GQS(KS,1,NX,NY)=GQS(KS,1,NX-1,NY)+CW
              GQS(KS,2,NX,NY)=CB(3)
              IF (NY.EQ.1) GQS(KS,3,NX,NY)=CB(5)
              IF (NY.NE.1) GQS(KS,3,NX,NY)=GQS(KS,3,NX,NY-1)+CD
            ELSE
              IF (NX.EQ.1) GQS(KS,1,NX,NY)=CB(1)
              IF (NX.NE.1) GQS(KS,1,NX,NY)=GQS(KS,1,NX-1,NY)+CW
              GQS(KS,3,NX,NY)=CB(5)
              IF (NY.EQ.1) GQS(KS,2,NX,NY)=CB(3)
              IF (NY.NE.1) GQS(KS,2,NX,NY)=GQS(KS,2,NX,NY-1)+CH
            END IF
120       CONTINUE
150     CONTINUE
        DO 200 K=4,NVAR
          DO 180 NY=1,NROW
            DO 160 NX=1,NCOLL
              IF (NY.EQ.1) THEN
                IF (NX.EQ.1) GQS(KS,K,NX,NY)=GTS(KS,K,1,1)
                IF (NX.EQ.NCOLL) GQS(KS,K,NX,NY)=GTS(KS,K,25,1)
                IF (NX.GT.1.AND.NX.LT.NCOLL) GQS(KS,K,NX,NY)
     &              =0.5*(GTS(KS,K,NX-1,NY)+GTS(KS,K,NX,NY))
              END IF
              IF (NY.EQ.NROW) THEN
                IF (NX.EQ.1) GQS(KS,K,NX,NY)=GTS(KS,K,1,NROW-1)
                IF (NX.EQ.NCOLL) GQS(KS,K,NX,NY)=GTS(KS,K,25,NROW-1)
                IF (NX.GT.1.AND.NX.LT.NCOLL) GQS(KS,K,NX,NY)
     &              =0.5*(GTS(KS,K,NX-1,NROW-1)+GTS(KS,K,NX,NROW-1))
              END IF
              IF (NY.GT.1.AND.NY.LT.NROW) THEN
                IF (NX.EQ.1) GQS(KS,K,NX,NY)
     &              =0.5*(GTS(KS,K,1,NY-1)+GTS(KS,K,1,NY))
                IF (NX.EQ.NCOLL) GQS(KS,K,NX,NY)
     &              =0.5*(GTS(KS,K,25,NY-1)+GTS(KS,K,25,NY))
                IF (NX.GT.1.AND.NX.LT.NCOLL) GQS(KS,K,NX,NY)
     &              =0.25*(GTS(KS,K,NX-1,NY-1)+GTS(KS,K,NX,NY-1)
     &              +GTS(KS,K,NX-1,NY)+GTS(KS,K,NX,NY))
              END IF
160         CONTINUE
180       CONTINUE
200     CONTINUE
        WRITE (4,*) 'ZONE T="DSMC3", I=',NSLC,' , J=',NROW,', F=POINT'
        DO 250 NY=1,NROW
          DO 220 NX=1,NCOLL
            WRITE (4,99001) (GQS(KS,K,NX,NY),K=1,NVAR)
220       CONTINUE
250     CONTINUE
99001   FORMAT (11E14.5)
300   CONTINUE
*
      CLOSE (4)
*
* generate flowfield sections at x=5, 10, 15, 20 cells from leading edge
*
      OPEN (4,FILE='DSMC3.TEC',FORM='FORMATTED')
      WRITE (4,*) 'TITLE = "DSMC3 OUTPUT"'
      WRITE (4,*) 'VARIABLES = Z, Y, N, T, U, V, W, M '
*
      DO 600 KS=1,4
 
        DO 350 NRW=1,18
          DO 340 NSL=1,18
            IF (KS.EQ.1) NCL=10
            IF (KS.EQ.2) NCL=15
            IF (KS.EQ.3) NCL=20
            IF (KS.EQ.4) NCL=25
 
            N=(NSL-1)*NCX*NCY+(NRW-1)*NCX+NCL
 
            A=FNUM/(CC*NSMP)
            SN=0.
            SM=0.
            DO 310 K=1,3
              SMU(K)=0.
310         CONTINUE
            SMCC=0.
            SRE=0.
            SRDF=0.
            DO 320 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 315 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
315           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
320         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 330 K=1,3
              VEL(K)=SMU(K)/SM
              SVEL(K,N)=VEL(K)
330         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)
 
 
 
            XC=0.5*(CG(4,N)+CG(5,N))
            YC=0.5*(CG(3,N)+CG(4,N))
*--XC,YC are the x,y coordinates of the midpoint of the cell
            GT(1,NSL,NRW)=XC
            GT(2,NSL,NRW)=YC
            GT(3,NSL,NRW)=DENN
            GT(4,NSL,NRW)=TEMP
            GT(5,NSL,NRW)=VEL(1)
            GT(6,NSL,NRW)=VEL(2)
            GT(7,NSL,NRW)=VEL(3)
            IF (TEMP.GT.1.) THEN
              GT(8,NSL,NRW)=SQRT(VEL(1)**2+VEL(2)**2+VEL(3)**2)
     &                      /SQRT(1.6667*BOLTZ*TEMP/SP(5,1))
            ELSE
              GT(8,NSL,NRW)=0.
            END IF
*  ABOVE ASSUMES SIMPLE MONTOMIC GAS
340       CONTINUE
350     CONTINUE
*
        DO 400 NY=1,NROW
          DO 360 NX=1,NSLC
            IF (NX.EQ.1) GQ(1,NX,NY)=CB(5)
            IF (NX.NE.1) GQ(1,NX,NY)=(NX-1)*CD
            IF (NY.EQ.1) GQ(2,NX,NY)=CB(3)
            IF (NY.NE.1) GQ(2,NX,NY)=(NY-1)*CH
360       CONTINUE
400     CONTINUE
        DO 450 K=3,NVAR
          DO 420 NY=1,NROW
            DO 410 NX=1,NSLC
              IF (NY.EQ.1) THEN
                IF (NX.EQ.1) GQ(K,NX,NY)=GT(K,1,1)
                IF (NX.EQ.NSLC) GQ(K,NX,NY)=GT(K,18,1)
                IF (NX.GT.1.AND.NX.LT.NSLC) GQ(K,NX,NY)
     &              =0.5*(GT(K,NX-1,NY)+GT(K,NX,NY))
              END IF
              IF (NY.EQ.NROW) THEN
                IF (NX.EQ.1) GQ(K,NX,NY)=GT(K,1,NROW-1)
                IF (NX.EQ.NSLC) GQ(K,NX,NY)=GT(K,18,NROW-1)
                IF (NX.GT.1.AND.NX.LT.NSLC) GQ(K,NX,NY)
     &              =0.5*(GT(K,NX-1,NROW-1)+GT(K,NX,NROW-1))
              END IF
              IF (NY.GT.1.AND.NY.LT.NROW) THEN
                IF (NX.EQ.1) GQ(K,NX,NY)=0.5*(GT(K,1,NY-1)+GT(K,1,NY))
                IF (NX.EQ.NSLC) GQ(K,NX,NY)
     &              =0.5*(GT(K,18,NY-1)+GT(K,18,NY))
                IF (NX.GT.1.AND.NX.LT.NSLC) GQ(K,NX,NY)
     &              =0.25*(GT(K,NX-1,NY-1)+GT(K,NX,NY-1)+GT(K,NX-1,NY)
     &              +GT(K,NX,NY))
              END IF
410         CONTINUE
420       CONTINUE
450     CONTINUE
        WRITE (4,*) 'ZONE T="DSMC2", I=',NSLC,' , J=',NROW,', F=POINT'
        DO 500 NY=1,NROW
          DO 460 NX=1,NSLC
            WRITE (4,99002) (GQ(K,NX,NY),K=1,NVAR)
460       CONTINUE
500     CONTINUE
99002   FORMAT (11E14.5)
 
600   CONTINUE
*
      CLOSE (4)
*
      RETURN
      END
