C	WESTERN MICHIGAN UNIVERSITY
C	CROSS.F4 (FILE NAME ON LIBRARY DECTAPE)
C	CROSS, 1.1.6 (CALLING NAME, SUBLST NO.)
C	CROSS TABULATION WITH CHI-SQUARE
C	THIS PROGRAM WAS TAKEN FROM THE BOOK "DATA PROCESSING"
C	 BY K. JANDA.  SEVERAL MODIFICATIONS WERE MADE TO MAKE
C	 IT COMPATIBLE AND SOME ADDITIONAL FEATURES WERE ADDED BY
C	 SAM ANEMA.
C	LIBRARY DECTAPE PROGS. USED:  USAGE.MAC
C	FORWMU PROGS. USED:  TTYPTY, EXISTS, PRINTS, DEVCHG
C	APLIB PROGS. USED:  IO, GETFOR
C	INTERNAL SUBR. USED:  COEFT
C	ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
      DIMENSION LAD(200,40), MAX(40), KIND(40),KONST(40),MBRAK(16)
      DIMENSION FMT(48), NTAB(101,21), NRTOT(101), NCLT(21),LOB(16)
      DIMENSIONTAB(101,21),RTOT(101), CLT(21), PCENT(21)
      DIMENSION NDEEP(72), NDIM(4,72), MXTEMP(4)
      DIMENSION STORE(40,12)
      DIMENSION CKK(16)
C
C
      EQUIVALENCE (NTAB,TAB)
      NOUTD=-1
      IND=-4
      INP=2
      IOUT=3
C---------------TTYPTY RETURNS 0 - TTY JOB, MINUS ONE - BATCH JOB
      CALL TTYPTY(ICODE)
      WRITE(NOUTD,8745)
8745  FORMAT(20X,'PDP-10 VERSION OF NUCROS')
    7 FORMAT(40I2)
C	CALL USAGE('CROSS')
C
C     READ AND PRINTTITLES
C---------------NDV1, NDV2 ARE RETURNED.  OTHER ARGS. ARE INPUT.
C--------------- 1 MEANS OUTPUT? PRINTS, 0 - INPUT? PRINTS
      CALL IO(IOUT,NDV1,NOUTD,IND,1,ICODE)
500   CALL IO(INP,NDV2,NOUTD,IND,0,ICODE)
C---------------FMT, ISTD ARE RETURNED.  OTHER ARGS. ARE INPUT.
C--------------- 48=NO. OF OBJ. TIME FORMAT WORDS (3 LINES).  
C--------------- 1 MEANS I TYPE FORMAT ONLY.
      CALL GETFOR(NOUTD,IND,FMT,ISTD,48,1)
      REWIND 20
      REWIND 21
      REWIND 22
  800 FORMAT(1H1)
      WRITE(NOUTD,1002)
1002  FORMAT(' ENTER IDENTIFICATION'/)
      READ(IND,10)(CKK(I),I=1,16)
      WRITE(IOUT,5)(CKK(I),I=1,16)
    5 FORMAT(1X,16A5)
C
C     READ CONTROL CARDS 4-8
230   WRITE(NOUTD,1003)
1003  FORMAT(' ENTER PROBLEM CONTROL PARAMETERS'/)
      READ(IND,1001)NOBS,NVAR,NTABLE,IDVAR,NCODE,NCHI,NPCT,
     1NTAU
1001  FORMAT(10I)
      IF(NVAR.LT.2.OR.NVAR.GT.40)GO TO 217
      IF(NTABLE .LT.1.OR.NTABLE.GT.72)GO TO 218
      GO TO 219
217   WRITE(NOUTD,220)
220   FORMAT(' ILLEGAL NUMBER OF VARIABLES SPECIFIED.'/)
      IF(ICODE.EQ.-1)CALL EXIT
      GO TO 230
218   WRITE(NOUTD,221)
221   FORMAT(' ILLEGAL NUMBER OF TABLES SPECIFIED'/)
      IF(ICODE.EQ.-1)CALL EXIT
      GO TO 230
