C This program was originally written by P. McCullagh (Dept. of
C Statistics, University of Chicago)
C It was adapted to Unix Fortran by X. Zhou (Dept. of Health Care Policy,
C Harvard Medical School)
C***********************************************************************
C*** PROGRAM TITLE
C      MASTER PLUM
C*********************************************************************
C   BEFORE LOADING CHECK ALL LINES PRECEDED BY C***
C**********************************************************************
C
      DOUBLE PRECISION RDATA,LOGLIK,Q,XBAR
      COMMON RDATA(5000),LOGLIK,XBAR(50),PROP(20),Q(400)
      REAL DATA(10000),NOTE,LINKNM
      EQUIVALENCE (RDATA(1),DATA(1))
      INTEGER R,ODD,PRINT
      COMMON/SPEC/T,K,KMAX,LINK,NV,LINKM(8),ODD,TITLE(80)
      LOGICAL BATCH,LREAD,ECHO,LFORM,LTIT
      INTEGER OUT6,OUTCH
      COMMON/PERIPH/IN5,OUT6,INCH,OUTCH,BATCH
      COMMON/ENVIR/NCHW,NCHL,NCHLM3,NCHLO
      COMMON/NO/NO(20),NOTE(80),NOBS
      COMMON/LOGIC/MAXIND,MAXCAT,LREAD,ECHO,IFL,LFORM
      COMMON/CONTRL/INDSSP,INDIND,MAXDAT,MAXCYC,PRINT,LTIT
      COMMON/LEFT/FRMT(20)
      LOGICAL START,RET
      INTEGER CHANL
      DATA BLANK/1H /
      DATA START/.FALSE./
C
      INCH=IN5
      OUTCH=OUT6
      IF(START)GOTO 20
  10  CALL INITIL
      START=.TRUE.
C*** PROMPT FOR NEW COMMAND-NOT NECESARY ON SOME SYSTEMS
  20  WRITE(*,501)
C***   NON STARNDARD FORMAT FOR PROMPT
 501  FORMAT(1X,1H?,'$')
      READ(*,502) DIR,(NOTE(I),I=1,NCHLM3)
 502  FORMAT(A3,80A1)
      IF(BATCH) WRITE(OUTCH,503) DIR,(NOTE(I),I=1,NCHLM3)
 503  FORMAT(1X,A3,80A1)
      CALL SYCODE(IFAIL,DIR,IDIR)
      IF (IFAIL.EQ.0) GO TO 30
      CALL ERROR(8)
      IF(BATCH)STOP
      GO TO 20
C   DECISION POINTER TO CORRECT SUBROUTINE
   30 GO TO (40,50,60,70,80,90,100,110,120,130,140,150,160,
     2170,180,190,200,210,220,230,240,250,260,270,280),IDIR
C  FORMAT DECLARATION
   40 N=NCHL/NCHW
      LFORM=.TRUE.
C*** PROMPT FOR FORMAT IF NECESSARY
      IF(INCH.EQ.IN5)WRITE(OUT6,508)
C*** NON-STANDARD FORMAT FOR PROMPT
  508 FORMAT(1X,1H*,'$')
      READ(*,504) (FRMT(I),I=1,N)
C*** FORMAT DEPENDS ON NUMBER OF CHARS/REAL WORD 20 A NCHW
 504  FORMAT(20A4)
      IF(BATCH.OR.ECHO) WRITE(OUTCH,505)(FRMT(I),I=1,N)
C*** FORMAT DEPENDS ON NUMBER OF CHARS/REAL WORD 20A NCHW
 505  FORMAT(1X,6HFORMAT,20A4)
      IF(BATCH.OR.ECHO) WRITE(OUTCH,506)
 506  FORMAT(1X,26HFORMAT SYNTAX NOT CHECKED )
      GO TO 20
C
C   INPUT AND OUTPUT CHANNEL DECLARATION
  50  CALL INTERP(NOTE,NCHLM3,NO,0,2,NOBS,IFL)
      IF(IFL.EQ.0) GO TO 55
      CALL ERROR(5)
      IF(BATCH)STOP
      GO TO 20
  55  INCH=IN5
      OUTCH=OUT6
      IF(NOBS.GE.1) INCH=NO(1)
      IF(NOBS.EQ.2) OUTCH=NO(2)
      IF(ECHO) WRITE(OUTCH,507) INCH,OUTCH
 507  FORMAT(1X,18HDATA INPUT CHANNEL ,I3/
     21X,19HDATA OUTPUT CHANNEL ,I3)
      GO TO 20
C
  60  CALL LINKS
      GO TO 20
C
  70  CALL READ
      GO TO 20
C
  80  CALL LOC
      GO TO 20
C
  90  CALL SCALE
      GO TO 20
C
C  FIT DIRECTIVE
 100  IF(LREAD) GO TO 103
      CALL ERROR(6)
      IF(BATCH)STOP
      GO TO 20
 103  CALL INTERP(NOTE,NCHLM3,NO,0,1,NOBS,IFL)
      IF(IFL.EQ.0) GO TO 106
      CALL ERROR(5)
      IF(BATCH) STOP
      GO TO 20
 106  PRINT=1
      IF(NOBS.GT.0)PRINT=NO(1)
      IF(K.GE.2)GOTO 108
      CALL ERROR(7)
      IF(BATCH)STOP
      GOTO 20
  108 CALL FIT
      GO TO 20
 110  CALL INITIL
      GOTO 20
 120  STOP
C
 130  ECHO=.NOT.ECHO
      GO TO 20
 140  BATCH=.NOT.BATCH
      GO TO 20
 150  CALL HELP
      GO TO 20
 160  GO TO 20
C NOTE ABOVE IS A NOTE DIRECTIVE ONLY
C  TITLE DIRECTIVE
 170  DO 175 I=1,NCHLM3
 175  TITLE(I)=NOTE(I)
      TITLE(1)=BLANK
      TITLE(2)=BLANK
      LTIT=.TRUE.
      GO TO 20
C
C  CANCEL FORMAT DECLARATION
 180  LFORM=.FALSE.
      GO TO 20
 190  CALL FACTOR
      GO TO 20
 200  CALL XVAR
      GO TO 20
 210  IFL=1
      CALL COMBIN(1,1)
      GO TO 20
 220  CALL RESET
      GO TO 20
 230  CALL LOOK
      GO TO 20
C JACKKNIFE COMMAND
  240 IF(LREAD)GOTO 242
      CALL ERROR(6)
      IF(BATCH)STOP
      GOTO 20
  242 CALL INTERP(NOTE,NCHLM3,NO,0,1,NOBS,IFL)
      IF(IFL.EQ.0)GOTO 246
  244 CALL ERROR(5)
      IF(BATCH)STOP
      GOTO 20
  246 IF(NOBS.EQ.0)NO(1)=10000
      PRINT=-1
      IF(NO(1).GT.K)GOTO 248
      CALL ERROR(15)
      IF(BATCH)STOP
      GOTO 20
  248 FRAC=1./FLOAT(NO(1))
      CALL JACK(FRAC)
      GOTO 20
  250 CALL INTERP(NOTE,NCHLM3,NO,0,1,NOBS,IFL)
      IF(IFL.NE.0)GOTO 244
      MAXCYC=15
      IF(NOBS.EQ.1)MAXCYC=NO(1)
      GOTO 20
  260 RET=.FALSE.
C DUMP OUTPUT
      CALL INTERP(NOTE,NCHLM3,NO,0,1,NOBS,IFL)
      IF(IFL.NE.0)GOTO 244
      IF(NOBS.EQ.1)GOTO 265
      CALL ERROR(12)
      IF(BATCH)STOP
      GOTO 20
  265 CHANL=NO(1)
      CALL DUMP(RET,CHANL)
      GOTO 20
  270 RET=.TRUE.
      CALL INTERP(NOTE,NCHLM3,NO,0,1,NOBS,IFL)
      IF(IFL.NE.0)GOTO 244
      IF(NOBS.EQ.1)GOTO 265
      CALL ERROR(12)
      IF(BATCH)STOP
      GOTO 20
  280 CALL INTERP(NOTE,NCHLM3,NO,0,1,NOBS,IFL)
      IF(IFL.NE.0)GOTO 244
      IF(NOBS.EQ.1)GOTO 285
      CALL ERROR(12)
      IF(BATCH) STOP
      GOTO 20
  285 CHANL=NO(1)
      REWIND CHANL
      GOTO 20
      END
C###############################################################
      BLOCK DATA
      REAL LINKNM
      LOGICAL LTIT,LFORM,LREAD,ECHO,LSCLP,BATCH
      INTEGER OUT6,OUTCH,T,ODD,PRINT
      COMMON/PERIPH/IN5,OUT6,INCH,OUTCH,BATCH
      COMMON/SPEC/T,K,KMAX,LINK,NV,LINKNM(8),ODD,TITLE(80)
      DOUBLE PRECISION BETA
      COMMON/PARAM/BETA(100),INDEX1(50),INDEX2(50),NLOCP,NSCLP,LSCLP
      COMMON/LOGIC/MAXIND,MAXCAT,LREAD,ECHO,IFL,LFORM
      COMMON/CONTRL/INDSSP,INDIND,MAXDAT,MAXCYC,PRINT,LTIT
      DOUBLE PRECISION ZERO,ONE,HALF,TWO,RT2PI,PINV,ETA,TOL
      COMMON/CONS/ZERO,ONE,HALF,TWO,RT2PI,PINV,ETA,TOL
      COMMON/ENVIR/NCHW,NCHL,NCHLM3,NCHLO
C*** NCHL AND NCHW ARE NUMBERS OF CHARACTERS PER INPUT LINE AND
C*** PER REAL WORD RESPECTIVELY.
      DATA NCHL/72/,NCHW/4/
      DATA ZERO/0.D0/,ONE/1.D0/,HALF/0.5D0/,TWO/2.D0/
      DATA RT2PI/0.398942280D0/,PINV/0.3183098862D0/
C*** ETA SHOULD BE SET TO THE SMALLEST POSITIVE NUMBER TO BE
C*** TAKEN AS NON-ZERO. TOL GOVERNS ACCURACY OF PARAMETER ESTIMATES.
      DATA ETA/1.0D-10/,TOL/1.D-4/
      DATA LINKNM(1)/4HLOGI/,LINKNM(2)/4HT   /,LINKNM(3)/4HPROB/,
     * LINKNM(4)/4HIT  /,LINKNM(5)/4HCAUC/,LINKNM(6)/4HHIT /,
     * LINKNM(7)/4HC-LO/,LINKNM(8)/4HGLOG/
      DATA MAXIND,MAXCAT,MAXDAT,MAXCYC/50,20,10000,15/
C*** SET DEFAULT INPUT AND OUTPUT CHANNELS
      DATA IN5,OUT6/5,6/
C*** BATCH VS. INTERACTIVE
      DATA BATCH/.FALSE./
      DATA ECHO/.TRUE./
C*** SET DEFAULT LINE LENGTH FOR OUTPUT
      DATA NCHLO/72/
      END
C##############################################################
      SUBROUTINE INITIL
C SET DEFAULT FOR NEW JOB OR NEW DATA SET
      REAL LINKNM
      DOUBLE PRECISION RDATA,LOGLIK,Q,XBAR
      LOGICAL LTIT,LFORM,LREAD,ECHO,LSCLP,BATCH
      INTEGER OUT6,OUTCH,T,ODD,PRINT
      COMMON RDATA(5000),LOGLIK,XBAR(50),PROP(20),Q(400)
      COMMON/PERIPH/IN5,OUT6,INCH,OUTCH,BATCH
      COMMON/SPEC/T,K,KMAX,LINK,NV,LINKNM(8),ODD,TITLE(80)
      DOUBLE PRECISION BETA
      COMMON/PARAM/BETA(100),INDEX1(50),INDEX2(50),NLOCP,NSCLP,LSCLP
      COMMON/LOGIC/MAXIND,MAXCAT,LREAD,ECHO,IFL,LFORM
      COMMON/CONTRL/INDSSP,INDIND,MAXDAT,MAXCYC,PRINT,LTIT
      DOUBLE PRECISION ZERO,ONE,HALF,TWO,RT2PI,PINV,ETA,TOL
      COMMON/CONS/ZERO,ONE,HALF,TWO,RT2PI,PINV,ETA,TOL
      COMMON/ENVIR/NCHW,NCHL,NCHLM3,NCHLO
