      SUBROUTINE MDL3(K,T,MATRIX,R,N,AVAHAT,PHI,THCHI,YRONE,HEADER)
      IMPLICIT DOUBLE PRECISION(A-G,O-Z)
      DOUBLE PRECISION INFO,NEXT,HIP,HIL,LOWP,LOWL,INVERS,LAMBDA,N,
     1 MATRIX,LOWF,HIF
      DOUBLE PRECISION MLS,LOWMLS,HIMLS
      INTEGER YEAR,YRONE,BANDS,S,SS,T,IMAT
      COMMON/PATCH/TEMP1,TEMP2  ,TEMP3
      DIMENSION MATRIX(21,21),R(20),N(20),THCHI(100)
      CHARACTER HEADER*80, FORMT*30
      COMMON /MDL3CM/ ISAVE(20),A(20),PD(2),AHAT(20),CHISQS(21
     1,21),INVERS(2,2),INFO(2,2),NEXT(2),EE(21,21),IMAT(20),
     2 PH(20), TMSPAC(1545)
      LAMBDA=(AVAHAT/100.)/(1.0-PHI)
      SUMDF=0.0
      DEGF=0.0
      ITER=0.0
      INVERS(1,1)=0.0
      INVERS(2,2)=0.0
      VARF=0.0
      BIGN=0.0
      BIGT=0.0
      Q=0.0
      NEXT(2)=0.0
      DO 5 I=1,T
      ISAVE(I)=0
      A(I)=0.0
      BIGT=BIGT+R(I)
      BIGN=BIGN+N(I)
      A(I)=N(I)
    5 CONTINUE
      DO 10 I=1,T
      DO 10 J=1,K
      Q=Q+(J-I)*MATRIX(I,J)
   10 CONTINUE
      EXPQ=Q
      EXPT=BIGT
      CALL MOD3W1
      WRITE (6,65) HEADER
   65 FORMAT (/' ',20X,A///)
      WRITE (6,70)
   70 FORMAT (' ','YEAR NUMBER',35X,'RECOVERY MATRIX')
      WRITE (6,75)
   75 FORMAT (' ',5X,'BANDED')
      WRITE (6,80)
   80 FORMAT ('+','____ ______',35X,'_______________'/)
      YEAR=YRONE-1
      DO 95 I=1,T
      BANDS=A(I)
      YEAR=YEAR+1
      DO 85 JJ=1,K
      IMAT(JJ)=MATRIX(I,JJ)
   85 CONTINUE
      IS=(I-1)*6+1
      WRITE(FORMT,90) IS
      WRITE (6,FORMT) YEAR,BANDS,(IMAT(JJ),JJ=I,K)
   90 FORMAT ('(1X,I4,1X,I5,',I2,'X,20I6)')
   95 CONTINUE
      WRITE (6,100)
  100 FORMAT (///' ','INTERMEDIATE STATISTICS'/)
      WRITE (6,105) BIGN,BIGT,Q
  105 FORMAT (' ','N = ',F6.0,5X,'T = ',F6.0,5X,'Q = ',F6.0/)
  110 CONTINUE
      IF (NEXT(2).GT.1.0) GO TO 310
      ITER=ITER+1
      IF (ITER.GT.20) GO TO 310
      SUM1=0.0
      SUM2=0.0
      SUM3=0.0
      SUM4=0.0
      SUM5=0.0
      DO 115 I=1,T
      SUM1=SUM1+(N(I)-R(I))/(1.-LAMBDA*(1.-PHI**(K-I+1)))
      SUM2=SUM2+(N(I)-R(I))*(K-I+1)*(PHI**(K-I))/(1.-LAMBDA*(1.-PHI**(K-
     1I+1)))
      SUM3=SUM3+(N(I)*(K-I+1)*PHI**(K-I))/(1.-LAMBDA*(1.-PHI**(K-I+1)))
      SUM4=SUM4+(N(I)/(1.-LAMBDA*(1.-PHI**(K-I+1))))
      SUM5=SUM5-LAMBDA*N(I)*(K-I+1)*(K-I)*(PHI**(K-I-1))
      AA=N(I)*LAMBDA*LAMBDA*(K-I+1)*(K-I+1)
      AA=AA*(PHI**(2.0*(K-I)))
      AA=AA/(1.0-LAMBDA*(1.0-PHI**(K-I+1)))
      SUM5=SUM5+AA
  115 CONTINUE
      SUM3=-1*SUM3
      SUM4=SUM4-BIGN
      PD(1)=(BIGN/LAMBDA)-(1./LAMBDA)*SUM1
      PD(2)=((-1*BIGT)/(1.-PHI))+(Q/PHI)+(LAMBDA*SUM2)
      INFO(1,2)=SUM3
      INFO(2,1)=SUM3
      INFO(1,1)=((1./LAMBDA)**2)*SUM4
      INFO(2,2)=(EXPT/((1.0-PHI)**2))+(EXPQ/PHI**2)+SUM5
C     COMPUTE THE INVERSE OF THE INFORMATION MATRIX
      DET=INFO(1,1)*INFO(2,2)-(INFO(2,1)*INFO(1,2))
      INVERS(1,1)=INFO(2,2)/DET
      INVERS(2,2)=INFO(1,1)/DET
      INVERS(2,1)=(-1.*INFO(2,1))/DET
      INVERS(1,2)=(-1.*INFO(1,2))/DET
C     THE INTERATION IS:
      NEXT(1)=LAMBDA+((INVERS(1,1)*PD(1))+(INVERS(1,2)*PD(2)))
      NEXT(2)=PHI+((INVERS(2,1)*PD(1))+(INVERS(2,2)*PD(2)))
      IF (DABS(NEXT(1)-LAMBDA).LT.0.00001.AND.DABS(NEXT(2)-PHI).LT.0.000
     101) GO TO 120
      LAMBDA=NEXT(1)
      OLD=PHI
      PHI=NEXT(2)
      IF (PHI.GT.1.0) PHI=(1.0+OLD)*0.5
      NEXT(2)=PHI
      GO TO 110
  120 NEXT(1)=NEXT(1)*100.
      NEXT(2)=NEXT(2)*100.
      STDEVP=TEMP2
      IF (INVERS(1,1).GT.0.0) STDEVL=DSQRT(INVERS(1,1))*100.
      CVP=STDEVP/NEXT(2)*100.
      CVL=STDEVL/NEXT(1)*100.
      LOWP=(NEXT(2)-(1.96*STDEVP))
      HIP=(NEXT(2)+(1.96*STDEVP))
      LOWL=(NEXT(1)-(1.96*STDEVL))
      HIL=(NEXT(1)+(1.96*STDEVL))
      F=(100.-NEXT(2))*(NEXT(1)/100.)
      COVAR=INVERS(2,1)
      F=F/100.
      AMORT=1.-(NEXT(2)/100.)
      VARF=(F/AMORT)**2*INVERS(2,2)-2.*F*COVAR+AMORT**2*INVERS(1,1)
C     COVFS=(-1.*(F/AMORT))*INVERS(2,2)+COVAR*AMORT
      COVFS=TEMP3
      F=F*100.
      STDEVF=TEMP1
      LOWF=F-(1.96*STDEVF)
      HIF=F+(1.96*STDEVF)
      CVF=(STDEVF/F)*100.
      WRITE (6,125)
  125 FORMAT (//' ',1X,'PARAMETER',31X,'ESTIMATE (%)',07X,'STD. ERR. (%)
     1',06X,'COEF. VARIAT. (%)',5X,'95% CONFIDENCE INTERVAL')
      WRITE (6,130)
  130 FORMAT (' ',124('-'))
      WRITE (6,135) NEXT(2),STDEVP,CVP,LOWP,HIP
  135 FORMAT (/' ','SURVIVAL RATE        (S)  ',18X,F5.2,13X,F5.2,15X,F5
     1.2,15X,F5.2,' -- ',F5.2/)
      WRITE (6,140) F,STDEVF,CVF,LOWF,HIF
  140 FORMAT (1X,'RECOVERY RATE        (F)             ',7X,F5.2,13X,
     1 F5.2,15X,F5.2,15X,F5.2,' -- ',F5.2/)
      WRITE (6,130)
      CORR=(INVERS(2,1))/((STDEVP*STDEVL)/10000.)
C     WRITE (6,160) CORR
C 160 FORMAT (//' ','CORRELATION(S,LAMBDA) = ',F10.8/)
      CORRFS=TEMP3/((TEMP2*TEMP1)/10000.0)
      WRITE (6,145) CORRFS
  145 FORMAT (' ','CORRELATION(S,F) = ',F10.8//)
      MLS=-1.0/DLOG(NEXT(2)/100.0)
      SDMLS=(MLS**2*(TEMP2/100.))/(NEXT(2)/100.)
      LOWMLS=-1.0/DLOG(LOWP/100.)
      HIMLS=-1.0/DLOG(HIP/100.)
      IF (HIP.GT.100.0) HIMLS=99999999.
      WRITE (6,150) MLS
  150 FORMAT (//' ','MEAN LIFE SPAN AS AN ADULT = ',F7.2)
      WRITE (6,155) SDMLS
  155 FORMAT (' ','STANDARD ERROR OF MEAN LIFE SPAN = ',F7.2)
      IF(HIMLS.NE.99999999.) GO TO 159
      WRITE (6,157) LOWMLS
  157 FORMAT(' ','95% CONFIDENCE INTERVAL OF LIFE SPAN ',
     1F7.2,' - *******'//)
      GO TO 162
  159 WRITE (6,160) LOWMLS,HIMLS
  160 FORMAT (' ','95% CONFIDENCE INTERVAL OF LIFE SPAN ',F7.2,' - ',F7.
     12//)
  162 CONTINUE
      S=T
      SS=T-1
      J=K
      DO 165 I=1,S
      AHAT(I)=F/100.
  165 CONTINUE
      DO 170 I=1,SS
      PH(I)=NEXT(2)/100.
  170 CONTINUE
      DO 175 JJ=1,J
      DO 175 I=1,S
      EE(I,JJ)=0.0
      CHISQS(I,JJ)=0.0
  175 CONTINUE
      DO 180 II=1,S
      DO 180 JJ=II,J
      EE(II,JJ)=A(II)*PH(1)**(JJ-II)*AHAT(1)
  180 CONTINUE
      DEGF=((S*J)-0.5*S*(S-1))-2
      DO 200 I=1,T
      DO 185 M=1,J
      IF (EE(I,J-M+1).GT.2.) GO TO 190
      EE(I,J)=EE(I,J)+EE(I,J-M)
      EE(I,J-M)=0.0
      IF (EE(I,J).GT.2.) GO TO 195
  185 CONTINUE
  190 ISAVE(I)=M-1
      GO TO 200
  195 ISAVE(I)=M
C     ISAVE(I) = NUMBER OF DEGREES OF FREEDOM LOST IN ROW I
  200 CONTINUE
      SUMDF=0.0
      DO 205 I=1,T
      SUMDF=SUMDF+ISAVE(I)
  205 CONTINUE
      DEGF=DEGF-SUMDF
      SS=S-1
      WRITE (6,210)
  210 FORMAT (///1X,'MATRIX OF EXPECTED VALUES -- ASSUMING CONSTANT SUR'
     1,'VIVAL AND RECOVERY RATES   (MODEL 3)'/)
      DO 220 I=1,S
      IS=(I-1)*6+7
      WRITE(FORMT,215) IS
      WRITE (6,FORMT) (EE(I,JJ),JJ=I,J)
  215 FORMAT ('(',I2,'X,20F6.1)')
  220 CONTINUE
      DO 230 I=1,T
      MM=ISAVE(I)
      DO 225 M=1,MM
      IF (ISAVE(I).NE.0) MATRIX(I,J)=MATRIX(I,J)+MATRIX(I,J-M)
      IF (ISAVE(I).NE.0) MATRIX(I,J-M)=0.0
  225 CONTINUE
  230 CONTINUE
      JPLUS=J+1
      DO 240 II=1,S
      POOL1=0.0
      POOL2=0.0
      DO 235 JJ=1,J
C     IF(J .EQ. 20) GO TO 1726
      POOL1=POOL1+EE(II,JJ)
      POOL2=POOL2+MATRIX(II,JJ)
  235 CONTINUE
      EE(II,J+1)=POOL1
      MATRIX(II,J+1)=POOL2
  240 CONTINUE
      CHIADD=0.0
      DO 245 I=1,S
      DO 245 JJ=I,JPLUS
      IF (JJ.NE.JPLUS.AND.EE(I,JJ).GT.0.) CHISQS(I,JJ)=((MATRIX(I,JJ)-EE
     1(I,JJ))**2)/EE(I,JJ)
      IF (JJ.EQ.JPLUS.AND.EE(I,JJ).GT.0.) CHISQS(I,JJ)=((MATRIX(I,JJ)-EE
     1(I,JJ))**2)/(A(I)-EE(I,JJ))
      IF (EE(I,JJ).GT..01) CHIADD=CHIADD+CHISQS(I,JJ)
  245 CONTINUE
      WRITE (6,250)
  250 FORMAT (//1X,'MATRIX OF CHI-SQUARE VALUES -- ASSUMING CONSTANT SU'
     1,'RVIVAL AND RECOVERY RATES   (MODEL 3)'/)
      DO 260 I=1,S
      IS=(I-1)*6+7
      WRITE(FORMT,255) IS
  255 FORMAT ('(',I2,'X,20F6.2)')
      WRITE (6,FORMT) (CHISQS(I,JJ),JJ=I,JPLUS)
  260 CONTINUE
      IF (SUMDF.GT.0.0) WRITE (6,265)
  265 FORMAT (/1X,'(FREQUENCIES WERE COMBINED WHERE EXPECTED VALUES WER'
     1,'E SMALL)'/)
      WRITE (6,270)
  270 FORMAT (//1X,'TEST OF THE NULL HYPOTHESIS THAT THE DATA FIT THE M'
     1,'ODEL ASSUMING CONSTANT SURVIVAL AND RECOVERY RATES'/)
      WRITE (6,275) CHIADD
  275 FORMAT (' ','   CHI-SQUARE VALUE (SAMPLE) = ',F8.2)
      WRITE (6,280)
  280 FORMAT (' ','   THEORETICAL CHI-SQUARE VALUE AT THE 5% LEVEL = ')
      JJ=DEGF
      IF (DEGF-100.0) 285,285,300
  285 JJ=DEGF
      WRITE (6,290) THCHI(JJ)
  290 FORMAT ('+',49X,F8.2)
      WRITE (6,295) JJ
  295 FORMAT (' ','   DEGREES OF FREEDOM = ',I5)
      GO TO 340
  300 THECHI=.5*(1.64+DSQRT((2.*DEGF)-1))**2
      WRITE (6,290) THECHI
      WRITE (6,295) JJ
      GO TO 340
  310 WRITE (6,315)
  315 FORMAT (1X,'THE ITERATIVE PROCEDURE FOR ESTIMATING THE CONSTANT S'
     1,'URVIVAL AND RECOVERY RATE IS NOT CONVERGING'//)
      WRITE (6,320) ITER
  320 FORMAT (1X,'AFTER',I6,2X,'ITERATIONS THE FOLLOWING RESULTS WERE O'
     1,'BTAINED:'/)
      PHI=PHI*100.
      WRITE (6,325) PHI
  325 FORMAT (' ',' S   = ',F10.5/)
      LAMBDA=LAMBDA*100.
      WRITE (6,330)
  330 FORMAT (1X,'THE INVERSE OF THE INFORMATION MATRIX (VARIANCE-COVAR'
     1,'IANCE) WAS:'/)
      WRITE (6,335) VARF,COVFS
      WRITE (6,335) COVFS,INVERS(2,2)
  335 FORMAT (' ',5X,2F20.14/)
  340 IDF=DEGF
      CALL CHI1(IDF,CHIADD)
      WRITE (6,345) ITER
  345 FORMAT (' ',I5,'  ITERATIONS')
      CALL DECIDE(THCHI)
      RETURN
      END
      SUBROUTINE MOD3W1
      WRITE (6,15)
   15 FORMAT ('0','MODEL 3')
      WRITE (6,20)
   20 FORMAT (' ','********'/)
      WRITE (6,25)
   25 FORMAT (1X,'ANALYSIS ASSUMING CONSTANT SURVIVAL AND RECOVERY RATE'
     1,'S')
      WRITE (6,30)
   30 FORMAT (1X,'(A GENERALIZATION AND EXTENSION OF THE MODELS DEVELOP'
     1,'ED BY CHAPMAN AND ROBSON (1960. BIOMETRICS) AND ')
      WRITE (6,35)
   35 FORMAT (1X,' HALDANE (1955. PROC. XI INT. ORN. CONGR.) -- SEE BOT'
     1,'TOM OF PAGE 245 OF BOOK BY SEBER))'/)
      WRITE (6,40)
   40 FORMAT (' ','SPECIFICALLY, THE MODEL STRUCTURE IS: '/)
      WRITE (6,45)
   45 FORMAT (' ','    F       SF       SSF       SSSF       SSSSF')
      WRITE (6,50)
   50 FORMAT (' ','             F        SF        SSF        SSSF')
      WRITE (6,55)
   55 FORMAT (' ','                       F         SF         SSF')
      WRITE (6,60)
   60 FORMAT (' ','                                  F          SF'///)
      RETURN
      END
