      SUBROUTINE MODLD2(TITLE,DCC)
      IMPLICIT DOUBLE PRECISION (A-Z)
      INTEGER IDIM
      PARAMETER (IDIM=50)
      INTEGER NITER,I,J,K,L,S,MAXITR,YR(IDIM),DCC,IPRNT
      CHARACTER*128 TITLE
      DIMENSION THETA(3),THETA1(3),H(3,3),G(3),DLT(IDIM)
      COMMON/BLK1/M(IDIM),N(IDIM),NY(IDIM),B(IDIM),R(IDIM),RY(IDIM),
     ,            NS(IDIM),SY(IDIM),YR
      COMMON /BLK3/S,FIMA,FIMD,FYMA,FYMD,EPSLON,MAXITR,IPRNT
      COMMON /BLK2/CA(IDIM,2,2),CM(IDIM,2,2),CN(IDIM,2,2),Z(IDIM,2,2)
      COMMON /BLK4/XAB(IDIM),XYB(IDIM),PB(IDIM),XA(IDIM),XY(IDIM),P
C
C      THIS SUBROUTINE COMPUTES ESTIMATES OF ADULT & YOUNG
C      SURVIVAL RATE (PHI & PHI') AND CAPTURE Probability ASSUMING
C      CONSTANT SURVIVAL RATES AND CAPTURE Probability.
C
      DCC=0
      WRITE(7,10)CHAR(12),TITLE
   10 FORMAT(A1,A128/
     //' Model D2 - Constant survival rate per unit time, ',
     2'constant capture probability per unit time.'/
     3/'  Parameters: PHI = Adult survival rate per unit time'/
     5 12X,'  PHI'' = Young survival rate per unit time'/
     6 14X,'p      = Adult capture probability per unit time'/)
C
C               FIRST, COMPUTE STARTING VALUE OF THETA
C
      XA(S)=.1D1
      DLT(S)=.0D0
      P=.0D0
      DO 20 J=2,S-1
   20   IF(R(J)*(M(J)+B(J)).GT..0D0)P=P+M(J)/(M(J)+NS(J)*B(J)/R(J))
      P=P/DBLE(FLOAT(S-2))
      IF(FIMA.LT..01)FIMA=.01
      IF(FIMA.GT..99)FIMA=.99
      IF(FYMA.LT..01)FYMA=.01
      IF(FYMA.GT..99)FYMA=.99
      THETA(1)=FIMA
      THETA(2)=P
      THETA(3)=FYMA
      FI=FIMA
      FY=FYMA
C
C             WRITE STARTING VALUES
C
      WRITE(7,30)THETA(1),THETA(3),THETA(2)
   30 FORMAT(/' Starting values of PHI (adult, young), and p:',3F10.4)
      NITER=0
C
C          START OF ITERATIVE PROCEDURE
C              SET SURV & CAPT. PROBS = THETA
C
   40 FI=THETA(1)
      P =THETA(2)
      Q =.1D1-P
      FY=THETA(3)
      DO 50 I=S-1,1,-1
        TMP=.1D1-Q*XA(I+1)
        XA(I)=.1D1-FI*TMP
        XY(I)=.1D1-FY*TMP
   50   DLT(I)=DLT(I+1)*FI*Q+(.1D1-XA(I))
C
C          COMPUTE GRADIANT VECTOR (G) &  INVERSE OF VAR-COV MATRIX(H)
C
      DO 60 I=1,3
        G(I)=.0D0
        DO 60 J=1,3
   60     H(J,I)=.0D0
      DO 70 I=1,S-1
C       G(1)=G(1)+R(I)+RY(I)+B(I)-DLT(I)*
C    (        ((NS(I)-R(I))/XA(I)+FY/FI*(SY(I)-RY(I))/XY(I))
        G(1)=G(1)+M(I+1)+B(I+1)-DLT(I)*
     (        ((NS(I)-R(I))/XA(I)+FY/FI*(SY(I)-RY(I))/XY(I))
        G(2)=G(2)+NS(I)-(NS(I)-R(I))/XA(I)
        G(3)=G(3)+SY(I)-(SY(I)-RY(I))/XY(I)
        X1A=.1D1-XA(I)
        X1Y=.1D1-XY(I)
        QXAI=.1D1-Q*XA(I)
        H(1,1)=H(1,1)+DLT(I)*DLT(I)/X1A*
     (    (NS(I)/XA(I)+M(I)*Q/QXAI)+
     +     Q*Q*FY*FY*DLT(I+1)*DLT(I+1)*SY(I)/(XY(I)*X1Y)
        H(2,2)=H(2,2)+((DLT(I)/Q-X1A/(P*Q))**2*
     *     (NS(I)/XA(I)+FY/FI*SY(I)/XY(I))+
     +     DLT(I)*DLT(I)*M(I)/(Q*QXAI))/X1A
        H(3,3)=H(3,3)+SY(I)*X1Y/XY(I)
        H(1,2)=H(1,2)+((DLT(I)-X1A/P)*(NS(I)/XA(I)*
     *     DLT(I)/Q+FY*DLT(I+1)*SY(I)/XY(I))+
     +     M(I)*DLT(I)*DLT(I)/QXAI)/X1A
        H(2,3)=H(2,3)+SY(I)/XY(I)*(DLT(I)-X1A/P)
   70   H(1,3)=H(1,3)+SY(I)*DLT(I+1)/XY(I)
      G(1)=G(1)/FI
      G(2)=G(2)/P/Q
      G(3)=G(3)/FY
      H(1,1)=H(1,1)/FI/FI
      H(3,3)=H(3,3)/FY/FY
      H(1,2)=-H(1,2)/FI
      H(2,1)=H(1,2)
      H(1,3)=H(1,3)*Q/FI
      H(3,1)=H(1,3)
      H(2,3)=-H(2,3)/Q/FI
      H(3,2)=H(2,3)
C
C                  COMPUTE NEW THETA = THETA + INV(H)*G
C
      CALL MATINV(3,H,3)
      CALL MMULT(H,G,THETA1,3,3,1,3,3,1)
      CALL MATADD(THETA1,THETA,THETA1,3)
C
C                 CHECK FOR CONVERGENCE
C
      DO 80 I=1,3
        IF(DABS(THETA(I)-THETA1(I)).GT.EPSLON)GO TO 90
   80   CONTINUE
      GO TO 140
C
C     NOT CONVERGED... INCREASE # ITERATIONS & SET THETA=NEW THETA
C
   90 NITER=NITER+1
      IF(NITER.GT.MAXITR)GO TO 110
      IF(THETA1(1).LT..0D0.OR.THETA1(3).LT..0D0)GO TO 110
      DO 100 I=1,3
        IF(DABS(THETA1(I)).GT..1D2)GO TO 110
  100   THETA(I)=THETA1(I)
      GO TO 40
C
C                 CONVERGENCE FAILURE... OUTPUT MESSAGE & QUIT
C
  110 WRITE(7,120)NITER,(THETA(I),I=1,3)
  120 FORMAT(/' Convergence failure after ',I6,' iterations'/
     /' THETA = '/(10F12.6))
      WRITE(7,130)(THETA1(I),I=1,3)
  130 FORMAT(/' THETA1 = '/(10F12.6))
      RETURN
C
C           CONVERGED... OUTPUT RESULTS
C
  140 DCC=1
      SE=XSQRT(H(1,1))
      C1=FI-SE*.196D1
      C2=FI+SE*.196D1
      WRITE(7,150)NITER,FI,H(1,1),SE,C1,C2,H(1,3)
  150 FORMAT(/' Final values after ',I6,' iterations'//
     /' Parameter  Estimate  Variance  Std.error',
     ,'   95% confidence interval     Covariance W/ PHI',3X,
     ,'Covariance W/ PHI',1H'/1X,110('-')/
     /5X,'PHI',4X,F8.4,F10.4,F11.4,F10.4,' - ',F8.4,31X,F13.8/)
      SE=XSQRT(H(3,3))
      C1=FY-SE*.196D1
      C2=FY+SE*.196D1
      WRITE(7,160)FY,H(3,3),SE,C1,C2,H(1,3)
  160 FORMAT(5X,4HPHI',3X,F8.4,F10.4,F11.4,F10.4,' - ',F8.4,10X,F13.8/)
      SE=XSQRT(H(2,2))
      C1=P-SE*.196D1
      C2=P+SE*.196D1
      WRITE(7,170)P,H(2,2),SE,C1,C2
  170 FORMAT('     p',6X,F8.4,F10.4,F11.4,F10.4,' - ',
     ,  F8.4/1X,62('-'))
      CALL D2MNB(TITLE,FI,FY,P,XA,XY,DLT,H)
      FIMD=FI
      FYMD=FY
      RETURN
      END