C
      MAXCYC=15
      LINK=1
      NCHLM3=NCHL-3
      LTIT=.FALSE.
      LFORM=.FALSE.
      LREAD=.FALSE.
      LSCLP=.FALSE.
      DO 10 I=1,MAXIND
  10  XBAR(I)=ZERO
      DO 20 I=1,100
  20  BETA(I)=ZERO
      DO 30 I=1,MAXCAT
  30  PROP(I)=0.0
      RETURN
      END
C  ##########################################
      SUBROUTINE SYCODE(IFAIL,DIR,IDIR)
C FIND DIRECTORY CODE
      REAL SYS(25),DIR
      DATA SYS(1)/3HFOR/,SYS(2)/3HI/O/,SYS(3)/3HLIN/,
     2SYS(4)/3HDAT/,SYS(5)/3HLOC/,SYS(6)/3HSCA/,SYS(7)/
     33HFIT/,SYS(8)/3HEND/,SYS(9)/3HSTO/,SYS(10)/3HECH/,
     4SYS(11)/3HBAT/,SYS(12)/3HHEL/,SYS(13)/3HNOT/,SYS(
     514)/3HTIT/,SYS(15)/3HFRE/,SYS(16)/3HFAC/
      DATA SYS(17)/3HXVA/,SYS(18)/3HCOM/,SYS(19)/3HRES/
      DATA SYS(20)/3HLOO/,SYS(21)/3HJAC/,SYS(22)/3HCYC/
      DATA SYS(23)/3HDUM/,SYS(24)/3HRET/,SYS(25)/3HREW/
      DATA NDIR/25/
C
      IFAIL=0
      DO 10 J=1,NDIR
      IDIR=J
      IF(DIR.EQ.SYS(J)) RETURN
  10  CONTINUE
      IFAIL=1
      RETURN
      END
C##############################################################
      SUBROUTINE DUMP(RET,CHANL)
      INTEGER CHANL
      LOGICAL RET
      REAL LINKNM
      DOUBLE PRECISION RDATA,LOGLIK,Q,XBAR
      LOGICAL LTIT,LFORM,LREAD,ECHO,LSCLP,BATCH
      INTEGER OUT6,OUTCH,T,ODD,PRINT
      COMMON RDATA(5000),LOGLIK,XBAR(50),PROP(20),Q(400)
      COMMON/PERIPH/IN5,OUT6,INCH,OUTCH,BATCH
      COMMON/SPEC/T,K,KMAX,LINK,NV,LINKNM(8),ODD,TITLE(80)
      DOUBLE PRECISION BETA
      COMMON/PARAM/BETA(100),INDEX1(50),INDEX2(50),NLOCP,NSCLP,LSCLP
      COMMON/LOGIC/MAXIND,MAXCAT,LREAD,ECHO,IFL,LFORM
      COMMON/CONTRL/INDSSP,INDIND,MAXDAT,MAXCYC,PRINT,LTIT
      DOUBLE PRECISION ZERO,ONE,HALF,TWO,RT2PI,PINV,ETA,TOL
      COMMON/CONS/ZERO,ONE,HALF,TWO,RT2PI,PINV,ETA,TOL
      COMMON/ENVIR/NCHW,NCHL,NCHLM3,NCHLO
      COMMON/LEFT/FRMT(20)
      IF(RET)GOTO 20
      IF(LREAD)GOTO 10
      CALL ERROR(6)
      IF(BATCH)STOP
      RETURN
   10 INDSSP=(INDIND+T*NV)/2+1
      WRITE(CHANL)INCH,OUTCH,BATCH,ECHO,LINK,LREAD,LTIT,LFORM,ODD
      IF(LTIT)WRITE(CHANL)TITLE
      IF(LFORM)WRITE(CHANL)FRMT
      WRITE(CHANL)XBAR,PROP
      WRITE(CHANL)INDSSP,(RDATA(J),J=1,INDSSP)
      WRITE(CHANL)T,K,KMAX,NV,LSCLP,NSCLP,NLOCP, INDIND
      WRITE(CHANL)BETA,INDEX1,INDEX2
      RETURN
   20 READ(CHANL)INCH,OUTCH,BATCH,ECHO,LINK,LREAD,LTIT,LFORM,ODD
      IF(LTIT)READ(CHANL)TITLE
      IF(LFORM)READ(CHANL)FRMT
      READ(CHANL)XBAR,PROP
      READ(CHANL)INDSSP,(RDATA(J),J=1,INDSSP)
      READ(CHANL)T,K,KMAX,NV,LSCLP,NSCLP,NLOCP, INDIND
      READ(CHANL)BETA,INDEX1,INDEX2
      RETURN
      END
      SUBROUTINE LINKS
C   LINK DECLARATION
      LOGICAL BATCH,ECHO,LTIT,LREAD,LFORM
      REAL NOTE,LINKNM
      INTEGER T,ODD,PRINT,OUT6,OUTCH
      COMMON/SPEC/T,K,KMAX,LINK,NV,LINKNM(8),ODD,TITLE(80)
      COMMON/PERIPH/IN5,OUT6,INCH,OUTCH,BATCH
      COMMON/ENVIR/NCHW,NCHL,NCHLM3,NCHLO
      COMMON/NO/NO(20),NOTE(80),NOBS
      COMMON/LOGIC/MAXIND,MAXCAT,LREAD,ECHO,IFL,LFORM
C
      CALL INTERP(NOTE,NCHLM3,NO,0,1,NOBS,IFL)
      IF(IFL.EQ.0) GO TO 10
      CALL ERROR(5)
      IF(BATCH)STOP
      RETURN
  10  IF(NOBS.EQ.1) LINK=NO(1)
      IF(NOBS.EQ.0) LINK=1
      IF(LINK.LE.4.AND.LINK.NE.0) GO TO 20
      CALL ERROR(10)
      IF(BATCH) STOP
      LINK=1
      RETURN
  20  L1=2*LINK-1
      L2=L1+1
      IF(ECHO) WRITE(OUTCH,501)LINK,LINKNM(L1),LINKNM(L2)
 501  FORMAT(1X ,5HLINK ,I2,3H = ,2A4)
      IF(LREAD)CALL RESET
      RETURN
      END
C###########################################################
      SUBROUTINE READ
C  DATA STRUCTURE DECLARATION AND READ, CUMULATIVE PROPORTIONS
C   AND CHECK FOR ZERO COLUMNS
      DOUBLE PRECISION RDATA,LOGLIK,Q,XBAR
      REAL WORK(800),NOTE,LINKNM
      COMMON RDATA(5000),LOGLIK,XBAR(50),PROP(20),Q(400)
      EQUIVALENCE(WORK(1),Q(1))
      REAL DATA(10000)
      EQUIVALENCE(RDATA(1),DATA(1))
      INTEGER T,ODD,PRINT
      CHARACTER*12 infname,lstring
      LOGICAL*1 ex
      COMMON/SPEC/T,K,KMAX,LINK,NV,LINKNM(8),ODD,TITLE(80)
      DOUBLE PRECISION BETA
      COMMON/PARAM/BETA(100),INDEX1(50),INDEX2(50),NLOCP,NSCLP,LSCL
      LOGICAL BATCH,ECHO,LSCLP,LTIT,LREAD,LFORM
      COMMON/CONTRL/INDSSP,INDIND,MAXDAT,MAXCYC,PRINT,LTIT
      COMMON/MODEL/NLOC,NSCL,NULL,ALAIS(100)
      COMMON/LOGIC/MAXIND,MAXCAT,LREAD,ECHO,IFL,LFORM
      INTEGER OUT6,OUTCH
      COMMON/NO/NO(20),NOTE(80),NOBS
      COMMON/PERIPH/IN5,OUT6,INCH,OUTCH,BATCH
      COMMON/ENVIR/NCHW,NCHL,NCHLM3,NCHLO
      COMMON/LEFT/FRMT(20)
C
      IF(.NOT.LREAD) GO TO 10
      CALL ERROR(11)
      IF(BATCH)STOP
      RETURN
  10  CALL INTERP(NOTE,NCHLM3,NO,0,3,NOBS,IFL)
      IF(IFL.EQ.0) GO TO 30
      CALL ERROR(5)
      IF(BATCH)STOP
      RETURN
   30 IF(NOBS.GE.2)GOTO 35
      CALL ERROR(12)
      IF(BATCH)STOP
      RETURN
   35 IF(NOBS.EQ.2)NO(3)=0
      IFAIL=0
      IF(NO(2).LE.1.OR.NO(2).GT.MAXCAT)IFAIL=1
      IF(NO(3).GT.MAXIND)IFAIL=1
      NTOT=NO(1)*(2*NO(2)+6+NO(3))+NO(2)*(NO(2)+1)+2
      IF(NTOT.LE.MAXDAT) GO TO 40
      CALL ERROR(14)
      IF(BATCH)STOP
      RETURN
  40  IF(IFAIL.EQ.0) GO TO 50
      CALL ERROR(15)
      IF(BATCH)STOP
      RETURN
  50  T=NO(1)
      K=NO(2)
      K1=K-1
      KMAX=K
      NV=NO(3)
      ODD=0
      IF(((T*K+T)/2)*2.NE.T*(K+1)) ODD=1
      NSCLP=0
      NSCL=0
      NLOCP=0
      INDIND=2*T*K+6*T+ODD+1
      IF(ECHO)WRITE(OUTCH,501)T,K,NV
 501  FORMAT(1X,8HUNITS = ,I2,14H CATEGORIES = ,I2/
     213H COVARIATES =,I3)
      LREAD=.TRUE.
      TOT=0.0
      L2=T*K
C     Give an option of either input data from file or keyboard
      write(*,300)
 300  format(1h1)
 301  write(*,*)'Do you want to input data from a file ?'
      write(*,*)' (Y/N)'
      read(5,302) lstring
 302  format(A1)
      if (lstring .eq. 'Y' .or. lstring .eq. 'y') then
 303  write(*,*)'Specify the name of input file:'
      read(5,304)infname
 304  format(a)
      ex=.false.
      inquire(file=infname,exist=ex)
      if(ex .eqv. .false.)then
         write(*,305)infname
 305     format(1x,'The input file (',a,') can not be located in',
     c              'your directory!')
         go to 303
         endif
         nin=2
         open(nin,file=infname,status='old')
      DO 360 I=1,T
      IF(NV.EQ.0) GO TO 306
      READ(nin,*)(DATA(J),J=I,L2,T),(WORK(J),J=1,NV)
 700  DO 800 J=1,NV
      M2=INDIND+I-1+(J-1)*T
      DATA(M2)=WORK(J)
      XBAR(J)=XBAR(J)+DBLE(WORK(J))
 800   continue
 900   VN=0.0
      DO 361 J=I,L2,T
      VN=VN+DATA(J)
 361  continue
      DATA(L2+I)=VN
      TOT=TOT+VN
      DO 362 J1=1,K
      PROP(J1)=PROP(J1)+DATA(I+(J1-1)*T)
 362  continue
      GO TO 360
 306  READ(nin,*)(DATA(J1),J1=I,L2,T)
      GO TO 900
 360  continue                      
      go to 1000 
      else
      DO 135 I=1,T
C*** PROMPT FOR NEW COMMAND - NOT NECESSARY ON SOME SYSTEMS
      IF(.NOT.BATCH.AND.IN5.EQ.INCH)WRITE(OUT6,502)
C*** NON-STANDARD FORMAT FOR PROMPT
 502  FORMAT(1X,1H*,'$')
      IF(LFORM) GO TO 120
      IF(NV.EQ.0) GO TO 110
