C	WESTERN MICHIGAN UNIVERSITY
C	BANOVA.F4 (FILENAME ON LIBRARY DECTAPE)
C	BANOVA, 1.9.1 (CALLING NAME, SUBLST #)
C	BALANCED ANALYSIS OF VARIANCE
C	ADAPTED FROM DECUS(NO. 10-101, IBM-SSP)
C	MODIFIED BY SAM ANEMA AND R.R. BARR III
C	LIBRARY DECTAPE PROGRAMS USED:  USAGE.MAC
C	FORWMU PROGS USED:  TTYPTY, PRINTS, DEVCHG, EXISTS
C	APLIB PROGS. USED:  FISHER, GETFOR, IO
C	INTERNAL SUBR. USED:  AVCAL, AVDAT, MEANQ, SETUP
C	ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
      DIMENSION F(63)
      DIMENSION ID(16),LFMT(96)
C---------------SEE COMMENTS FOR SUBR. AVCAL RELATING TO X(3000).
       DIMENSION X(3000)
      DIMENSION HEAD(6),LEVEL(6),ISTEP(6),KOUNT(6),LASTS(6)
      DIMENSION SUMSQ(63),NDF(63),SMEAN(63)
      DIMENSION FMT(15)
1     FORMAT(A4,A2,I2,A4,3X,11(A1,I4)/(A1,I4,A1,I4,A1,
     1I4,A1,I4,A1,I4))
2     FORMAT(' ','ANALYSIS OF VARIANCE............',A4,A2//)
3     FORMAT(' ','LEVELS OF FACTORS'/(3X,A1,7X,I4))
4     FORMAT(//' ','GRAND MEAN',F20.5////)
5     FORMAT(' ','SOURCE OF',18X,'SUMS OF',10X,'DEGREES OF',9X,
     1'MEAN'/' VARIATION',18X,'SQUARES',11X,'FREEDOM',10X,'SQUARES'/)
6     FORMAT(' ',12A1,F20.5,10X,I6,F20.5)
7     FORMAT(' ','TOTAL',7X,F20.5,10X,I6)
8     FORMAT(12F6.0)
      IOUT=3
      INP=2
      NOUTD=-1
      IND=-4
      WRITE(NOUTD,203)
203   FORMAT(//20X,'WMU'/5X,'BALANCED ANALYSIS OF VARIANCE'/)
C	CALL USAGE('BANOVA')
C---------------TTYPTY RETURNS ZERO - TTY JOB, MINUS ONE - BATCH JOB
      CALL TTYPTY(ICODE)
C---------------NDV1, NDV2 ARE RETURNED.  OTHER ARGS ARE INPUT.  1 CAUSES
C--------------- OUTPUT? TO PRINT, 0 CAUSES INPUT? TO PRINT.  ICODE
C--------------- COMES FROM CALL TTYPTY.
      CALL IO(IOUT,NDV1,NOUTD,IND,1,ICODE)
201   CALL IO(INP,NDV2,NOUTD,IND,0,ICODE)
C---------------LFMT, ISTD ARE RETURNED.  OTHER ARGS ARE INPUT.  96=NO.
C--------------- OF FMT. WORDS FOR OBJ. TIME FORMAT (6 LINES)
C--------------- 2 MEANS F - TYPE ONLY FOR FORMAT
      CALL GETFOR(NOUTD,IND,LFMT,ISTD,96,2)
C---------------USER CHOSE STANDARD FORMAT IF ISTD=1
      IF(ISTD.EQ.1)LFMT(1)='(10F)'
      BLANK='    '
      PR='    '
      PR1='  '
      WRITE(NOUTD,220)
220   FORMAT(' ','ENTER NUMBER OF FACTORS.'/)
      READ(IND,230)K
230   FORMAT(I)
      WRITE(NOUTD,240)
240   FORMAT(' ','ENTER NAMES OF FACTORS. '
     1,'(SINGLE ALPHAMERIC CHARACTERS)'/)
      READ(IND,250)(HEAD(I),I=1,K)
250   FORMAT(10(A1,1X))
      WRITE(NOUTD,260)
260   FORMAT(' ','ENTER NUMBER OF LEVELS FOR EACH FACTOR.'/)
      READ(IND,270)(LEVEL(I),I=1,K)
270   FORMAT(10I)
      N=LEVEL(1)
      DO 280 I=2,K
280   N=N*LEVEL(I)
      LINES=IFIX(FLOAT(N)/10.0+.999)
1002  WRITE(NOUTD,1003)
1003  FORMAT(' ','ENTER IDENTIFICATION'/)
      READ(IND,1004)(ID(I),I=1,16)
1004  FORMAT(16A5)
      IDC=INP
      IF(NDV2.NE.'TTY')GO TO 739
      IF(ISTD.EQ.1)WRITE(NOUTD,740)
740   FORMAT(' ENTER DATA (AT MOST 10 PER LINE)'/)
      IF(ISTD.NE.1)WRITE(NOUTD,330)
      CALL SETUP(LINES,NOUTD,IND)
      IDC=1
      CALL IFILE(1,'JOB')
      GO TO 337
739   WRITE(NOUTD,741)
741   FORMAT(' YOUR DATA IS BEING READ.'/)
330   FORMAT(' ','ENTER DATA.'/)
337   READ(IDC,LFMT)(X(I),I=1,N)
210   WRITE(IOUT,2)PR,PR1
      WRITE(IOUT,1008)(ID(I),I=1,16)
1008  FORMAT(' ',16A5)
      WRITE(IOUT,3)(HEAD(I),LEVEL(I),I=1,K)
      WRITE(IOUT,401)
401   FORMAT(///' MEANS ON LEVELS'//' FACTOR',5X,'LEVEL',9X,'MEAN'/)
      KA=1
      KF=N
      DO 490 I=1,K
      IF(I.EQ.1)GO TO 491
      KA=KA*LEVEL(I-1)
491   KF=KF/LEVEL(I)
      XNUM=FLOAT(N/LEVEL(I))
      DO 480 J=1,LEVEL(I)
      SUM=0.0
      DO 481 KK=1,KF
      KSTART=(KK-1)*LEVEL(I)*KA+(J-1)*KA+1
      KLAST=KSTART+KA-1
      DO 482 KL=KSTART,KLAST
482   SUM=SUM+X(KL)
481   CONTINUE
      SUM=SUM/XNUM
      WRITE(IOUT,483)HEAD(I),J,SUM
483   FORMAT(4X,A1,6X,I4,F17.4)
480   CONTINUE
490   CONTINUE
      CALL AVDAT (K,LEVEL,N,X,L,ISTEP,KOUNT)
      CALL AVCAL (K,LEVEL,X,L,ISTEP,LASTS)
      CALL MEANQ (K,LEVEL,X,GMEAN,SUMSQ,NDF,SMEAN,ISTEP,KOUNT,LASTS)
      WRITE(IOUT,4)GMEAN
      WRITE(IOUT,5)
      LL=(2**K)-1
      ISTEP(1)=1
      DO 105 I=2,K
105   ISTEP(I)=0
      DO 110 I=1,15
110   FMT(I)=BLANK
      NN=0
      SUM=0.
120   NN=NN+1
      L=0
      DO 140 I=1,K
      FMT(I)=BLANK
      IF(ISTEP(I))130,140,130
130   L=L+1
      FMT(L)=HEAD(I)
140   CONTINUE
      WRITE(IOUT,6)(FMT(I),I=1,12),SUMSQ(NN),NDF(NN),SMEAN(NN)
      SUM=SUM+SUMSQ(NN)
      IF(NN-LL)145,170,170
145   DO 160 I=1,K
      IF(ISTEP(I))147,150,147
147   ISTEP(I)=0
      GO TO 160
150   ISTEP(I)=1
      GO TO 120
160   CONTINUE
170   N=N-1
      WRITE(IOUT,7)SUM,N
      WRITE(NOUTD,1009)
1009   FORMAT(///' WOULD YOU LIKE TO POOL OVER THE LAST FACTOR?'/)
      READ(IND,1010)IP
      IF(IP.NE.'YES')GO TO 201
1010   FORMAT(A3)
19    NST=2**(K-1)
      ISTEP(1)=1
      DO 40 I=2,K
40    ISTEP(I)=0
      ERROR=0.
      NDFT=0
      NN=0
      DO 20 I=NST,LL
      ERROR=ERROR+SUMSQ(I)
20    NDFT=NDFT+NDF(I)
      SMEAN(NST)=ERROR/FLOAT(NDFT)
      SUMSQ(NST)=ERROR
      DO 27 I=1,NST-1
27    F(I)=SMEAN(I)/SMEAN(NST)
      NDF(NST)=NDFT
      WRITE(IOUT,51)
28    NN=NN+1
      L=0
      DO 23 I=1,K
      FMT(I)=BLANK
      IF(ISTEP(I))22,23,22
22    L=L+1
      FMT(L)=HEAD(I)
23    CONTINUE
	FPRB=FISHER(NDF(NN),NDF(NST),F(NN))
      WRITE(IOUT,29)(FMT(I),I=1,8),SUMSQ(NN),NDF(NN),SMEAN(NN),
     1F(NN),FPRB
      IF(NN-(NST-1))31,32,32
31    DO 33 I=1,K
      IF(ISTEP(I))34,35,34
34    ISTEP(I)=0
      GO TO 33
35    ISTEP(I)=1
      GO TO 28
33    CONTINUE
32    WRITE(IOUT,38)SUMSQ(NST),NDF(NST),SMEAN(NST)
      WRITE(IOUT,39)SUM,N
      WRITE(NOUTD,204)
204   FORMAT(///)
      GO TO 201
51    FORMAT(////' ANALYSIS OF VARIANCE WITH LAST FACTOR',
     1' POOLED AND CALLED ERROR'//' SOURCE OF',7X,'SUMS OF',
     13X,'DEGREES OF',6X,'MEAN',9X,'F',9X,'F'/' VARIATION',7X,'SQUARES',
     15X,'FREEDOM',6X,'SQUARES',4X,'VALUES',5X,'PROB'/)
29    FORMAT(' ',8A1,F17.5,I9,F15.5,F10.3,F10.5)
38    FORMAT(' ERROR   ',F17.5,I9,F15.5)
39    FORMAT(' TOTAL   ',F17.5,I9)

      END
C
C     ..................................................................
C
C        SUBROUTINE AVCAL
C
C        PURPOSE
C           PERFORM THE CALCULUS OF A FACTORIAL EXPERIMENT USING
C           OPERATOR SIGMA AND OPERATOR DELTA.  THIS SUBROUTINE IS
C           PRECEDED BY SUBROUTINE ADVAT AND FOLLOWED BY SUBROUTINE
C           MEANQ IN THE PERFORMANCE OF ANALYSIS OF VARIANCE FOR A
C           COMPLETE FACTORIAL DESIGN.
C
C        USAGE
C           CALL AVCAL (K,LEVEL,X,L,ISTEP,LASTS)
C
C        DESCRIPTION OF PARAMETERS
C           K     - NUMBER OF VARIABLES (FACTORS). K MUST BE .GT. ONE.
C           LEVEL - INPUT VECTOR OF LENGTH K CONTAINING LEVELS (CATE-
C                   GORIES) WITHIN EACH VARIABLE.
C           X     - INPUT VECTOR CONTAINING DATA.  DATA HAVE BEEN PLACED
C                   IN VECTOR X BY SUBROUTINE AVDAT.  THE LENGTH OF X
C                   IS (LEVEL(1)+1)*(LEVEL(2)+1)*...*(LEVEL(K)+1).
C           L     - THE POSITION IN VECTOR X WHERE THE LAST INPUT DATA
C                   IS LOCATED.  L HAS BEEN CALCULATED BY SUBROUTINE
C                   AVDAT.
C           ISTEP - INPUT VECTOR OF LENGTH K CONTAINING STORAGE CONTROL
C                   STEPS WHICH HAVE BEEN CALCULATED BY SUBROUTINE
C                   AVDAT.
C           LASTS - WORKING VECTOR OF LENGTH K.
C
C        REMARKS
C           THIS SUBROUTINE MUST FOLLOW SUBROUTINE AVDAT.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           THE METHOD IS BASED ON THE TECHNIQUE DISCUSSED BY H. O.
C           HARTLEY IN 'MATHEMATICAL METHODS FOR DIGITAL COMPUTERS',
C           EDITED BY A. RALSTON AND H. WILF, JOHN WILEY AND SONS,
C           1962, CHAPTER 20.
C
C     ..................................................................
C
C---------------K, LEVEL, X, L, ISTEP ARE INPUT.  LASTS IS
C--------------- OUTPUT.  X IS MODIFIED.
      SUBROUTINE AVCAL (K,LEVEL,X,L,ISTEP,LASTS)
      DIMENSION LEVEL(1),X(1),ISTEP(1),LASTS(1)
C
C        ...............................................................
C
C        IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
C        C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
C        STATEMENT WHICH FOLLOWS.
C
C     DOUBLE PRECISION X,SUM
C
C        THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
C        APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
C        ROUTINE.
C
C        ...............................................................
C
C     CALCULATE THE LAST DATA POSITION OF EACH FACTOR
C
      LASTS(1)=L+1
      DO 145 I=2,K
  145 LASTS(I)=LASTS(I-1)+ISTEP(I)
C
C     PERFORM CALCULUS OF OPERATION
C
  150 DO 175 I=1,K
      L=1
      LL=1
      SUM=0.0
      NN=LEVEL(I)
      FN=NN
      INCRE=ISTEP(I)
      LAST=LASTS(I)
C
C     SIGMA OPERATION
C
  155 DO 160 J=1,NN
      SUM=SUM+X(L)
  160 L=L+INCRE
      X(L)=SUM
C
C     DELTA OPERATION
C
      DO 165 J=1,NN
      X(LL)=FN*X(LL)-SUM
  165 LL=LL+INCRE
      SUM=0.0
      IF(L-LAST) 167, 175, 175
  167 IF(L-LAST+INCRE) 168, 168, 170
  168 L=L+INCRE
      LL=LL+INCRE
      GO TO 155
  170 L=L+INCRE+1-LAST
      LL=LL+INCRE+1-LAST
      GO TO 155
  175 CONTINUE
      RETURN
      END
C
C     ..................................................................
C
C        SUBROUTINE AVDAT
C
C        PURPOSE
C           PLACE DATA FOR ANALYSIS OF VARIANCE IN PROPERLY DISTRIBUTED
C           POSITIONS OF STORAGE.  THIS SUBROUTINE IS NORMALLY FOLLOWED
C           BY CALLS TO AVCAL AND MEANQ SUBROUTINES IN THE PERFORMANCE
C           OF ANALYSIS OF VARIANCE FOR A COMPLETE FACTORIAL DESIGN.
C
C        USAGE
C           CALL AVDAT (K,LEVEL,N,X,L,ISTEP,KOUNT)
C
C        DESCRIPTION OF PARAMETERS
C           K     - NUMBER OF VARIABLES (FACTORS). K MUST BE .GT. ONE.
C           LEVEL - INPUT VECTOR OF LENGTH K CONTAINING LEVELS (CATE-
C                   GORIES) WITHIN EACH VARIABLE.
C           N     - TOTAL NUMBER OF DATA POINTS READ IN.
C           X     - WHEN THE SUBROUTINE IS CALLED, THIS VECTOR CONTAINS
C                   DATA IN LOCATIONS X(1) THROUGH X(N).  UPON RETURNING
C                   TO THE CALLING ROUTINE, THE VECTOR CONTAINS THE DATA
C                   IN PROPERLY REDISTRIBUTED LOCATIONS OF VECTOR X.
C                   THE LENGTH OF VECTOR X IS CALCULATED BY (1) ADDING
C                   ONE TO EACH LEVEL OF VARIABLE AND (2) OBTAINING THE
C                   CUMULATIVE PRODUCT OF ALL LEVELS.  (THE LENGTH OF
C                   X = (LEVEL(1)+1)*(LEVEL(2)+1)*...*(LEVEL(K)+1).)
C           L     - OUTPUT VARIABLE CONTAINING THE POSITION IN VECTOR X
C                   WHERE THE LAST INPUT DATA IS STORED.
C           ISTEP - OUTPUT VECTOR OF LENGTH K CONTAINING CONTROL STEPS
C                   WHICH ARE USED TO LOCATE DATA IN PROPER POSITIONS
C                   OF VECTOR X.
C           KOUNT - WORKING VECTOR OF LENGTH K.
C        REMARKS
C           INPUT DATA MUST BE ARRANGED IN THE FOLLOWING MANNER.
C           CONSIDER THE 3-VARIABLE ANALYSIS OF VARIANCE DESIGN, WHERE
C           ONE VARIABLE HAS 3 LEVELS AND THE OTHER TWO VARIABLES HAVE
C           2 LEVELS.  THE DATA MAY BE REPRESENTED IN THE FORM X(I,J,K),
C           I=1,2,3  J=1,2  K=1,2.  IN ARRANGING DATA, THE INNER
C           SUBSCRIPT, NAMELY I, CHANGES FIRST.  WHEN I=3, THE NEXT
C           INNER SUBSCRIPT, J, CHANGES AND SO ON UNTIL I=3, J=2, AND
C           K=2.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           THE METHOD IS BASED ON THE TECHNIQUE DISCUSSED BY H. O.
C           HARTLEY IN 'MATHEMATICAL METHODS FOR DIGITAL COMPUTERS',
C           EDITED BY A. RALSTON AND H. WILF, JOHN WILEY AND SONS,
C           1962, CHAPTER 20.
C
C     ..................................................................
C
C---------------K, N, X, ISTEP ARE INPUT.  L, KOUNT ARE OUTPUT.
C--------------- X, ISTEP ARE MODIFIED.
      SUBROUTINE AVDAT (K,LEVEL,N,X,L,ISTEP,KOUNT)
      DIMENSION LEVEL(1),X(1),ISTEP(1),KOUNT(1)
C
C        ...............................................................
C
C        IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
C        C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
C        STATEMENT WHICH FOLLOWS.
C
C     DOUBLE PRECISION X
C
C        THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
C        APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
C        ROUTINE.
C
C        ...............................................................
C
C     CALCULATE TOTAL DATA AREA REQUIRED
C
      M=LEVEL(1)+1
      DO 105 I=2,K
  105 M=M*(LEVEL(I)+1)
C
C     MOVE DATA TO THE UPPER PART OF THE ARRAY X
C     FOR THE PURPOSE OF REARRANGEMENT
C
      N1=M+1
      N2=N+1
      DO 107 I=1,N
      N1=N1-1
      N2=N2-1
  107 X(N1)=X(N2)
C
C     CALCULATE MULTIPLIERS TO BE USED IN FINDING STORAGE LOCATIONS FOR
C     INPUT DATA
C
      ISTEP(1)=1
      DO 110 I=2,K
  110 ISTEP(I)=ISTEP(I-1)*(LEVEL(I-1)+1)
      DO 115 I=1,K
  115 KOUNT(I)=1
C
C     PLACE DATA IN PROPER LOCATIONS
C
      N1=N1-1
      DO 135 I=1,N
      L=KOUNT(1)
      DO 120 J=2,K
  120 L=L+ISTEP(J)*(KOUNT(J)-1)
      N1=N1+1
      X(L)=X(N1)
      DO 130 J=1,K
      IF(KOUNT(J)-LEVEL(J)) 124, 125, 124
  124 KOUNT(J)=KOUNT(J)+1
      GO TO 135
  125 KOUNT(J)=1
  130 CONTINUE
  135 CONTINUE
      RETURN
      END
C
C     ..................................................................
C
C        SUBROUTINE MEANQ
C
C        PURPOSE
C           COMPUTE SUM OF SQUARES, DEGREES OF FREEDOM, AND MEAN SQUARE
C           USING THE MEAN SQUARE OPERATOR.  THIS SUBROUTINE NORMALLY
C           FOLLOWS CALLS TO AVDAT AND AVCAL SUBROUTINES IN THE PER-
C           FORMANCE OF ANALYSIS OF VARIANCE FOR A COMPLETE FACTORIAL
C           DESIGN.
C
C        USAGE
C           CALL MEANQ (K,LEVEL,X,GMEAN,SUMSQ,NDF,SMEAN,MSTEP,KOUNT,
C                        LASTS)
C
C        DESCRIPTION OF PARAMETERS
C           K     - NUMBER OF VARIABLES (FACTORS). K MUST BE .GT. ONE.
C           LEVEL - INPUT VECTOR OF LENGTH K CONTAINING LEVELS (CATE-
C                   GORIES) WITHIN EACH VARIABLE.
C           X     - INPUT VECTOR CONTAINING THE RESULT OF THE SIGMA AND
C                   DELTA OPERATORS. THE LENGTH OF X IS
C                   (LEVEL(1)+1)*(LEVEL(2)+1)*...*(LEVEL(K)+1).
C           GMEAN - OUTPUT VARIABLE CONTAINING GRAND MEAN.
C           SUMSQ - OUTPUT VECTOR CONTAINING SUMS OF SQUARES.  THE
C                   LENGTH OF SUMSQ IS 2 TO THE K-TH POWER MINUS ONE,
C                   (2**K)-1.
C           NDF   - OUTPUT VECTOR CONTAINING DEGREES OF FREEDOM.  THE
C                   LENGTH OF NDF IS 2 TO THE K-TH POWER MINUS ONE,
C                   (2**K)-1.
C           SMEAN - OUTPUT VECTOR CONTAINING MEAN SQUARES.  THE
C                   LENGTH OF SMEAN IS 2 TO THE K-TH POWER MINUS ONE,
C                   (2**K)-1.
C           MSTEP - WORKING VECTOR OF LENGTH K.
C           KOUNT - WORKING VECTOR OF LENGTH K.
C           LASTS - WORKING VECTOR OF LENGTH K.
C
C        REMARKS
C           THIS SUBROUTINE MUST FOLLOW SUBROUTINE AVCAL
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           THE METHOD IS BASED ON THE TECHNIQUE DISCUSSED BY H. O.
C           HARTLEY IN 'MATHEMATICAL METHODS FOR DIGITAL COMPUTERS',
C           EDITED BY A. RALSTON AND H. WILF, JOHN WILEY AND SONS,
C           1962, CHAPTER 20.
C
C     ..................................................................
C
C---------------K, LEVEL, X ARE INPUT.  GMEAN, SUMSQ, NDF, SMEAN,
C--------------- KOUNT, LASTS, MSTEP ARE OUTPUT.
      SUBROUTINE MEANQ (K,LEVEL,X,GMEAN,SUMSQ,NDF,SMEAN,MSTEP,KOUNT,
     1                  LASTS)
      DIMENSION LEVEL(1),X(1),SUMSQ(1),NDF(1),SMEAN(1),MSTEP(1),
     1          KOUNT(1),LASTS(1)
C
C        ...............................................................
C
C        IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
C        C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
C        STATEMENT WHICH FOLLOWS.
C
C     DOUBLE PRECISION X,GMEAN,SUMSQ,SMEAN,FN1
C
C        THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
C        APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
C        ROUTINE.
C
C        ...............................................................
C
C     CALCULATE TOTAL NUMBER OF DATA
C
      N=LEVEL(1)
      DO 150 I=2,K
  150 N=N*LEVEL(I)
C
C     SET UP CONTROL FOR MEAN SQUARE OPERATOR
C
      LASTS(1)=LEVEL(1)
      DO 178 I=2,K
  178 LASTS(I)=LEVEL(I)+1
      NN=1
C
C     CLEAR THE AREA TO STORE SUMS OF SQUARES
C
      LL=(2**K)-1
      MSTEP(1)=1
      DO 180 I=2,K
  180 MSTEP(I)=MSTEP(I-1)*2
      DO 185 I=1,LL
  185 SUMSQ(I)=0.0
C
C     PERFORM MEAN SQUARE OPERATOR
C
      DO 190 I=1,K
  190 KOUNT(I)=0
  200 L=0
      DO 260 I=1,K
      IF(KOUNT(I)-LASTS(I)) 210, 250, 210
  210 IF(L) 220, 220, 240
  220 KOUNT(I)=KOUNT(I)+1
      IF(KOUNT(I)-LEVEL(I)) 230, 230, 250
  230 L=L+MSTEP(I)
      GO TO 260
  240 IF(KOUNT(I)-LEVEL(I)) 230, 260, 230
  250 KOUNT(I)=0
  260 CONTINUE
      IF(L) 285, 285, 270
  270 SUMSQ(L)=SUMSQ(L)+X(NN)*X(NN)
      NN=NN+1
      GO TO 200
C
C     CALCULATE THE GRAND MEAN
C
  285 FN=N
      GMEAN=X(NN)/FN
C
C     CALCULATE FIRST DIVISOR REQUIRED TO FORM SUM OF SQUARES AND SECOND
C     DIVISOR, WHICH IS EQUAL TO DEGREES OF FREEDOM, REQUIRED TO FORM
C     MEAN SQUARES
C
      DO 310 I=2,K
  310 MSTEP(I)=0
      NN=0
      MSTEP(1)=1
  320 ND1=1
      ND2=1
      DO 340 I=1,K
      IF(MSTEP(I)) 330, 340, 330
  330 ND1=ND1*LEVEL(I)
      ND2=ND2*(LEVEL(I)-1)
  340 CONTINUE
      FN1=N*ND1
      FN2=ND2
      NN=NN+1
      SUMSQ(NN)=SUMSQ(NN)/FN1
      NDF(NN)=ND2
      SMEAN(NN)=SUMSQ(NN)/FN2
      IF(NN-LL) 345, 370, 370
  345 DO 360 I=1,K
      IF(MSTEP(I)) 347, 350, 347
  347 MSTEP(I)=0
      GO TO 360
  350 MSTEP(I)=1
      GO TO 320
  360 CONTINUE
  370 RETURN
      END
C---------------SETUP ALLOWS CORRECTION OF INPUT DATA FROM TERMINAL.
C--------------- SEE WRITE UP BOTTOM PAGE 6.  ALL ARGS. ARE INPUT.
      SUBROUTINE SETUP(LINES,NOUTD,IND)
      DIMENSION X(16)
       DOUBLE PRECISION FIL
      DATA FIL/'JOB.DAT'/
      CALL OFILE (1,'JOB')
      DO 3 J=1,LINES
1     READ(IND,10,END=20)(X(I),I=1,16)
10    FORMAT(16A5)
3     WRITE(1,10)(X(I),I=1,16)
20    ENDFILE 1
C---------------TTYPTY RETURNS ZERO - TTY JOB, MINUS ONE - BATCH JOB
      CALL TTYPTY(IWHER)
      IF(IWHER.EQ.-1)GO TO 26
200   WRITE(NOUTD,21)
21    FORMAT(' OK? '/)
      READ(IND,22)HUH
22    FORMAT(A3)
      IF(HUH.EQ.'NO')GO TO 25
      IF(HUH.EQ.'YES')GO TO 26
      IF(HUH.NE.'OK')GO TO 200
      GO TO 26
25    CALL IFILE(1,'JOB')
      K=1
40    READ(1,10,END=41)(X(I),I=1,16)
      DO 27 I=1,16
      IF(X(I).NE.' ')GO TO 28
27    CONTINUE
28    DO 29 J=I,16
      IF(X(J).EQ.' ')GO TO 30
29    CONTINUE
      J=17
30    WRITE(NOUTD,31)K,(X(I),I=1,J-1)
31    FORMAT(' ',I3,2X,16A5)
      K=K+1
      GO TO 40
41    WRITE(NOUTD,45)
45    FORMAT(//)
      REWIND 1
      CALL DEFINE FILE(1,80,NEVER,FIL,0,0)
101   WRITE(NOUTD,42)
42    FORMAT('+CORRECTIONS? ',$)
      READ(IND,43)K,(X(I),I=1,16)
43    FORMAT(I,16A5)
      IF(K.EQ.0)GO TO 100
      WRITE(1#K,10)(X(I),I=1,16)
      REWIND 1
      GO TO 101
100   CALL RELEAS(1)
26    RETURN
      END