219   LCHI=0
      IF(NCHI.GT.0.OR.NTAU.GT.0)LCHI=1
   10 FORMAT(16A5)
      IF(NDV1.NE.'TTY')WRITE(IOUT,690)(FMT(I),I=1,48)
690   FORMAT(' FORMAT OF DATA'/(' ',16A5))
      LPCT = MIN0(NPCT,1)
      WRITE (IOUT,185) NOBS
      WRITE(NOUTD,1006)
1006  FORMAT(//' ENTER MAXIMUM VALUES'/)
      READ (IND,1012) (MAX(J),J=1,NVAR)
      WRITE (IOUT,9)(J,J=1,NVAR)
      WRITE(IOUT,423) (MAX(J),J=1,NVAR)
    6 FORMAT(I4,9I2)
  185 FORMAT(/20H TOTAL CASES USED =  I4)
  780 FORMAT (1H0)
  423 FORMAT(11H MAX.VALUE  ,40I3)
    9 FORMAT(//11H VARIABLE  ,  40I3)
C
C     READ + PRINT CONTOL CDS 9,10,11,12.. DETERMINE ALL TABLE DEPTHS
      WRITE(NOUTD,1007)
1007  FORMAT(//' ENTER TABLE CONTROL PARAMETERS'/)
1012  FORMAT(30I2)
      DO 1000 I=1,4
 1000 READ (IND,1012) (NDIM(I,K),K=1,NTABLE)
      DO472  K=1,NTABLE
      NDEEP(K) =2
      DO 472 I=3,4
      IF (NDIM(I,K))472,472,470
  470 NDEEP(K)=I
  472 CONTINUE
      WRITE (IOUT,424)
      IF(NTABLE-36)9001,9001,9002
 9001 WRITE(IOUT,781)
      GO TO 9003
 9002 WRITE (IOUT,426)
 9003 CONTINUE
      LTABLE = MIN0(NTABLE,36)
      DO 4 L=1,LTABLE
      K=L
      NDEPTH = NDEEP(K)
      WRITE (IOUT,8) K,(NDIM(I,K),I=1,NDEPTH)
      DO 550 I = 1,NDEPTH
      J = NDIM(I,K)
  550 MXTEMP(I) = MAX(J)
      WRITE(IOUT,425) (MXTEMP(I),I=1,NDEPTH)
      K=L+36
      IF (K .GT. NTABLE) GO TO 427
      NDEPTH = NDEEP(K)
      WRITE  (IOUT,428) K,(NDIM(I,K),I=1,NDEPTH)
      DO 429I=1,NDEPTH
      J=NDIM(I,K)
  429 MXTEMP(I) = MAX(J)
      WRITE (IOUT,430) (MXTEMP(I),I=1,NDEPTH)
      GO TO 4
  427 CONTINUE
    4 CONTINUE
    8 FORMAT (1X,9XI6,7X,4I4)
  424 FORMAT (//1H+,9X9HTABLE NO.,6X9HVARIABLES,7X14HMAXIMUM VALUES)
  425 FORMAT  (1H+,38X,4I4)
  426 FORMAT ('+'75X9HTABLE NO.,6X9HVARIABLES,7X14HMAXIMUM VALUES)
  428 FORMAT (1H+,74X,I6,7X,4I4)
  430 FORMAT ('+'104X,4I4)
  781 FORMAT (' ')
C
C     READ AND PRINT RECODE-CONTROL-CARDS, IF RECODING OPTED
      IF (NCODE) 20,20,19
19    WRITE(NOUTD,1014)
1014  FORMAT(' ENTER TWO RECODING LINES'/)
      READ (IND,1012) (KIND(J),J=1,NVAR)
      READ (IND,1012) (KONST(J),J=1,NVAR)
      DO 400 J=1,NVAR
      NOPER = KIND(J)
      IF (NOPER-3) 400,400,403
403   WRITE(NOUTD,1015)
1015  FORMAT(' ENTER GROUP BOUNDARIES'/)
      READ(IND,1019) MBRAK
1019  FORMAT(16I)
      IF(NOPER.EQ.8)GO TO 650
      DO 651 I=0,15
651   LOB(I+1)=I
      GO TO 652
650   WRITE(NOUTD,653)
653   FORMAT(' ENTER CATEGORIZATION'/)
      READ(IND,1019)LOB
652   WRITE (21) MBRAK,LOB
  400 CONTINUE
      WRITE(IOUT,186)
      REWIND 21
      DO 2200 J=1,NVAR
      NOPER=KIND(J)
      IF(NOPER)2200,2200,2202
2202  IF(NOPER-3)2201,2201,2203
2201  WRITE(IOUT,188)J,NOPER,KONST(J)
      GO TO 2200
2203  READ(21)MBRAK,LOB
      IF(NOPER.EQ.4)WRITE(IOUT,187)J,NOPER,MBRAK
      IF(NOPER.GT.4.AND.NOPER.LT.8)WRITE(IOUT,9004)J,NOPER,KONST(J),MBRA
     1K
      IF(NOPER.EQ.8)WRITE(IOUT,660)J,NOPER,MBRAK,LOB
2200  CONTINUE
  186 FORMAT (27H1RECODING OPTION CALLED FOR // 4X,67H VAR.NO. CODE  
     X CONSTANT     RECODED CATEGORIES (UPPER LIMITS)     )
  188 FORMAT (8XI2,5XI1,4XI3,10X4HNONE)
  187 FORMAT ( 8XI2,5XI1,4X4HNONE,16I3)
 9004 FORMAT(8XI2,5XI1,4XI3,1X,16I3)
660   FORMAT(8X,I2,5X,I1,4X,'NONE',16I3/24X,16I3)
      REWIND 21
C
C     READ AND PRINT VAR NAMES IF OPTED
   20 IF (IDVAR) 480,480,495
495   WRITE(NOUTD,1008)
1008  FORMAT(' ENTER VARIABLE NAMES'/)
      DO 490 I=1,NVAR
490   READ(IND,82)(STORE(I,J),J=1,12)
      WRITE(IOUT,83)
      DO 1009 I=1,NVAR
1009  WRITE(IOUT,57)I,(STORE(I,J),J=1,12)
   83 FORMAT (8H1VAR.NO.//)
   82 FORMAT(12A5)
   57 FORMAT(I5,2X,12A5)
C
C     INITIALIZE FOR FIRST BATCH
  480 NBATCH = MIN0 (NOBS,200)
      NOBS = NOBS- NBATCH
      IBATCH = 1
      IF(NDV2.EQ.'TTY')GO TO 214
      WRITE(NOUTD,3576)
      GO TO 215
214   IF(ISTD.EQ.1)WRITE(NOUTD,216)
216   FORMAT(/' ENTER DATA (AT MOST 10 PER LINE)'/)
      IF(ISTD.EQ.0)WRITE(NOUTD,3577)
3576  FORMAT(/' YOUR DATA IS BEING READ. PLEASE WAIT.'/)
3577  FORMAT(/' ENTER DATA'/)
215   IOSW = 24
      LINE = 60
C
C
C      READ A BATCH OF DATA CARDS
299   DO 3817 N=1,NBATCH
3817  READ(INP,FMT)(LAD(N,J),J=1,NVAR)
C
C     RECODE VARIABLES IF OPTED
  302 IF (NCODE) 485,485,483
  483 DO 32 J = 1,NVAR
      NOPER = KIND(J)
      IF(NOPER) 32,32,176
  176  GO TO (23,24,25,52,52,52,52,52 ), NOPER
C
C     DO TYPE 1,2,OR 3 RECODING
   23 DO 1201 I=1,NBATCH
 1201 LAD(I,J) = LAD(I,J)+KONST (J)
      GO TO 32
   24 DO1203 I = 1,NBATCH
 1203 LAD(I,J) = MAX0(LAD(I,J)-KONST(J),0)
      GO TO 32
   25 DO1202 I= 1,NBATCH
 1202 LAD(I,J)= LAD(I,J) /KONST(J)
      GO TO 32
C
C     DO TYPE 4 REC DING
   52 READ (21)  MBRAK,LOB
      DO 29 I = 1,NBATCH
      DO 28 K = 1,16
      IF(LAD(I,J)-MBRAK(K)) 27,27,28
   28 CONTINUE
      LAD(I,J) = 20
      GO TO 29
   27 LAD(I,J) =LOB(K)
   29 CONTINUE
      IF(NOPER.EQ.8)GO TO 32
      IF(NOPER-4) 32,32,7638
 7638 NOPER= NOPER-4
      GO TO 176
   32 CONTINUE
      REWIND 21
C
C     DEFINE TABLE DEPTH,MAX VAL, VARS TO BE USED...FOR EACH TABLE
  485 DO 115 KT=1,NTABLE
      NDEPTH = NDEEP(KT)
      J1 = NDIM(1,KT)
      J2 = NDIM(2,KT)
      MAXA = MAX(J1) +1
      MAXB = MAX(J2) +1
      MAXC = 0
      MAXD = 0
      LREQ = 10 + MAXB + 2*LCHI
      LREQPC = LPCT * (MAXB-1+6)
      GO TO (36,36,33,34),NDEPTH
   34 J4 = NDIM(4,KT)
      MAXD = MAX(J4)
   33 J3 = NDIM(3,KT)
      MAXC = MAX(J3)
      LREQ = LREQ +2
C
   36 DO 115 L3 = 0,MAXC
      DO 115 L4 = 0,MAXD
      IF (IBATCH-1) 41,41,42
C     ZERO OUT TABLE AREA IF FIRST BATCH
   41 DO 43 I=1,MAXB
      DO 43 J=1,MAXA
   43 NTAB(I,J) = 0
      GO TO 44
C     OTHERWISE READ INTO TABLE-AREA STORED-RESULTS FROM PREVIOUS BATCH(ES)
   42 IF (IOSW-24) 46,46,47
   46 READ (20) ((NTAB(I,J),J=1,MAXA),I=1,MAXB)
      GO TO 44
   47 READ (22) ((NTAB(I,J),J=1,MAXA),I=1,MAXB)
C
C     DO A TABLE.  TEST 3RD AND 4TH-DIMENSION VARIABLES FOR APPROP VALUE
C     BEFORE COUNTING THE OBSERVATION (WHERE 3RD OR 4TH DIM IS USED)
   44 DO 45 N=1,NBATCH
      IF (NDEPTH-3) 40,35,39
   39 IF (LAD(N,J4) -L4) 45,35,45
   35 IF (LAD(N,J3) - L3) 45,40,45
   40 L1 = LAD(N,J1) +1
      L2 = LAD(N,J2) +1
	IF (L1 .LT. 1 .OR. L2 .LT. 1) GO TO 45
      IF (L1.GT.MAXA  .OR.  L2.GT.MAXB) GO TO 45
      NTAB(L2,L1) = NTAB(L2,L1) +1
   45 CONTINUE
C
C     IF MORE OBSERVATIONS REMAIN, STORE THIS TABLE (RESULTS TO DATE)
      IF (NOBS) 75,75,48
   48 IF (IOSW-24) 50,50,60
   50 WRITE (22) ((NTAB(I,J),J=1,MAXA),I=1,MAXB)
      GO TO 115
   60 WRITE (20) ((NTAB(I,J),J=1,MAXA),I=1,MAXB)
      GO TO 115
C
C
C     IF NO MORE OBSERVATIONS WRITE OUT THIS TABLE
   75 LINE= LINE + LREQ
      IF(LINE+2*LREQPC-58)431,431,811
  431 WRITE(IOUT,780)
      GO TO 76
  811 WRITE (IOUT,800)
      LINE=LREQ-1
   76 GO TO (822,822,823,824),NDEPTH
C     TWO-DIMENSION TITLE
  822 IF(IDVAR)522,522,523
  523 WRITE(IOUT,524)KT,J1,(STORE(J1,J),J=1,12),J2,(STORE(J2,J),J=1,12)
     1  ,MAXA,MAXB
  524 FORMAT(7X,11H TABLE NO.  ,I3,10X,13HVARIABLE NO.  ,I2,2X,  12A5,/
     1  25X,6HVRS.  ,13HVARIABLE NO.  ,I2,2X,12A5,  / 7X,13HTABLE SIZE =
     2  I3,5H  BY ,I3)
      GO TO 86
  522 WRITE (IOUT,80) KT,J1,J2,MAXA,MAXB
   80 FORMAT (    10X,11H TABLE NO. I3,5X,13HVARIABLE NO.  I2,3X,
     1  9HVRS.  NO. ,  I2//13X,13HTABLE SIZE  =  I3,5H BY   I3)
      GO TO 86
C     THREE DIMENSIOANL TITLE
  823 IF (IDVAR)525, 525,526
  526 WRITE(IOUT,527) KT,J1,(STORE (J1,J),J=1,12),J2,(STORE(J2,J),J=1,12
     1),J3,(STORE(J3,J),J=1,12), MAXA,MAXB,L3
  527 FORMAT(7X,11H TABLE NO. ,  I3,10X,13HVARIABLE NO.  I2,2X,12A5,
     1 2(/25X,6HVRS   ,13HVARIABLE NO.  I2,2X,12A5)   //7X,
     213HTABLE SIZE =  I3,5H  BY I3,10X,16H 3RD DIMENSION = I3)
      GO TO 86
  525 WRITE (IOUT,77) KT,J1,J2,J3, MAXA,MAXB,L3
   77 FORMAT(10X,11HTABLE NO.  I3,5X,13HVARIABLE NO  I2,3X,9HVRS NO.  
     1I2,3X,9HVRS. NO. I2,///13X,13HTABLE SIZE =  I3,5H BY  I3//13X,
     216H3RD DIMENSION =  I3)
      GO TO 86
C     FOUR DIMENSION TITLE
  824 IF(IDVAR)528,528,529
  529 WRITE(IOUT,530)KT,J1,(STORE(J1,J),J=1,12),J2,(STORE(J2,J),J=1,12),
     1J3,(STORE(J3,J),J=1,12),J4,(STORE(J4,J),J=1,12),MAXA,MAXB,L3,L4
  530 FORMAT(7X,11H TABLE NO.  I3,10X,13HVARIABLE NO.  I2,2X, 12A5,
     13(/25X,6HVRS.  ,13HVARIABLE NO. I2,2X, 12A5) //7X,
     213HTABLE SIZE =  I3,5H BY  I3,10X,16H3RD DIMENSION =  I3,5X,
     316H4TH DIMENSION =   I3)
      GO TO 86
  528 WRITE(IOUT,72) KT,J1,J2,J3,J4,MAXA,MAXB,L3,L4
   72 FORMAT (10X,11H TABLE NO.  I3,5X,13HVARIABLE NO.  I2,3X,9HVRS. 
     1NO.  I2,3X,9HVRS NO.  I2,3X,9HVRS NO.  I2//13X,13HTABLE SIZE = 
     2I3,5H BY  I3//13X,16H3RD DIMENSION = I3,5X,16H4TH DIMENSION = I3)
   86 CONTINUE
C
C
C  PRINT BODY OF MAP,COMPUTE AND PRINT COLUMN AND ROW TOTALS
   84 MSUB = MAXA-1
      WRITE (IOUT,91)  (IH,IH=0,MSUB)
      WRITE (IOUT,781)
      NGRTO=0
      DO 85 J=1,MAXA
   85 NCLT(J)=0
      DO 100 I=1,MAXB
      ISUB = I-1
      NRTOT(I)=0
      DO  90 J=1,MAXA
      NRTOT(I)= NRTOT(I) + NTAB(I,J)
   90 NCLT(J)= NCLT(J)+NTAB (I,J)
      WRITE(IOUT,105) ISUB,NRTOT(I),(NTAB(I,J),J=1,MAXA)
  100 NGRTO = NGRTO+NRTOT(I)
      WRITE (IOUT,106) NGRTO,(NCLT(J),J=1,MAXA)
   91 FORMAT (/17X,3HTOT,3X,21I5)
  105 FORMAT(9X,I3,I8,3X,21I5)
  106 FORMAT (/9X5HTOTAL,I6,3X,21I5)
C
C
C DO CHISQARES +/OR PERCENTS +/OR COEFFICIENTS IF OPTED
C GET FLOATING _PT EQUIVALENTS..DEDUCT ZEROROW,COLUMN FROM TOTALS
      IF(NCHI +NPCT + NTAU) 115,115,404
404   DO 407 I=1,MAXB
      RTOT(I)=NRTOT(I)
      DO 407 J=1,MAXA
  407 TAB(I,J) = NTAB (I,J)
      DO 405 I=1,MAXA
405   CLT(I)=NCLT(I)
      GRTO=NGRTO
C
      IF(NCHI+NTAU) 410,410,605
  605 IF (NCHI) 606,606,408
  606 WRITE (IOUT,781)
      GO TO 304
  408 CHISQ = 0.0
      IF(NGRTO.NE.0)GO TO 1071
      CHISQ=0.
      CRXC=0.
      GO TO 1070
1071  DO 409 I=1,MAXB
      DO409 J=1,MAXA
      A= (CLT(J)*RTOT(I)) /GRTO
      B= (TAB(I,J)-A )**2
      IF(A.NE.0.)GO TO 409
      A=1.
      B=0.
409   CHISQ = CHISQ + B/A
      CRXC=(  CHISQ/(CHISQ+GRTO)) **.5
1070  WRITE(IOUT,608) CHISQ,CRXC
  608 FORMAT(/1H+,6X,12HCHI SQUARE =  ,F8.3,3X,4HC =  ,F5.3)
      IF(MAXB.NE.2.OR.MAXA.NE.2)GO TO 631
      A=TAB(1,1)
      B=TAB(1,2)
      C=TAB(2,1)
      D=TAB(2,2)
      C22=GRTO*(ABS(A*D-B*C)-GRTO/2)**2/((A+B)*(C+D)*(A+C)*(B+D))
      WRITE(IOUT,632)C22
632   FORMAT(4X,'2X2 CORRECTED CHI-SQUARE = ',F8.3)
631   IF(NTAU) 604,604,304
  604 WRITE(IOUT,781)
      GO TO 410
C
C
C      DO ALL CORRELATION COEFFICIENTS
C     NTAU IS NOT USED BY THE SUBROUTINE BUT MAY BE USED TO
C      EXERCISE OPTIONS IN ALTERNATE SUBROUTINES
  304 CALL COEFT (TAB,MAXB,MAXA,NTAU)
C
C
C     DO COLUMN PERCENTAGES
  410 IF (NPCT) 115,115,411
  411 LINE = LINE +LREQPC
      IF (LINE-58) 412,412,413
  412 WRITE (IOUT,108)
      GO TO 414
  413 WRITE (IOUT,800)
      WRITE (IOUT,108)
      WRITE (IOUT,113) (IH,IH=0,MSUB)
      WRITE (IOUT,781)
      LINE = LREQPC + 3
  108 FORMAT(//7X,40HPERCENTS BY COLUMN FROM THE ABOVE MATRIX     //)
  113 FORMAT (17X,3HTOT,8X,20I5 )
C
  414 DO 416 I = 1,MAXB
      ISUB = I-1
      DO 415 J = 1,MAXA
      IF(CLT(J).NE.0.)GO TO 1080
      PCENT(J)=0.
      GO TO 415
1080  PCENT(J) = TAB(I,J) / CLT(J)
415   CONTINUE
  416 WRITE(IOUT,109) ISUB,NRTOT(I),  (PCENT(J),J = 1,MAXA)
      WRITE (IOUT,111) NGRTO,(NCLT(J), J=1,MAXA )
  109 FORMAT (9X,I3,I8,9X,2P20F5.1)
  111 FORMAT(/9X,5HTOTAL, I6,8X,20I5 )
C
C     DO ROW PERCENTS
      LINE = LINE+ LREQPC
      IF(LINE-58) 417,417,418
  417 WRITE( IOUT,114)
      GO TO 419
  418 WRITE( IOUT,800)
      WRITE (IOUT,114)
      WRITE (IOUT,113) (IH,IH=0,MSUB)
      WRITE (IOUT,781)
      LINE = LREQPC + 3
  114 FORMAT(//7X,38HPERCENTS BY ROW FROM THE ABOVE MATRIX        //)
C
  419 DO 421 I= 1,MAXB
      ISUB = I-1
      DO 420 J= 1,MAXA
      IF(RTOT(I).NE.0.)GO TO 1081
      PCENT(J)=0.
      GO TO 420
1081  PCENT(J) = TAB (I,J) / RTOT(I)
420   CONTINUE
  421 WRITE (IOUT,109) ISUB,NRTOT(I),(PCENT(J),J=1,MAXA)
      WRITE (IOUT,111) NGRTO, (NCLT(J),J=1,MAXA)
  115 CONTINUE
C
C
C     PREPARE PROCESS NEXT BATCH IF ANY
      IF (NOBS) 500,500,130
  130 NBATCH = MIN0 (NOBS,200)
      NOBS = NOBS-NBATCH
      IBATCH = IBATCH + 1
      IF(IOSW-24) 12,12,13
   12 IOSW = 42
      GO TO 15
   13 IOSW = 24
   15 REWIND 20
      REWIND 22
      GO TO 299
       END
C---------------NTAU APPARENTLY NOT USED.  OTHER ARGS. ARE INPUT.
      SUBROUTINE COEFT(A,II,JJ,NTAU)
      DIMENSION A(101,21), TI(101), TJ(21)
      IOUT=3
      SUM = 0.
      DO 2  I = 1,II
      DO 2 J = 1,JJ
    2 SUM=SUM+A(I,J)
      P=0.
      IX = II-1
      JX = JJ-1
      DO 10 I=1,IX
      DO 10 J=1,JX
      IY=I+1
      JY=J+1
      DO 10 K=IY,II
      DO 10 L=JY,JJ
   10 P=P+A(I,J)*A(K,L)
      Q=0.
      DO 11 I=1,IX
      DO 11J=2,JJ
      IN=I+1
      JY=J-1
      DO 11 K=IN,II
      DO 11 L=1,JY
   11 Q = Q+A(I,J)*A(K,L)
      S=P-Q
      GAMMA=(P-Q)/(P+Q)
  100 DO 21 I=1,II
      TI(I)=0.
      DO21 J=1,JJ
   21 TI(I)=TI(I)+A( I,J)
      DO 22 J=1,JJ
      TJ(J)=0.
      DO22 I=1,II
   22 TJ(J)=TJ(J)+A(I,J)
      I1=0
      DO 105 I=1,II
105   IF(TI(I).NE.0.0)I1=I1+1
      J1=0
      DO 106 J=1,JJ
106   IF(TJ(J).NE.0.0)J1=J1+1
   82 IF (I1-J1) 86,85,86
   85 IPT = 2
      GO TO 84
   86 IPT = 3
   84 CONTINUE
      XU=0.
      DO55 J=1,JJ
   55 XU=XU + TJ(J)**2
      XU=.5*(SUM**2-XU)
      YU=0.
      DO56 I=1,II
   56 YU=YU+TI(I)**2
      YU=.5*(SUM**2-YU)
      DYX = S/XU
      DXY=S/YU
      IF(IPT-2) 57,57,102
   57 CONTINUE
      T=0.
      DO 23 I=1,II
   23 T=T+.5*TI(I)*(TI(I)-1.)
      U=0.
      DO24 J=1,JJ
   24 U=U+.5*TJ(J)*(TJ(J)-1.)
      BI=(SUM*(SUM-1.) *.5-T)**(.5)
      BJ=(SUM*(SUM-1.) *.5-U)**(.5)
      TAUB=S/(BI*BJ)
  104 WRITE (IOUT,301) TAUB,GAMMA,DXY,DYX
      RETURN
  102 IF(I1-J1) 5,5,6
    5 AM=I1
      GO TO 12
    6 AM=J1
   12 TAUC=S/((.5*SUM**2*(AM-1.)/AM))
      WRITE(IOUT,302) TAUC,GAMMA,DXY,DYX
      RETURN
  301 FORMAT(1H+,41X,8HTAU-B =   ,F6.3,3X,8HGAMMA =  ,F6.3,3X,6HDXY = ,
     1F6.3,3X, 6HDYX =  ,F6.3)
  302 FORMAT(1H+,41X,8HTAU-C =   ,F6.3,3X,8HGAMMA =  ,F6.3,3X,6HDXY = ,
     1F6.3,3X, 6HDYX =  ,F6.3)
      END