C*** FREE FORMAT DATA INPUT
      READ(INCH,*)(DATA(J),J=I,L2,T),(WORK(J),J=1,NV)
      IF(BATCH)WRITE(OUTCH,504)(DATA(J),J=I,L2,T),(WORK(J),J=1,NV)
  504 FORMAT(1X,17F7.2)
  70  DO 80 J=1,NV
      M2=INDIND+I-1+(J-1)*T
      DATA(M2)=WORK(J)
  80  XBAR(J)=XBAR(J)+DBLE(WORK(J))
  90  VN=0.0
      DO 100 J=I,L2,T
 100  VN=VN+DATA(J)
      DATA(L2+I)=VN
      TOT=TOT+VN
      DO 105 J1=1,K
 105  PROP(J1)=PROP(J1)+DATA(I+(J1-1)*T)
      GO TO 135
C*** FREEE FORMAT DATA INPUT
 110  READ(INCH,*)(DATA(J1),J1=I,L2,T)
      IF(BATCH)WRITE(OUTCH,505)(DATA(J),J=I,L2,T)
  505 FORMAT(1X,20F6.1)
      GO TO 90
C    FORMATTED DATA INPUT
 120  IF(NV.EQ.0) GO TO 130
      READ(INCH,FRMT)(DATA(J),J=I,L2,T),(WORK(J),J=1,NV)
      IF(BATCH)WRITE(OUTCH,504)(DATA(J),J=I,L2,T),(WORK(J),J=1,NV)
      GO TO 70
C    FORMATTED DATA INPUT
 130  READ(INCH,FRMT)(DATA(J),J=I,L2,T)
      IF(BATCH)WRITE(OUTCH,505)(DATA(J),J=I,L2,T)
      GO TO 90
  135 CONTINUE
      IF(LFORM.AND.ECHO.AND..NOT.BATCH)WRITE(OUTCH,503)
  503 FORMAT(1X,50HFORMAT REMAINS IN FORCE. CANCEL OR CHANGE IF REQD.)
C
      endif
 1000 IF(NV.EQ.0) GO TO 150
      VN=FLOAT(T)
      DO 140 J=1,NV
 140  XBAR(J)=XBAR(J)/DBLE(VN)
 150  VN=0.0
      DO 160 J=1,K
      VN=VN+PROP(J)
 160  PROP(J)=VN/TOT
C  CHECK FOR ZERO COLUMNS
 170  N1=1
      IF(PROP(1).EQ.0.) GO TO 190
      K1=K-1
      DO 180 J=1,K1
      N1=J
      IF(PROP(J).EQ.PROP(J+1)) GO TO 190
 180  CONTINUE
      CALL RESET
      RETURN
 190  N2=N1+1
      IFL=0
      IF(K.GT.2)GOTO 200
      CALL ERROR(7)
      IF(BATCH)STOP
      CALL COMBIN(N1,N2)
      RETURN
  200 CALL COMBIN(N1,N2)
      GOTO 170
      END
C
C####################################################
      SUBROUTINE LOC
C  LOCATION MODEL SPECIFICATION
      LOGICAL LREAD,ECHO,LFORM,BATCH
      REAL NOTE
      INTEGER OUT6,OUTCH
      COMMON/SPEC/T,K,KMAX,LINK,NV,LINKNM(8),ODD,TITLE(80)
      COMMON/MODEL/NLOC,NSCL,NULL,ALIAS(100)
      DOUBLE PRECISION BETA
      COMMON/PARAM/BETA(100),INDEX1(50),INDEX2(50),NLOCP,NSCLP,LSCLP
      COMMON/PERIPH/IN5,OUT6,INCH,OUTCH,BATCH
      COMMON/ENVIR/NCHW,NCHL,NCHLM3,NCHLO
      COMMON/NO/NO(20),NOTE(80),NOBS
      COMMON/LOGIC/MAXIND,MAXCAT,LREAD,ECHO,IFL,LFORM
      COMMON/LEFT/FRMT(20)
C
      IF(LREAD) GO TO 10
      CALL ERROR(6)
      IF(BATCH)STOP
      RETURN
   10 IF(NV.NE.0)GOTO 15
      CALL ERROR(9)
      IF(BATCH)STOP
      RETURN
   15 CALL INTERP(NOTE,NCHLM3,NO,0,NV,NLOCP,IFL)
      IF(IFL.EQ.0) GO TO 20
      CALL ERROR(5)
      IF(BATCH) STOP
      NLOCP=0
      RETURN
  20  IF(NLOCP.GT.0) GO TO 30
      IF(ECHO)WRITE(OUTCH,501)
 501  FORMAT(1X,19HNULL LOCATION MODEL)
      RETURN
   30 NTEMP=NLOCP
      DO 40 J=1,NTEMP
      IF(NO(J).LE.NV.AND.NO(J).GT.0) GO TO 40
      CALL ERROR(18)
      NLOCP=0
      IF(BATCH)STOP
      IF(ECHO)WRITE(OUTCH,501)
      RETURN
  40  INDEX1(J)=NO(J)
      IF(ECHO)WRITE(OUTCH,502)(INDEX1(J),J=1,NLOCP)
 502  FORMAT(1X,23HVARS. IN LOCATION MODEL,20I3)
      RETURN
      END
C  ####################################################
      SUBROUTINE SCALE
C  SCALE MODEL SPECIFICATION
      REAL NOTE
      LOGICAL LREAD,LSCLP,ECHO,LFORM,BATCH
      INTEGER OUT6,OUTCH
      COMMON/SPEC/T,K,KMAX,LINK,NV,LINKNM(8),ODD,TITLE(80)
      COMMON/MODEL/NLOC,NSCL,NULL,ALIAS(100)
      COMMON/PERIPH/IN5,OUT6,INCH,OUTCH,BATCH
      COMMON/ENVIR/NCHW,NCHL,NCHLM3,NCHLO
      COMMON/NO/NO(20),NOTE(80),NOBS
      COMMON/LOGIC/MAXIND,MAXCAT,LREAD,ECHO,IFL,LFORM
      DOUBLE PRECISION BETA
      COMMON/PARAM/BETA(100),INDEX1(50),INDEX2(50),NLOCP,NSCLP,LSCLP
      COMMON/LEFT/FRMT(20)
C
      IF(LREAD) GO TO 10
      CALL ERROR(6)
      IF(BATCH)STOP
      RETURN
  10  IF(NV.NE.0) GO TO 20
      CALL ERROR(9)
      RETURN
  20  CALL INTERP(NOTE,NCHLM3,NO,0,NV,NSCLP,IFL)
      IF(IFL.EQ.0) GO TO 40
      CALL ERROR(5)
      IF(BATCH)STOP
  30  LSCLP=.FALSE.
      NSCLP=0
      NSCL=0
      IF(ECHO) WRITE(OUTCH,501)
      RETURN
  40  IF(NSCLP.GT.0) GO TO 50
 501  FORMAT(1X,16HNULL SCALE MODEL)
      GO TO 30
  50  DO 60 J=1,NSCLP
      IF(NO(J).LE.NV.AND.NO(J).GT.0) GO TO 60
      CALL ERROR(18)
      IF(BATCH)STOP
      GO TO 30
  60  INDEX2(J)=NO(J)
      NSCL=NSCLP
      LSCLP=.TRUE.
      IF(ECHO)WRITE(OUTCH,502)(INDEX2(J),J=1,NSCLP)
 502  FORMAT(1X,20HVARS. IN SCALE MODEL ,20I3)
      RETURN
      END
C  ########################################################
      SUBROUTINE FACTOR
C  FACTOR DECLARATIONS
      DOUBLE PRECISION RDATA,LOGLIK,Q,XBAR
      COMMON RDATA(5000),LOGLIK,XBAR(50),PROP(20),Q(400)
      DIMENSION NO(800)
      EQUIVALENCE (NO(1),Q(1))
      REAL DATA(10000),NOTE,LINKNM
      EQUIVALENCE(RDATA(1),DATA(1))
      INTEGER T,ODD,PRINT,OUT6,OUTCH
      COMMON/SPEC/T,K,KMAX,LINK,NV,LINKNM(8),ODD,TITLE(80)
      LOGICAL LTIT,LREAD,BATCH,ECHO,LFORM
      COMMON/CONTRL/INDSSP,INDIND,MAXDAT,MAXCYC,PRINT,LTIT
      COMMON/PERIPH/IN5,OUT6,INCH,OUTCH,BATCH
      COMMON/ENVIR/NCHW,NCHL,NCHLM3,NCHLO
      COMMON/NO/NTO(20),NOTE(80),NOBS
      COMMON/LOGIC/MAXIND,MAXCAT,LREAD,ECHO,IFL,LFORM
      COMMON/LEFT/FRMT(20)
      DATA BLANK/1H /
C
      IF(LREAD) GO TO 10
      CALL ERROR(6)
      IF(BATCH) STOP
      RETURN
  10  CALL INTERP(NOTE,NCHLM3,NO,0,2,NOBS,IFL)
      IF(IFL.EQ.0) GO TO 20
      CALL ERROR(5)
      IF(BATCH)STOP
      RETURN
  20  IF(NOBS.EQ.1) GO TO 90
      IF(NOBS.EQ.2) GO TO 30
      CALL ERROR(16)
      IF(BATCH)STOP
      RETURN
  30  NLEVS=NO(1)
      IF(NLEVS.GT.1.AND.NLEVS.LE.MAXIND-NV.AND.NLEVS.LE.T)
     2 GO TO 40
      CALL ERROR(17)
      IF(BATCH)STOP
      RETURN
  40  NREP=NO(2)
      IF(NREP.NE.0) GO TO 45
      CALL ERROR(19)
      IF(BATCH)STOP
      RETURN
  45  IF((NLEVS-1)*NREP.LT.T) GO TO 50
      CALL ERROR(19)
      IF(BATCH)STOP
      RETURN
  50  NVAR1=NV
      NV=NV+NLEVS
      L1=INDIND+NVAR1*T
      IREP=1
      ICT=1
C AUTOMATIC GENERATION OF FACTOR LEVELS
      DO 80 IT=1,T
      DO 70 I=1,NREP
      L=L1
      DO 60 J=1,NLEVS
      DATA(L)=0.0
      IF(J.EQ.IREP)DATA(L)=1.0
      XBAR(NVAR1+J)=XBAR(NVAR1+J)+DATA(L)
  60  L=L+T
      L1=L1+1
      ICT=ICT+1
      IF(ICT.GT.T) GO TO 170
   70 CONTINUE
      IREP=IREP+1
  80  IF(IREP.GT.NLEVS) IREP=1
C READ IN FACTOR LEVELS
  90  NLEVS=NO(1)
      IF(NLEVS.GT.1.AND.NLEVS.LE.MAXIND-NV.AND.NLEVS.
     2LE.T) GO TO 100
      CALL ERROR(17)
      IF(BATCH)STOP
      RETURN
 100  IF(ECHO.AND.LFORM) WRITE(OUTCH,501)T
 501  FORMAT(1X,6HINPUT ,I2,24H FORMATTED FACTOR LEVELS )
      IF(ECHO.AND..NOT.LFORM)WRITE(OUTCH,502)T
 502  FORMAT(1X,6HINPUT ,I2,26H UNFORMATTED FACTOR LEVELS)
      IF(LFORM) GO TO 110
      NOBS=0
C  READ FACTOR LEVELS IN A-FORMAT
      DO 105 J3=1,T
C***PROMPT FOR NEW COMMAND - NOT NECESSARY ON SOME SYSTEMS
      IF(.NOT.BATCH.AND.INCH.EQ.IN5) WRITE(OUT6,507)
C***  NON-STANDARD FORMAT FOR PROMPT
 507  FORMAT(1X,1H*,'$')
      READ(INCH,503)(NOTE(I),I=2,NCHL)
      NOTE(1)=BLANK
 503  FORMAT(80A1)
      IF(BATCH) WRITE(OUTCH,504)(NOTE(I),I=1,NCHL)
 504  FORMAT(1X,80A1)
      IN=NOBS
      CALL INTERP(NOTE,NCHL,NO,IN,T,NOBS,IFL)
      IF(IFL.NE.0) GO TO 107
      IF(NOBS.EQ.T) GO TO 120
 105  CONTINUE
 107  CALL ERROR(5)
      IF(BATCH) STOP
      RETURN
 110  READ(INCH,FRMT)(NO(I),I=1,T)
      IF(ECHO.AND..NOT.BATCH) WRITE(OUT6,505)
  505 FORMAT(1X,50HFORMAT REMAINS IN FORCE. CANCEL OR CHANGE IF REQD.)
C
C  CHECK FACTOR LEVELS ARE VALID
C
 120  DO 140 I=1,T
      IF(NO(I).LE.0.OR.NO(I).GT.NLEVS) GO TO 130
      GO TO 140
 130  CALL ERROR(20)
      IF(BATCH)STOP
      RETURN
 140  CONTINUE
      NVAR1=NV
      NV=NV+NLEVS
      L1=INDIND+NVAR1*T
      DO 150 I=1,T
      L=L1+I-1
      DO 160 J=1,NLEVS
      DATA(L)=0.0
      IF(NO(I).EQ.J) DATA(L)=1.0
      XBAR(NVAR1+J)=XBAR(NVAR1+J)+DBLE(DATA(L))
 160  L=L+T
 150  CONTINUE
 170  NVAR1=NVAR1+1
      DO 180 J=NVAR1,NV
 180  XBAR(J)=XBAR(J)/DBLE(FLOAT(T))
      IF(ECHO)WRITE(OUTCH,506) NVAR1,NV
 506  FORMAT(1X,32HFACTOR LEVELS STORED IN VARIATES,I3,3H TO,I3)
      RETURN
      END
C  ############################################################
      SUBROUTINE XVAR
C  XVARIATE DECLARED WITH POSSIBLE OVERWRITE
      DOUBLE PRECISION RDATA,LOGLIK,Q,XBAR
      COMMON RDATA(5000),LOGLIK,XBAR(50),PROP(20),Q(400)
      REAL DATA(10000),NOTE,LINKNM
      EQUIVALENCE(RDATA(1),DATA(1))
      INTEGER T,OUT6,OUTCH,PRINT,ODD
      COMMON/SPEC/T,K,KMAX,LINK,NV,LINKNM(8),ODD,TITLE(80)
      COMMON/CONS/ZERO,ONE,HALF,TWO,RT2PI,PINV,ETA,TOL
      LOGICAL BATCH,ECHO,LREAD,LFORM,LTIT
      COMMON/CONTRL/INDSSP,INDIND,MAXDAT,MAXCYC,PRINT,LTIT
      COMMON/PERIPH/IN5,OUT6,INCH,OUTCH,BATCH
      COMMON/LOGIC/MAXIND,MAXCAT,LREAD,ECHO,IFL,LFORM
      COMMON/LEFT/FRMT(20)
      COMMON/NO/NO(20),NOTE(80),NOBS
      COMMON/ENVIR/NCHW,NCHL,NCHLM3,NCHLO
C
      IF(LREAD) GO TO 5
      CALL ERROR(6)
      IF(BATCH)STOP
      RETURN
  5   CALL INTERP(NOTE,NCHLM3,NO,0,1,NOBS,IFL)
      IF(IFL.EQ.0) GO TO 10
      CALL ERROR(5)
      IF(BATCH)STOP
      RETURN
  10  N2=NV+1
      IF(NOBS.EQ.1) N2=NO(1)
      IF(N2.LE.MAXIND.AND.N2.GT.0) GO TO 20
      CALL ERROR(19)
      IF(BATCH)STOP
      RETURN
  20  M=1
      N3=NV+1
      IF(N2.LE.N3) GO TO 30
      CALL ERROR(19)
      IF(BATCH)STOP
      RETURN
  30  IF(N2.LT.N3) M=0
      L1=INDIND+(N2-1)*T
      L2=L1+T-1
      IF(LFORM) GO TO 40
C** PROMPT FOR NEW DATA NOT NECESSARY ON SOME SYSTEMS
      IF(.NOT.BATCH.AND.INCH.EQ.IN5) WRITE(OUT6,503)
C*** NON-STARDARD FORMAT FOR PROMPT
 503  FORMAT(1X,1H*,'$')
C*** FREE FORMAT DATA INPUT
      READ(INCH,*)(DATA(J),J=L1,L2)
      IF(BATCH)WRITE(OUTCH,504)(DATA(J),J=L1,L2)
  504 FORMAT(1X,17F7.2)
      GO TO 50
C    FORMATTED INPUT
  40  READ(INCH,FRMT)(DATA(J),J=L1,L2)
      IF(BATCH)WRITE(OUTCH,504)(DATA(J),J=L1,L2)
      IF(ECHO.AND..NOT.BATCH)WRITE(OUTCH,501)
  501 FORMAT(1X,50HFORMAT REMAINS IN FORCE. CANCEL OR CHANGE IF REQD.)
  50  NV=NV+M
      XBAR(N2)=ZERO
      DO 60 I=1,T
  60  XBAR(N2)=XBAR(N2)+DBLE(DATA(L1+I-1))
      XBAR(N2)=XBAR(N2)/DBLE(FLOAT(T))
      IF(ECHO) WRITE(OUTCH,502) N2
 502  FORMAT(1X,19HVARIATE STORED IN X,I2)
      RETURN
      END
C  #################################################################
      SUBROUTINE COMBIN(N1,N2)
C  MERGES COLUMN N1 AND N2
      DOUBLE PRECISION RDATA,LOGLIK,Q,XBAR
      COMMON RDATA(5000),LOGLIK,XBAR(50),PROP(20),Q(400)
      LOGICAL LSCL,BATCH,LREAD,ECHO
      REAL DATA(1000),NOTE,LINKNM
      EQUIVALENCE(RDATA(1),DATA(1))
      INTEGER T,ODD,OUT6,OUTCH
      COMMON/SPEC/T,K,KMAX,LINK,NV,LINKNM(8),ODD,TITLE(80)
      DOUBLE PRECISION BETA
      COMMON/PARAM/BETA(100),INDEX1(50),INDEX2(50),NLOCP,NSCLP,LSCL
      COMMON/PERIPH/IN5,OUT6,INCH,OUTCH,BATCH
      COMMON/ENVIR/NCHW,NCHL,NCHLM3,NCHLO
      COMMON/LOGIC/MAXIND,MAXCAT,LREAD,ECHO,IFL,LFORM
      COMMON/NO/NO(20),NOTE(80),NOBS
C
      IF(LREAD)GOTO 4
      CALL ERROR(6)
      IF(BATCH)STOP
      RETURN
    4 IF(IFL.EQ.0) GO TO 8
      CALL INTERP(NOTE,NCHLM3,NO,0,2,NOBS,IFL)
      IF(IFL.EQ.0) GO TO 5
      CALL ERROR(5)
      IF(BATCH)STOP
      RETURN
   5  IF(NOBS.EQ.2) GO TO 6
      CALL ERROR(12)
      IF(BATCH) STOP
      RETURN
   6  N1=NO(1)
      N2=NO(2)
      IF(N1.LE.K.AND.N2.LE.K) GO TO 10
    7 CALL ERROR(15)
      IF(BATCH)STOP
      RETURN
   8  N3=N1
      N4=N2
      GO TO 20
  10  N3=MIN0(N1,N2)
      N4=MAX0(N1,N2)
      IF(N3.EQ.0)GOTO 7
      IF(N4-N3.EQ.1) GO TO 20
      CALL ERROR(21)
      IF(BATCH)STOP
      RETURN
  20  L=(N3-1)*T
      M1=L+T
      DO 30 J=1,T
  30  DATA(L+J)=DATA(L+J)+DATA(M1+J)
      K=K-1
      M2=K*T+T
C  MOVE UP OBSERVATION AND TOTALS
      M1=M1+1
      DO 40 J=M1,M2
  40  DATA(J)=DATA(J+T)
C ADJUST PROPORTIONS
      K1=K-1
      DO 50 J=N3,K1
      PROP(J)=PROP(J+1)
  50  BETA(J)=BETA(J+1)
      IF(ECHO) WRITE(OUTCH,501) N3,N4
 501  FORMAT(1X,10HCATEGORIES,2I4,7H MERGED)
      RETURN
      END
C  ######################################################
      SUBROUTINE RESET
C
      REAL NOTE,LINKNM
      DOUBLE PRECISION RDATA,LOGLIK,Q,XBAR
      COMMON RDATA(5000),LOGLIK,XBAR(50),PROP(20),Q(400)
      INTEGER T,ODD,PRINT
      LOGICAL LREAD,BATCH
      COMMON/SPEC/T,K,KMAX,LINK,NV,LINKNM(8),ODD,TITLE(80)
      DOUBLE PRECISION BETA
      COMMON/PARAM/BETA(100),INDEX1(50),INDEX2(50),NLOCP,NSCLP,LSCL
      COMMON/NO/NO(20),NOTE(80),NOBS
      COMMON/LOGIC/MAXIND,MAXCAT,LREAD,ECHO,IFL,LFORM
      COMMON/PERIPH/IN5,OUT6,INCH,OUTCH,BATCH
      DOUBLE PRECISION ZERO,PI,ONE,GAMA1,GAMA2
      DOUBLE PRECISION PPND
C
      DATA ZERO/0.0D0/,PI/3.1415927D0/,ONE/1.0D0/
      IF(LREAD)GOTO 5
      CALL ERROR(6)
      IF(BATCH)STOP
      RETURN
    5 DO 10 J=K,100
  10  BETA(J)=ZERO
  20  K1=K-1
      DO 30 J=1,K1
  30  BETA(J)=DBLE(PROP(J))
      GO TO (40,55,70,90),LINK
  40  DO 50 J=1,K1
  50  BETA(J)=DLOG(BETA(J)/(ONE-BETA(J)))
      RETURN
   55 DO 60 J=1,K1
  60  BETA(J)=PPND(BETA(J),IFL)
      RETURN
  70  DO 80 J=1,K1
      GAMA1=DSIN(BETA(J)*PI)
      GAMA2=DCOS(BETA(J)*PI)
  80  BETA(J)=-GAMA2/GAMA1
      RETURN
  90  DO 100 J=1,K1
      BETA(J)=-DLOG(ONE-BETA(J))
 100  BETA(J)=DLOG(BETA(J))
      RETURN
      END
C  ###############################################
      DOUBLE PRECISION FUNCTION PPND(P,IFAULT)
C
C  ALGORITHM AS 111  APPL. STATIST. (1977) P.120
C
      DOUBLE PRECISION ZERO,SPLIT,HALF,ONE
      DOUBLE PRECISION A0,A1,A2,A3,B1,B2,B3,B4,C0,C1,C2,C3,D1,D2
      DOUBLE PRECISION P,Q,R
      DATA ZERO/0.0D0/,HALF/0.5D0/,ONE/1.0D0/
      DATA SPLIT/0.42D0/
      DATA A0/2.50662823884D0/,A1/-18.61500062529D0/
      DATA A2/41.39119773534D0/,A3/-25.44106049637D0/
      DATA B1/-8.47351093090D0/,B2/23.08336743743D0/
      DATA B3/-21.06224101826D0/,B4/3.13082909833D0/
      DATA C0/-2.78718931138D0/,C1/-2.29796479134D0/
      DATA C2/4.85014127135D0/,C3/2.32121276858D0/
      DATA D1/3.54388924762D0/,D2/1.63706781897D0/
C
      IFAULT=0
      Q=P-HALF
      IF(DABS(Q).GT.SPLIT)GOTO 1
      R=Q*Q
      PPND=Q*(((A3*R+A2)*R+A1)*R+A0)/((((B4*R+B3)*R+B2
     1 )*R+B1)*R+ONE)
      RETURN
    1 R=P
      IF(Q.GT.ZERO)R=ONE-P
      IF(R.LE.ZERO)GOTO 2
      R=DSQRT(-DLOG(R))
      PPND=(((C3*R+C2)*R+C1)*R+C0)/((D2*R+D1)*R+ONE)
      IF(Q.LT.ZERO)PPND=-PPND
      RETURN
    2 IFAULT=1
      PPND=ZERO
      RETURN
      END
      SUBROUTINE INTERP(NOTE,NCHL,IARG,IN,INTEXP,INTOBS,IFAIL)
      REAL NOTE,NUM
      DIMENSION NOTE(NCHL),IARG(INTEXP),NUM(10)
      DATA NUM(1)/1H0/,NUM(2)/1H1/,NUM(3)/1H2/,NUM(4)/1H3/,
     * NUM(5)/1H4/,NUM(6)/1H5/,NUM(7)/1H6/,NUM(8)/1H7/,
     * NUM(9)/1H8/,NUM(10)/1H9/
      DATA BLANK/1H /,PLUS/1H+/,EX/1HX/
      IFAIL=0
      INTOBS=IN
      I=0
    1 I=I+1
      IF(I.GT.NCHL)RETURN
      CHAR=NOTE(I)
      IF(CHAR.NE.BLANK)GOTO 1
C BLANK CHARACTER IS SEPARATER - FIND NEXT SIGNIFICANT CHAR
    2 I=I+1
      IF(I.GT.NCHL)RETURN
      CHAR=NOTE(I)
      IF(CHAR.EQ.BLANK.OR.CHAR.EQ.EX.OR.CHAR.EQ.PLUS)GOTO 2
C CHECK CHARACTER IS AN INTEGER & CONSTRUCT NUMBER
      INT=0
    3 DO 4 J=1,10
      JJ=J
      IF(CHAR.EQ.NUM(J))GOTO 5
    4 CONTINUE
C NON INTEGER CHARACTER
      IFAIL=1
      RETURN
    5 INT=INT*10+JJ-1
      I=I+1
      IF(I.GT.NCHL)GOTO 6
      CHAR=NOTE(I)
      IF(CHAR.EQ.PLUS)GOTO 6
      IF(CHAR.NE.BLANK)GOTO 3
C END OF INTEGER FIELD
    6 INTOBS=INTOBS+1
      IARG(INTOBS)=INT
      IF(INTOBS.GE.INTEXP)RETURN
      GOTO 2
      END
       SUBROUTINE HELP
      INTEGER OUT6,OUTCH
      COMMON/PERIPH/IN5,OUT6,INCH,OUTCH,BATCH
      WRITE(OUTCH,500)
  500 FORMAT(1X,26HLIST OF SYSTEM DIRECTIVES  /
     * 1X,40HBATCH    COMBINE  DATA     ECHO     END /
     * 1X,40HFACTOR   FIT      FORMAT   FREE     HELP/
     * 1X,40HI/O      LINK     LOCATION LOOK     NOTE/
     * 1X,40HRESET    SCALE    STOP     TITLE    XVAR/
     * 1X,42HJACK     CYCLE    DUMP     RETRIEVE REWIND)
       RETURN
       END
C###############################################################
      SUBROUTINE ERROR(NERROR)
      INTEGER OUT6,OUTCH
      COMMON/PERIPH/IN5,OUT6,INCH,OUTCH,BATCH
      GOTO(10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,
     * 170,180,190,200,210,220),NERROR
   10 WRITE(OUT6,510)
  510 FORMAT(1X,30HERROR - NEGATIVE FITTED VALUE  /
     *1X,43HTYPE RESET AND BEGIN WITH FEWER COVARIATES  )
      RETURN
   20 WRITE(OUT6,520)
  520 FORMAT(1X,37HERROR - INSUFFICIENT SPACE AVAILABLE  )
      RETURN
   30 WRITE(OUT6,530)
  530 FORMAT(1X,34HERROR - INDEFINITE WORKING MATRIX  )
      RETURN
   40 WRITE(OUT6,540)
  540 FORMAT(1X,29HERROR - NEGATIVE OBSERVATION  )
      RETURN
   50 WRITE(OUT6,550)
  550 FORMAT(1X,42HERROR - ILLEGAL CHARACTER IN INTEGER FIELD )
      RETURN
   60 WRITE(OUT6,560)
  560 FORMAT(1X,15HERROR - NO DATA )
      RETURN
   70 WRITE(OUT6,570)
  570 FORMAT(1X,30HERROR - ONLY ONE CATEGORY LEFT  )
      RETURN
   80 WRITE(OUT6,580)
  580 FORMAT(1X,33HERROR - UNRECOGNISED SYSTEM WORD  )
      RETURN
   90 WRITE(OUT6,590)
  590 FORMAT(1X,31HERROR - NO COVARIATES DECLARED  )
      RETURN
  100 WRITE(OUT6,600)
  600 FORMAT(1X,28HERROR - INVALID LINK NUMBER  )
      RETURN
  110 WRITE(OUT6,610)
  610 FORMAT(1X,41HERROR - ILLEGAL ATTEMPT TO OVERWRITE DATA /
     *1X,14HTYPE END FIRST )
      RETURN
  120 WRITE(OUT6,620)
  620 FORMAT(1X,37HERROR - INSUFFICIENT NO. OF ARGUMENTS  )
      RETURN
  130 WRITE(OUT6,630)
  630 FORMAT(1X,49HERROR - INVALID JACKKNIFE ARGUMENT. INCREASE ARG.)
      RETURN
  140 WRITE(OUT6,640)
  640 FORMAT(1X,23HERROR - SPACE EXCEEDED  )
      RETURN
  150 WRITE(OUT6,650)
  650 FORMAT(1X,33HERROR - ARGUMENT(S) OUT OF RANGE  )
      RETURN
  160 WRITE(OUT6,660)
  660 FORMAT(1X,44HERROR - NO ARGUMENT GIVEN WHEN ONE EXPECTED  )
      RETURN
  170 WRITE(OUT6,670)
  670 FORMAT(1X,37HERROR - INVALID NO. OF FACTOR LEVELS  )
      RETURN
  180 WRITE(OUT6,680)
  680 FORMAT(1X,25HERROR - ILLEGAL COVARIATE  )
      RETURN
  190 WRITE(OUT6,690)
  690 FORMAT(1X,24HERROR - INVALID ARGUMENT  )
      RETURN
  200 WRITE(OUT6,700)
  700 FORMAT(1X,29HERROR - INVALID FACTOR LEVEL  )
      RETURN
  210 WRITE(OUT6,710)
  710 FORMAT(1X,50HERROR - ATTEMPT TO COMBINE NON-ADJACENT CATEGORIES)
      RETURN
  220 WRITE(OUT6,720)
  720 FORMAT(1X,28HERROR - TOO MANY COVARIATES  )
      RETURN
      END
C##############################################################
       SUBROUTINE LOOK
       DOUBLE PRECISION RDATA,LOGLIK,Q,XBAR
       COMMON RDATA(5000),LOGLIK,XBAR(50),PROP(20),Q(400)
       REAL DATA(10000),LINKNM,NOTE
       EQUIVALENCE (RDATA(1),DATA(1))
       INTEGER T,ODD,PRINT
       COMMON/SPEC/T,K,KMAX,LINK,NV,LINKNM(8),ODD,TITLE(80)
       LOGICAL LTIT,BATCH,LREAD,ECHO,LFORM
       COMMON/CONTRL/INDSSP,INDIND,MAXDAT,MAXCYC,PRINT,LTIT
       INTEGER OUT6,OUTCH
       COMMON/PERIPH/IN5,OUT6,INCH,OUTCH,BATCH
       COMMON/LOGIC/MAXIND,MAXCAT,LREAD,ECHO,IFL,LFORM
       COMMON/NO/NO(20),NOTE(80),NOBS
       COMMON/ENVIR/NCHW,NCHL,NCHLM3,NCHLO
      DIMENSION INDEX(20)
C
      DATA TOTAL/4HTOTL/,CODE/4H COL/
      NV1=NV+1
      IF(LREAD)GOTO 3
      CALL ERROR(6)
      IF(BATCH)STOP
      RETURN
    3 CALL INTERP(NOTE,NCHLM3,NO,0,NV1,NOBS,IFL)
      IF(IFL.EQ.0)GOTO 5
      CALL ERROR(5)
      IF(BATCH)STOP
      RETURN
    5 IF(NOBS.GT.0)GOTO 60
C PRINT DATA ARRAY
   10 WRITE(OUTCH,500)
  500 FORMAT(8X,11HDATA MATRIX )
      MAXL=NCHLO/7
      KP1=K+1
      IF(KP1.GT.MAXL)GOTO 30
C      WRITE(OUTCH,501)((CODE,J),J=1,K),TOTAL
  501 FORMAT(1X,4HUNIT,18(1X,A4,I2))
      KTT=K*T+T
      DO 20 I=1,T
   20 WRITE(OUTCH,502)I,(DATA(J),J=I,KTT,T)
  502 FORMAT(1X,I4,18F7.1)
      IF(NOBS.GT.1)GOTO 70
      RETURN
   30 MAXL=MAXL-1
C      WRITE(OUTCH,501)((CODE,J),J=1,MAXL)
      KTT=MAXL*T
      DO 40 I=1,T
   40 WRITE(OUTCH,502)I,(DATA(J),J=I,KTT,T)
      L1=MAXL*T
      MAXL=MAXL+1
      L2=K*T+T
C      WRITE(OUTCH,501)((CODE,J),J=MAXL,K),TOTAL
      DO 50 I=1,T
      L1=L1+1
   50 WRITE(OUTCH,502)I,(DATA(J),J=L1,L2,T)
      IF(NOBS.GT.1)GOTO 70
      RETURN
   60 IF(NO(1).EQ.0)GOTO 10
C PRINT X-VARIATES
   70 MAXL=NCHLO/11
      II=0
      DO 80 I=1,NOBS
      IF(NO(I).LE.0.OR.NO(I).GT.NV)GOTO 80
      II=II+1
      INDEX(II)=NO(I)
      IF(II.GT.MAXL)GOTO 90
   80 CONTINUE
   90 NOBS=II
      IF(NOBS.EQ.0)RETURN
      WRITE(OUTCH,503)(INDEX(I),I=1,NOBS)
  503 FORMAT(1X,10HCOVARIATES/1X,4HUNIT,11(I5,6X))
      DO 100 I=1,NOBS
  100 INDEX(I)=(INDEX(I)-1)*T+INDIND
      DO 110 I=1,T
      WRITE(OUTCH,504)I,(DATA(INDEX(J)),J=1,NOBS)
  504 FORMAT(1X,I3,1X,12G11.4)
      DO 110 J=1,NOBS
  110 INDEX(J)=INDEX(J)+1
      RETURN
      END
      SUBROUTINE CYCL(IFAIL)
C COMPUTES SSP MATRIX AND ADJUSTED DEPENDENT VARIATE  P.MC C.
      DOUBLE PRECISION RDATA,LOGLIK,WKSP,XBAR
      COMMON RDATA(5000),LOGLIK,XBAR(50),PROP(20),WKSP(400)
      REAL DATA(10000)
      EQUIVALENCE (DATA(1),RDATA(1))
      DOUBLE PRECISION ZERO,ONE,HALF,TWO,RT2PI,PINV,ETA,TOL
      COMMON/CONS/ZERO,ONE,HALF,TWO,RT2PI,PINV,ETA,TOL
      INTEGER T,ODD,PRINT
      REAL LINKNM
      COMMON/SPEC/T,K,KMAX,LINK,NV,LINKNM(8),ODD,TITLE(80)
      DOUBLE PRECISION BETA
      LOGICAL BATCH,LSCLP,LTIT
      COMMON/PARAM/BETA(100),INDEX1(50),INDEX2(50),NLOCP,NSCLP,LSCLP
      COMMON/CONTRL/INDSSP,INDIND,MAXDAT,MAXCYC,PRINT,LTIT
      COMMON/MODEL/NLOC,NSCL,NULL,ALIAS(100)
      INTEGER OUT6,OUTCH
      COMMON/PERIPH/IN5,OUT6,INCH,OUTCH,BATCH
      COMMON/ENVIR/NCHW,NCHL,NCHLM3,NCHLO
C LOCAL VARIABLES
      DOUBLE PRECISION ALNORM
      LOGICAL UPPER
      DOUBLE PRECISION Q(400)
      EQUIVALENCE (Q(1),WKSP(1))
      INTEGER R,S,INDX1(50),INDX2(50)
      DOUBLE PRECISION VLOC,DEV,VN,VSCL,LLINC,P,F,AINC,BINC,XTEMP1,
     1 XTEMP2,SUM1,SUM2,SUM3,SUM4,SUM5,V,VM1,FM1,W
      DATA UPPER/.FALSE./
C
      K1=K-1
      K2=K*2
      K3=K*3
      K4=K*4
      K5=K*5
      K6=K*6
      K7=K*7
      K8=K*8
      K9=K*9
      IF(NLOCP.EQ.0)GOTO 2
      DO 1 I=1,NLOCP
    1 INDX1(I)=INDIND+T*(INDEX1(I)-1)
    2 IF(.NOT.LSCLP)GOTO 4
      DO 3 I=1,NSCL
    3 INDX2(I)=INDIND+T*(INDEX2(I)-1)
    4 ILOC=(KMAX*T+T+1)/2+1
      ISCL=ILOC+T
      ITOT=T*K+1
      IMU=(ISCL+T-1)*2+1
      ILL=IMU+T*K
      LOGLIK=ZERO
      L1=INDSSP
      NP=NLOC+NSCL
      L2=L1+(NP*NP+NP)/2+NP-1
      DO 5 I=L1,L2
    5 RDATA(I)=ZERO
      Q(K3)=ONE
      Q(K4)=ZERO
      Q(K5)=ZERO
      Q(K6)=ZERO
      I=0
C LOOP OVER UNITS
    6 I=I+1
      LOCD=I
C COMPUTE LOCATION VALUE
      VLOC=ZERO
      IF(NLOCP.EQ.0)GOTO 8
      DO 7 J=1,NLOCP
      JJ=INDX1(J)
      JK=INDEX1(J)
      DEV=DBLE(DATA(JJ))-XBAR(JK)
    7 VLOC=VLOC+BETA(K1+J)*DEV
    8 RDATA(ILOC)=VLOC
C COMPUTE SCALE VALUE
      VSCL=ZERO
      IF(.NOT.LSCLP)GOTO 10
      DO 9 J=1,NSCL
      JJ=INDX2(J)
      JK=INDEX2(J)
      DEV=DBLE(DATA(JJ))-XBAR(JK)
    9 VSCL=VSCL-BETA(NLOC+J)*DEV
   10 W=DEXP(VSCL)
      RDATA(ISCL)=W
      VN=DBLE(DATA(ITOT))
      Q(1)=DBLE(DATA(LOCD))/VN
      DO 16 J=1,K1
      LOCD=LOCD+T
      Q(J+1)=Q(J)+DBLE(DATA(LOCD))/VN
      Q(K+J)=W*(BETA(J)-VLOC)
      GOTO (11,12,13,14),LINK
C LOGIT LINK
   11 Q(K2+J)=ONE/(ONE+DEXP(-Q(K+J)))
      Q(K3+J)=Q(K2+J)*(ONE-Q(K2+J))
      GOTO 16
C PROBIT LINK
   12 Q(K2+J)=ALNORM(Q(K+J),UPPER)
      Q(K3+J)=RT2PI*DEXP(-HALF*Q(K+J)**2)
      GOTO 16
C CAUCHIT LINK
   13 Q(K2+J)=PINV*DATAN(Q(K+J))+HALF
      Q(K3+J)=PINV/(ONE+Q(K+J)**2)
      GOTO 16
C COMPLEMENTARY LOG-LOG LINK
   14 Q(K2+J)=ONE-DEXP(-DEXP(Q(K+J)))
      Q(K3+J)=DEXP(Q(K+J)-DEXP(Q(K+J)))
   16 CONTINUE
C
      LLINC=ZERO
      IF(Q(1).GT.ZERO)LLINC=DLOG(Q(K2+1)/Q(1))*Q(1)
      DO 18 J=1,K1
      P=Q(J+1)-Q(J)
      F=Q(K2+J+1)-Q(K2+J)
      IF(P.LT.ZERO)GOTO 1001
      IF(F.LE.ZERO)GOTO 1000
      IF(P.EQ.ZERO)GOTO 17
      LLINC=LLINC+P*DLOG(F/P)
   17 Q(K4+J)=Q(K2+J)-Q(K2+J)**2/Q(K2+J+1)
      IFIT=IMU+T*J
   18 DATA(IFIT)=SNGL(VN*F)
      DATA(IMU)=SNGL(VN*Q(K2+1))
      DATA(ILL)=-2.*SNGL(VN*LLINC)
      SUM1=ZERO
      SUM2=ZERO
      SUM3=ZERO
      SUM4=ZERO
      SUM5=ZERO
      F=ZERO
      V=ONE
      DO 19 J=1,K1
      FM1=F
      F=Q(K2+J)/Q(K2+J+1)
      VM1=V
      V=Q(K4+J)
      Q(K5+J)=Q(K3+J)-F*Q(K3+J+1)
      Q(K6+J)=Q(K+J)*Q(K3+J)-F*Q(K+J+1)*Q(K3+J+1)
      Q(K7+J)=Q(K3+J)*(Q(K5+J)/V-FM1*Q(K5+J-1)/VM1)
      Q(K8+J)=Q(K3+J)*(Q(K6+J)/V-FM1*Q(K6+J-1)/VM1)
      Q(K9+J)=Q(K6+J)*(ONE+VSCL)+Q(J)-F*Q(J+1)
      SUM1=SUM1+Q(K5+J)**2/V
      SUM2=SUM2+Q(K5+J)*Q(K6+J)/V
      SUM3=SUM3+Q(K6+J)**2/V
      SUM4=SUM4+Q(K5+J)*Q(K9+J)/V
      SUM5=SUM5+Q(K6+J)*Q(K9+J)/V
   19 CONTINUE
C
      LOCA=INDSSP
      LOCB=INDSSP+(NP*NP+NP)/2
C AUGMENT SSP MATRIX AND ADJUSTED DEPENDENT VARIATE
      DO 32 R=1,NP
      IF(R.GE.K)GOTO 23
C SECOND DERIVATIVES W.R.T. THETA
      DO 22 S=1,R
      IF(S.LT.R-1)GOTO 22
      IF(R.NE.S)GOTO 20
C R=S
      AINC=ONE/Q(K4+S)
      IF(R.NE.1)AINC=AINC+(Q(K2+S-1)/Q(K2+S))**2/Q(K4+S-1)
      AINC=AINC*(W*Q(K3+S))**2
      GOTO 21
C S=R-1
   20 AINC=-W**2*Q(K3+S)*Q(K3+R)*Q(K2+S)/Q(K2+R)/Q(K4+S)
   21 RDATA(LOCA)=RDATA(LOCA)+VN*AINC
   22 LOCA=LOCA+1
C INCREMENT B
      BINC=Q(K9+R)/Q(K4+R)
      IF(R.NE.1)BINC=BINC-Q(K9+R-1)*Q(K2+R-1)/Q(K2+R)/Q(K4+R-1)
      RDATA(LOCB)=RDATA(LOCB)+VN*W*Q(K3+R)*BINC
      LOCB=LOCB+1
      GOTO 32
   23 IF(R.GT.NLOC)GOTO 27
      JJ=INDX1(R-K1)
      JK=INDEX1(R-K1)
      XTEMP1=XBAR(JK)-DBLE(DATA(JJ))
      DO 26 S=1,R
      IF(S.GE.K)GOTO 24
C MIXED DERIVATIVES - THETA AND LOCATION
      AINC=XTEMP1*W*W*Q(K7+S)
      GOTO 25
C DERIVATIVES W.R.T. LOCATION PARAMETERS
   24 JJ=INDX1(S-K1)
      JK=INDEX1(S-K1)
      XTEMP2=XBAR(JK)-DBLE(DATA(JJ))
      AINC=XTEMP1*XTEMP2*W*W*SUM1
   25 RDATA(LOCA)=RDATA(LOCA)+VN*AINC
   26 LOCA=LOCA+1
C INCREMENT B
      BINC=W*XTEMP1*SUM4
      RDATA(LOCB)=RDATA(LOCB)+VN*BINC
      LOCB=LOCB+1
      GOTO 32
C SCALE PARAMETERS
   27 JJ=INDX2(R-NLOC)
      JK=INDEX2(R-NLOC)
      XTEMP1=XBAR(JK)-DBLE(DATA(JJ))
      DO 31 S=1,R
      IF(S.GE.K)GOTO 28
C MIXED DERIVATIVES - THETA AND SCALE
      AINC=XTEMP1*W*Q(K8+S)
      GOTO 30
   28 IF(S.GT.NLOC)GOTO 29
C MIXED DERIVATIVES - LOCATION AND SCALE
      JJ=INDX1(S-K1)
      JK=INDEX1(S-K1)
      XTEMP2=XBAR(JK)-DBLE(DATA(JJ))
      AINC=XTEMP1*XTEMP2*W*SUM2
      GOTO 30
C DERIVATIVES W.R.T. SCALE PARAMETERS
   29 JJ=INDX2(S-NLOC)
      JK=INDEX2(S-NLOC)
      XTEMP2=XBAR(JK)-DBLE(DATA(JJ))
      AINC=XTEMP1*XTEMP2*SUM3
   30 RDATA(LOCA)=RDATA(LOCA)+VN*AINC
   31 LOCA=LOCA+1
C INCREMENT B
      BINC=XTEMP1*SUM5
      RDATA(LOCB)=RDATA(LOCB)+VN*BINC
      LOCB=LOCB+1
   32 CONTINUE
C INCREMENT INDICES
      IF(NLOCP.EQ.0)GOTO 34
      DO 33 J=1,NLOCP
   33 INDX1(J)=INDX1(J)+1
   34 IF(.NOT.LSCLP)GOTO 36
      DO 35 J=1,NSCL
   35 INDX2(J)=INDX2(J)+1
   36 ILOC=ILOC+1
      ISCL=ISCL+1
      ILL=ILL+1
      IMU=IMU+1
      ITOT=ITOT+1
      LOGLIK=LOGLIK-TWO*LLINC*VN
      IF(I.LT.T)GOTO 6
      IFAIL=0
      RETURN
C NEGATIVE FITTED VALUE
 1000 CALL ERROR(1)
      IF(BATCH)STOP
      IFAIL=1
      RETURN
C NEGATIVE OBSERVATION
 1001 CALL ERROR(4)
      IF(BATCH)STOP
      IFAIL=1
      RETURN
      END
      SUBROUTINE FIT
C FITS MODEL AND PRINTS RESULTS         P.MC C.  JUNE 1978
      DOUBLE PRECISION RDATA,LOGLIK,WKSP,XBAR
      COMMON RDATA(5000),LOGLIK,XBAR(50),PROP(20),WKSP(400)
      REAL DATA(10000)
      EQUIVALENCE (DATA(1),RDATA(1))
      DOUBLE PRECISION ZERO,ONE,HALF,TWO,RT2PI,PINV,ETA,TOL
      COMMON/CONS/ZERO,ONE,HALF,TWO,RT2PI,PINV,ETA,TOL
      INTEGER T,ODD,PRINT
      REAL LINKNM
      COMMON/SPEC/T,K,KMAX,LINK,NV,LINKNM(8),ODD,TITLE(80)
      DOUBLE PRECISION BETA
      LOGICAL BATCH,LSCLP,LTIT,ERR
      COMMON/ERRB/ERR
      COMMON/PARAM/BETA(100),INDEX1(50),INDEX2(50),NLOCP,NSCLP,LSCLP
      COMMON/CONTRL/INDSSP,INDIND,MAXDAT,MAXCYC,PRINT,LTIT
      COMMON/MODEL/NLOC,NSCL,NULL,ALIAS(100)
      INTEGER OUT6,OUTCH
      COMMON/PERIPH/IN5,OUT6,INCH,OUTCH,BATCH
      COMMON/ENVIR/NCHW,NCHL,NCHLM3,NCHLO
C LOCAL VARIABLES
      DOUBLE PRECISION X,VN,BIG,FAC,GSQ
      DIMENSION X(100),RX(800),SX(200)
      EQUIVALENCE (RX(1),WKSP(1)),(SX(1),RX(601))
      INTEGER CYCLE,SPACE
      LOGICAL CONV
      DATA BIG/1.D20/
C
      GSQ=BIG
      ERR=.FALSE.
      INDSSP=(INDIND+T*NV)/2+1
      IF(LTIT.AND.PRINT.GE.0)WRITE(OUTCH,100)(TITLE(I),I=1,NCHLM3)
  100 FORMAT(1X,80A1)
      K1=K-1
      CYCLE=0
      NLOC=NLOCP+K1
      NP=NLOC+NSCL
C CALCULATE TOTAL SPACE REQUIRED
      SPACE=NP*NP+3*NP+(INDSSP-1)*2
      IF(SPACE.LE.MAXDAT)GOTO 1
      CALL ERROR(2)
      IF(BATCH)STOP
      RETURN
    1 CALL CYCL(IFAIL)
C##
      IF(PRINT.NE.9)GOTO 503
      WRITE(5,499)(X(I),I=1,NP)
  499 FORMAT(1X,4HX = ,7F9.4)
      WRITE(5,500)(WKSP(J),J=1,40)
  500 FORMAT(1X,4G12.4)
      NP1=NP+1
      L2=INDSSP-1
      DO 501 I=1,NP1
      L1=L2+1
      L2=L2+I
  501 WRITE(5,502)(RDATA(J),J=L1,L2)
  502 FORMAT(1X,6D12.4)
  503 CONTINUE
C##
      ERR=IFAIL.NE.0
      IF(IFAIL.NE.0)RETURN
      IF(LOGLIK.LT.GSQ)GOTO 24
      CYCLE=CYCLE+1
      IF(CYCLE.GT.MAXCYC)GOTO 6
      FAC=DSQRT(GSQ/LOGLIK)*HALF
      DO 2 I=1,NP
    2 BETA(I)=FAC*BETA(I)+(ONE-FAC)*X(I)
      GOTO 1
   24 DO 26 I=1,NP
   26 X(I)=BETA(I)
      GSQ=LOGLIK
      CALL SOLVE(INDSSP,BETA,NP,NULL,ALIAS,1,IFAIL)
      ERR=IFAIL.NE.0
      IF(IFAIL.EQ.0)GOTO 3
C NON POSITIVE SEMI-DEFINITE MATRICES ARE THOUGHT TO BE
C IMPOSSIBLE UNDER THE MODEL
      CALL ERROR(3)
      RETURN
    3 IF(PRINT.LT.6)GOTO 4
      WRITE(OUTCH,101)LOGLIK,CYCLE
  101 FORMAT(1X,11HDEVIANCE = ,D14.6,10H  CYCLE = ,I2)
      IF(PRINT.EQ.6)GOTO 4
      IF(NCHLO.LE.80)WRITE(OUTCH,102)(X(I),I=1,NP)
      IF(NCHLO.GT.80)WRITE(OUTCH,103)(X(I),I=1,NP)
  102 FORMAT(1X,5D14.6)
  103 FORMAT(1X,9D14.6)
C CHECK FOR CONVERGENCE
    4 CYCLE=CYCLE+1
      CONV=.TRUE.
      IF(MAXCYC.EQ.0)GOTO 6
      DO 5 I=1,NP
      IF(DABS(BETA(I)).LT.ETA)GOTO 5
      IF(DABS(X(I)-BETA(I)).GT.TOL*DABS(BETA(I)))CONV=.FALSE.
    5 CONTINUE
      IF(CONV)GOTO 7
      IF(CYCLE.LE.MAXCYC)GOTO 1
    6 LOGLIK=GSQ
      DO 64 I=1,NP
   64 BETA(I)=X(I)
      WRITE(OUTCH,104)CYCLE
  104 FORMAT(1X,24HNO CONVERGENCE BY CYCLE   ,I2)
    7 IDF=T*K1-NP+NULL
      L1=2*LINK-1
      L2=L1+1
      IF(PRINT.LT.0)RETURN
      WRITE(OUTCH,105)
  105 FORMAT(1X,28H DEVIANCE    DF. CYCLE  LINK   )
      WRITE(OUTCH,106)LOGLIK,IDF,CYCLE,LINKNM(L1),LINKNM(L2)
  106 FORMAT(1X,D11.4,I5,I6,2X,2A4)
      IF(PRINT.LE.0)RETURN
C OBTAIN MATRIX INVERSE
      CALL SYMINV(INDSSP,NP)
      WRITE(OUTCH,107)
  107 FORMAT(1X,44H PARAMETER   ESTIMATE   ST.ERROR ALIAS  XBAR )
C COMPUTE STANDARD DEVIATIONS
      LDIAG=INDSSP-1
      DO 8 I=1,NP
      LDIAG=LDIAG+I
    8 X(I)=DSQRT(RDATA(LDIAG))
      DO 9 I=1,K1
    9 WRITE(OUTCH,108)I,BETA(I),X(I),ALIAS(I)
  108 FORMAT(4X,2HTH,I2,2X,D12.4,D11.3,3X,A1)
      IF(NLOCP.EQ.0)GOTO 11
      DO 10 I=K,NLOC
      J=I-K+1
      JJ=INDEX1(J)
   10 WRITE(OUTCH,110)JJ,BETA(I),X(I),ALIAS(I),XBAR(JJ)
  110 FORMAT(4X,2HLX,I2,2X,D12.4,D11.3,3X,A1,G9.3)
   11 IF(.NOT.LSCLP)GOTO 13
      NLOC1=NLOC+1
      DO 12 I=NLOC1,NP
      J=I-NLOC
      JJ=INDEX2(J)
   12 WRITE(OUTCH,111)JJ,BETA(I),X(I),ALIAS(I),XBAR(JJ)
  111 FORMAT(4X,2HSX,I2,2X,D12.4,D11.3,3X,A1,G9.3)
   13 IF(PRINT.LE.1)RETURN
C COMPUTE AND PRINT LOCATION, SCALE, DEV/UNIT AND RESIDUALS
      WRITE(OUTCH,112)
  112 FORMAT(1X,52HUNIT LOCATION SCALE DEV/UNIT  STANDARDISED RESIDUALS)
C SET UP INDICES FOR LOCATION AND SCALE
      L1=(T*(KMAX+1)+1)/2+1
      L2=L1+T
C INDICES FOR DEV/UNIT AND FITTED VALUES IN SINGLE PRECISION
      L3=T*(K+KMAX+5)+ODD+1
      L4=T*(KMAX+5)+ODD+1
      KT=K*T
      DO 15 I=1,T
      X1=SNGL(RDATA(L1))
      X2=1./SNGL(RDATA(L2))
      X3=DATA(L3)
      L5=I
      L=L4
      TOT=DATA(KT+I)
      DO 14 J=1,K
      O=DATA(L5)
      F=DATA(L)
      RX(L5)=(O-F)/SQRT(F*(1.-F/TOT))
      L5=L5+T
   14 L=L+T
      L1=L1+1
      L2=L2+1
      L3=L3+1
      L4=L4+1
C PRINT ONE LINE OF RESULTS
   15 WRITE(OUTCH,113)I,X1,X2,X3,(RX(J),J=I,KT,T)
  113 FORMAT(1X,I4,F8.3,F7.3,F9.3,12F6.2)
      IF(PRINT.LE.2)RETURN
C COMPUTE CORRELATION MATRIX
      L2=INDSSP-1
      NCL=(NCHLO)/6-1
      NBL=(NP-1)/NCL+1
      WRITE(OUTCH,116)
  116 FORMAT(1X,25HCORRELATIONS OF ESTIMATES )
      DO 18 IBL=1,NBL
      I1=(IBL-1)*NCL+1
      DO 17 I=I1,NP
      DO 155 J=I1,I
  155 SX(J)=0.0
      L1=L2+I1
      L2=L2+I
      L=L1-1
      ICOL=I
      IF(I-I1+1.GT.NCL)ICOL=NCL+I1-1
      IF(X(I).EQ.ZERO)GOTO 17
      DO 16 J=I1,ICOL
      L=L+1
      IF(X(J).EQ.ZERO)GOTO 16
      SX(J)=SNGL(RDATA(L)/X(I)/X(J))
   16 CONTINUE
   17 WRITE(OUTCH,114)I,(SX(J),J=I1,ICOL)
      L2=INDSSP+(NCL*IBL)*(NCL*IBL+1)/2-1
   18 WRITE(OUTCH,115)(J,J=I1,ICOL)
  115 FORMAT(3X,20(I4,2X))
  114 FORMAT(1X,I2,20F6.3)
      IF(PRINT.LE.3)RETURN
C PRINT OBSERVED, FITTED AND STANDARDISED RESIDUALS
      IF(K.EQ.2)GOTO 20
      WRITE(OUTCH,117)
  117 FORMAT(1X,26HOBSERVED AND FITTED VALUES  )
      L1=T*(KMAX+5)+ODD+1
      L2=L1+KT-1
      KT1=K*T+T
      DO 19 I=1,T
      WRITE(OUTCH,119)(DATA(J),J=I,KT1,T)
      L=L1+I-1
      WRITE(OUTCH,119)(DATA(J),J=L,L2,T)
   19 WRITE(OUTCH,120)
  119 FORMAT(1X,15F8.2)
  120 FORMAT(1X)
      GOTO 22
C BINARY DATA
   20 L1=T*(KMAX+5)+ODD+1
      L2=T*K+1
      WRITE(OUTCH,121)
  121 FORMAT(7X,30HN OBSERVED    FITTED  RESIDUAL  )
      DO 21 I=1,T
      TOT=DATA(L2)
      O=DATA(I)
      F=DATA(L1)
      P=F/TOT
      RES=(O-F)/SQRT(TOT*P*(1.-P))
      RX(I)=RES
      WRITE(OUTCH,122)TOT,O,F,RES
  122 FORMAT(1X,2F8.1,2F10.2)
      L1=L1+1
   21 L2=L2+1
C OPTION PRINT=5 RESERVED FOR RANKIT AND OTHER PLOTS IS
C NOT YET AVAILABLE
   22 RETURN
      END
      SUBROUTINE SOLVE(LOC,BETA,N,NULL,ALIAS,IA,IFL)
C  ALGORITHM A.S. 6   APPL.STATIST. (1968) VOL.17 NO.2
C  ADAPTED P.MC C. APRIL 1978
      DOUBLE PRECISION RDATA,LOGLIK,Q,XBAR
      COMMON RDATA(5000),LOGLIK,XBAR(50),P(20),Q(400)
      DOUBLE PRECISION ZERO,ONE,HALF,TWO,RT2PI,PINV,ETA,TOL
      COMMON/CONS/ZERO,ONE,HALF,TWO,RT2PI,PINV,ETA,TOL
      DOUBLE PRECISION BETA(N),W,W1
      REAL ALIAS(100)
      DATA BLANK,STAR/1H ,1H*/
C
      IFL=1
      IF(N.LE.0)GOTO 100
      IFL=2
      NULL=0
      N1=N+1
      J=LOC
      K=LOC-1
      DO 10 ICOL=1,N1
      L=LOC-1
      DO 11 IROW=1,ICOL
      IF(IROW.EQ.N1)GOTO 40
      K=K+1
      W=RDATA(K)
      M=J
      DO 12 I=1,IROW
      L=L+1
      IF(I.EQ.IROW)GOTO 13
      W=W-RDATA(L)*RDATA(M)
      M=M+1
   12 CONTINUE
   13 IF(IROW.EQ.ICOL)GOTO 14
      IF(RDATA(L).EQ.ZERO)GOTO 21
      RDATA(K)=W/RDATA(L)
      GOTO 11
   21 RDATA(K)=ZERO
   11 CONTINUE
   14 IF(DABS(W).LE.DABS(ETA*RDATA(K)))GOTO 20
      IF(W.LT.ZERO)GOTO 100
      RDATA(K)=DSQRT(W)
      ALIAS(IA+ICOL-1)=BLANK
      GOTO 15
   20 RDATA(K)=ZERO
      NULL=NULL+1
      ALIAS(IA+ICOL-1)=STAR
   15 J=J+ICOL
   10 CONTINUE
  100 RETURN
   40 IFL=0
C BACK SUBSTITUTION
      J=LOC+(N*N+N)/2-1
      K=J+N
      ICOL=N
      W=RDATA(J)
      BETA(N)=ZERO
      IF(W.GT.ZERO)BETA(N)=RDATA(K)/W
      IF(N.EQ.1)GOTO 32
   29 L=J-1
      K=K-1
      J=J-ICOL
      W=RDATA(J)
      W1=RDATA(K)
      BETA(ICOL-1)=ZERO
      IF(W.EQ.ZERO)GOTO 31
      DO 30 JCOL=ICOL,N
      W1=W1-RDATA(L)*BETA(JCOL)
   30 L=L+JCOL
      BETA(ICOL-1)=W1/W
   31 ICOL=ICOL-1
      IF(ICOL.GT.1)GOTO 29
   32 RETURN
      END
      DOUBLE PRECISION FUNCTION ALNORM(X,UPPER)
C ALGORITHM A.S.66 APPLIED STATISTICS (1973) VOL22 NO.3
      DOUBLE PRECISION ZERO,ONE,HALF,TWO,RT2PI,PINV,ETA,TOL
      COMMON/CONS/ZERO,ONE,HALF,TWO,RT2PI,PINV,ETA,TOL
      DOUBLE PRECISION CON,Z,Y,X
      DOUBLE PRECISION P,Q,R,A1,A2,A3,B1,B2,C1,C2,C3,C4,C5,C6
      DOUBLE PRECISION D1,D2,D3,D4,D5
      LOGICAL UPPER,UP
C*** MACHINE DEPENDENT CONSTANTS
      DOUBLE PRECISION LTONE,UTZERO
      DATA LTONE/7.0D0/,UTZERO/18.66D0/
      DATA CON/1.28D0/
      DATA P/0.398942280444D0/,Q/0.39990348504D0/,R/0.398942280385D0/
      DATA A1/5.75885480458D0/,A2/2.62433121679D0/,A3/5.92885724438D0/
      DATA B1/-29.8213557807D0/,B2/48.6959930692D0/
      DATA C1/-3.8052D-8/,C2/3.98064794D-4/,C3/-0.151679116635D0/
      DATA C4/4.8385912808D0/,C5/0.742380924027D0/,C6/3.99019417011D0/
      DATA D1/1.00000615302D0/,D2/1.98615381364D0/,D3/5.29330324926D0/
      DATA D4/-15.1508972451D0/,D5/30.789933034D0/
C
      UP=UPPER
      Z=X
      IF(Z.GE.ZERO)GOTO 10
      UP=.NOT.UP
      Z=-Z
   10 IF(Z.LE.LTONE.OR.UP.AND.Z.LE.UTZERO)GOTO 20
      ALNORM=ZERO
      GOTO 40
   20 Y=HALF*Z*Z
      IF(Z.GT.CON) GOTO 30
C
      ALNORM=HALF-Z*(P-Q*Y/(Y+A1+B1/(Y+A2+B2/(Y+A3))))
      GOTO 40
   30 ALNORM=R*DEXP(-Y)/(Z+C1+D1/(Z+C2+D2/(Z+C3+D3/(Z+C4+D4/(Z+C5+D5/
     2   (Z+C6))))))
   40 IF(.NOT.UP)ALNORM=ONE-ALNORM
      RETURN
      END
      SUBROUTINE JACK(FRAC)
C JACKKNIFES MAXIMUM LIKELIHOOD ESTIMATES IN DECLARED MODEL
      DOUBLE PRECISION RDATA,LOGLIK,Q,XBAR
      COMMON RDATA(5000),LOGLIK,XBAR(50),PROP(20),Q(400)
      REAL DATA(10000)
      EQUIVALENCE (DATA(1),RDATA(1))
      DOUBLE PRECISION ZERO,ONE,HALF,TWO,RT2PI,PINV,ETA,TOL
      COMMON/CONS/ZERO,ONE,HALF,TWO,RT2PI,PINV,ETA,TOL
      INTEGER T,ODD,PRINT
      REAL LINKNM
      COMMON/SPEC/T,K,KMAX,LINK,NV,LINKNM(8),ODD,TITLE(80)
      DOUBLE PRECISION BETA
      LOGICAL BATCH,LSCLP,LTIT,ERR
      COMMON/ERRB/ERR
      COMMON/PARAM/BETA(100),INDEX1(50),INDEX2(50),NLOCP,NSCLP,LSCLP
      COMMON/CONTRL/INDSSP,INDIND,MAXDAT,MAXCYC,PRINT,LTIT
      COMMON/MODEL/NLOC,NSCL,NULL,ALIAS(100)
      INTEGER OUT6,OUTCH
      COMMON/PERIPH/IN5,OUT6,INCH,OUTCH,BATCH
      COMMON/ENVIR/NCHW,NCHL,NCHLM3,NCHLO
      DOUBLE PRECISION F,FAC,SUM(50),SUMSQ(50),WT,TOTWT
      EQUIVALENCE (SUM(1),Q(301)),(SUMSQ(1),Q(351))
C
      K1=K-1
      MAXCYC=15
      NLOC=NLOCP+K-1
      TOL=TOL*0.01D0
      CALL FIT
      IF(ERR)RETURN
      NP=NLOCP+NSCLP+K-1
      LOCN=INDSSP+NP+NP*(NP+1)/2
      LOCNP=LOCN+NP-1
      DO 10 I=LOCN,LOCNP
      RDATA(I)=BETA(I-LOCN+1)
      RDATA(NP+I)=ZERO
   10 RDATA(2*NP+I)=ZERO
      ITOT=K*T+1
      TOTWT=ZERO
      DO 40 I=1,T
      LOCD=I
      RED=1.0
      IF(DATA(ITOT).GT.1./FRAC)RED=FRAC*DATA(ITOT)
      FAC=DBLE((DATA(ITOT)-1.)/DATA(ITOT))
      DATA(ITOT)=DATA(ITOT)-RED
      DO 12 J=1,NP
      SUM(J)=ZERO
      SUMSQ(J)=ZERO
   12 CONTINUE
      DO 30 J=1,K
      IF(DATA(LOCD).LE.0.0)GOTO 30
      IF(DATA(LOCD).GT.RED)GOTO 15
      CALL ERROR(13)
      IF(BATCH)STOP
      DATA(ITOT)=DATA(ITOT)+RED
      RETURN
   15 F=DBLE(DATA(LOCD))
      DATA(LOCD)=DATA(LOCD)-RED
      CALL FIT
      IF(.NOT.ERR)GOTO 17
      DATA(LOCD)=DATA(LOCD)+RED
      DATA(ITOT)=DATA(ITOT)+RED
      RETURN
   17 DO 20 JJ=1,NP
      SUM(JJ)=SUM(JJ)+F*BETA(JJ)
   20 SUMSQ(JJ)=SUMSQ(JJ)+F*BETA(JJ)**2
      DATA(LOCD)=DATA(LOCD)+RED
   30 LOCD=LOCD+T
      DATA(ITOT)=DATA(ITOT)+RED
      WT=DBLE(DATA(ITOT))
      DO 35 JJ=LOCN,LOCNP
      IJ=JJ-LOCN+1
      RDATA(NP+JJ)=RDATA(NP+JJ)+SUM(IJ)*FAC/DBLE(RED)
   35 RDATA(2*NP+JJ)=RDATA(2*NP+JJ)+FAC*(SUMSQ(IJ)-SUM(IJ)**2/WT)
     * /DBLE(RED)**2
      ITOT=ITOT+1
   40 TOTWT=TOTWT+(WT-ONE)/DBLE(RED)
      DO 50 I=1,NP
      I1=LOCN+I-1
      I2=I1+NP
      I3=I2+NP
      Q(I)=RDATA(I1)*(TOTWT+ONE)-RDATA(I2)
   50 Q(NP+I)=DSQRT(RDATA(I3))
C
      WRITE(OUTCH,500)
  500 FORMAT(1X,41HJACKKNIFED ESTIMATES AND ESTIMATED S.E.S /
     *  1X,32H PARAMETER   ESTIMATE   ST.ERROR  )
      DO 100 I=1,K1
  100 WRITE(OUTCH,501)I,Q(I),Q(NP+I)
  501 FORMAT(4X,2HTH,I2,2X,D12.4,D11.3)
      IF(NLOCP.EQ.0)GOTO 120
      DO 110 I=K,NLOC
      J=I-K+1
      JJ=INDEX1(J)
  110 WRITE(OUTCH,502)JJ,Q(I),Q(NP+I)
  502 FORMAT(4X,2HLX,I2,2X,D12.4,D11.3)
  120 IF(.NOT.LSCLP)GOTO 140
      NLOC1=NLOC+1
      DO 130 I=NLOC1,NP
      J=I-NLOC
      JJ=INDEX2(J)
  130 WRITE(OUTCH,503)JJ,Q(I),Q(NP+I)
  503 FORMAT(4X,2HSX,I2,2X,D12.4,D11.3)
C
  140 DO 150 I=1,NP
  150 BETA(I)=Q(I)
      TOL=TOL*100.0D0
      RETURN
      END
      SUBROUTINE SYMINV(INDSSP,N)
C  ALGORITHM A.S. 7   J.R.S.S.(C), (1968) VOL.17  NO.2.
C  ADAPTED P.MC. C.  13/4/78
      DOUBLE PRECISION RDATA,LOGLIK,W,XBAR
      COMMON RDATA(5000),LOGLIK,XBAR(50),P(20),W(400)
      DOUBLE PRECISION X,ZERO,ONE
      DATA ZERO/0.0D0/,ONE/1.0D0/
C
      NROW=N
      NN=(NROW*NROW+NROW)/2+INDSSP-1
      IROW=NROW
      NDIAG=NN
   16 IF(RDATA(NDIAG).EQ.ZERO)GOTO 11
      L=NDIAG
      DO 10 I=IROW,NROW
      W(I)=RDATA(L)
      L=L+I
   10 CONTINUE
      ICOL=NROW
      JCOL=NN
      MDIAG=NN
   15 L=JCOL
      X=ZERO
      IF(ICOL.EQ.IROW)X=ONE/W(IROW)
      K=NROW
   13 IF(K.EQ.IROW)GOTO 12
      X=X-W(K)*RDATA(L)
      K=K-1
      L=L-1
      IF(L.GT.MDIAG)L=L-K+1
      GOTO 13
   12 RDATA(L)=X/W(IROW)
      IF(ICOL.EQ.IROW)GOTO 14
      MDIAG=MDIAG-ICOL
      ICOL=ICOL-1
      JCOL=JCOL-1
      GOTO 15
   11 L=NDIAG
      DO 17 J=IROW,NROW
      RDATA(L)=ZERO
      L=L+J
   17 CONTINUE
   14 NDIAG=NDIAG-IROW
      IROW=IROW-1
      IF(IROW.NE.0)GOTO 16
      RETURN
      END


