C
C
C
C      PROGRAM MAIN (INPUT,OUTPUT,LPOUT,TAPE5=INPUT,TAPE6=OUTPUT,
C
C     1              TAPE10=LPOUT)
C
C*****************************************************************
C*      THIS VERSION OF INVAR (COMPLETED AUGUST 1982) INCLUDES   *
C*      ALL OF THE ORIGINAL 18 SUBROUTINES (WITH MINOR ADJUST-   *
C*      MENTS) BY RALL AND FUNDERLIC (MAY 1980) EXCEPT FOR PLOT. *
C*      THE SUBROUTINE PLOT HAS BEEN REPLACED BY THE SUBROU-     *
C*      TINES SETUP1, SETUP2, SPLT, SPLTH, PLTMAIN, ORPLOT,      *
C*      VPLT, VPLTER, VPMAIN AND VPMAINH, WHICH WERE TAKEN       *
C*      FROM STATLIB (AGAIN WITH MINOR ADJUSTMENTS) BY P.V.      *
C*      TRYON AND J.R. DONALDSON.                                *
C*****************************************************************
C
C
C   THE FOLLOWING PROGRAM AND SET OF SUBROUTINES
C   AN INTERACTIVE APPROACH FOR SOLVING EITHER
C   LINEAR OR NONLINEAR LEAST SQUARES PROBLEMS.
C
C   THE PACKAGE CONSISTS OF A MAIN PROGRAM AND 30 
C   SUBROUTINES:  CH, PR, VARPRO, ORFAC1, ORFAC2, 
C   DPA, BACSUB, POSTPR, COV, VARERR, ADA, STAT,
C   ARRAYD, SETUP1, SETUP2, SPLT, SPLTH, PLTMN,
C   PLTERR, ORPLOT, VPLT, VPLTER, VPMAIN, VMAINH, 
C   PLOT1 AND PLOT2.  THE USER IS ASKED TO SUPPLY 
C   5 SUBROUTINES:  CALC, DERIV, WT, TRANS
C   AND BAKTRN.
C
C   THE METHOD USED FOR SOLVING THE LEAST SQUARES 
C   PROBLEM IS DESCRIBED IN THE COMMENTS PRECEDING
C   SUBROUTINE VARPRO.  ALSO INCLUDED ARE A
C   DESCRIPTION OF THE PARAMETERS, COMMENTS ON
C   CONVERGENCE FAILURES AND A LIST OF REFERENCES.
C
C   FOR FURTHER INFORMATION, ONE SHOULD REFER TO
C   THE REPORT "INTERACTIVE VARPRO (INVAR), A NON-
C   LINEAR LEAST SQUARES PROGRAM", BY J.E.RALL
C   (ORAU,UT),AND R.E.FUNDERLIC, (CSD,EXT.4-3148).
C
C
C***********************************************************************
C*
C*
C*
C*                 DRIVER PROGRAM
C*
C*
C*
C***********************************************************************
C
      PARAMETER (MMAX=300, MMAX2=600, MIV=6, MAXN=10, MAXNP2=33,
     2           MAXN1=11, MAXNKG=18, MAXKG=8, MAXP=21 )
      DIMENSION A(MMAX,MAXNP2),T(MMAX,MIV),Y(MMAX),W(MMAX),WEIGHT(MMAX),
     2          ALF(MAXKG),C(MAXN),PH(MAXN1),DPH(MAXN1,MAXKG),ETA(MMAX),
     3          R(MMAX),CR(MAXNKG,MAXNKG),BY(MMAX),BKETA(MMAX),
     4          BT(MMAX,MIV)
      REAL NDI(MAXP)
      DIMENSION ISYM(MMAX2),TPL(MMAX2),ETAPL(MMAX2)
      INTEGER IA(MAXN1),INC(MAXKG,MAXN1)
      INTEGER P
      CHARACTER*1 TITLE(50),NAM(40)
      CHARACTER*1 IANS
      CHARACTER*4 ICODE
      CHARACTER*26 FNAME
      LOGICAL IP,IQ,NC,IWT,ITB,ICV,ICR,WRNBEF,TS,ITRN,ITP,EX,DF
      COMMON /BLOKIP/ IP
      COMMON /BLKNDI/ NDI
      COMMON /BLKINC/ INC
      COMMON /IPIBLK/ IPI
      COMMON /NCBLK/ NC
      COMMON /BLKPH/ PH
      COMMON /BLKDPH/ DPH
      COMMON /BLKTS/ TS
      COMMON /BLKITP/ ITP
      COMMON /BLOCKP/ P
      DATA ITMAX /50/ ,EPS1 /1.0E-5/ ,IPRINT /1/
      DATA IP,IQ,NC,IWT,ITB,ICV,ICR,WRNBEF,TS,ITRN /10*.FALSE./
      DATA ITP /.TRUE./
      DATA IND /4HIND / ,VAR /4HVAR /
      DATA IPR6 /6/ ,IPR10 /10/
C
C************INITIALIZATION SEGMENT OF PROGRAM**************************
C
      PRINT 10
10    FORMAT (/' ENTER TITLE (MAXIMUM OF 50 CHARACTERS):  ')
      READ 20, (TITLE(I),I=1,50)
20    FORMAT (50A1) 
150   PRINT 160
160   FORMAT (' DO YOU WANT LINE PRINTER OUTPUT?',
     1        ' (TYPE "Y" FOR YES)  ')
      READ 170, IANS
170   FORMAT (A1)
      IF (IANS.EQ.'Y') IP=.TRUE.
      IF (IP) WRITE (10,1120) TITLE
      PRINT 180
180   FORMAT (' ARE DATA, WEIGHTS, INITIAL VALUES, ETC. '/
     1        '            ALREADY IN A VARPRO-DEFINED FILE ?  ')
      READ 170, IANS
      IF (IANS.EQ.'Y') GO TO 640
190   PRINT 200
200   FORMAT (' NO. OF DATA POINTS:  ') 
      READ (5,*) M
      IF (IP) WRITE (10,210) M
210   FORMAT (/' NO. OF DATA POINTS: ',I3)
      IF (M.LT.1.OR.M.GT.MMAX) GOTO 190 
      IF (IQ) GO TO 770
220   PRINT 230
230   FORMAT (' NO. OF LINEAR PARAMETERS:  ')
      READ (5,*) N
      IF (N.LT.0.OR.N.GT.MAXN) GOTO 220 
      IF (IP) WRITE (10,240) N
240   FORMAT (' NO. OF LINEAR PARAMETERS: ',I2)
      IF (IQ) GO TO 770
250   PRINT 260
260   FORMAT (' NO. OF NONLINEAR PARAMETERS:  ')
      READ (5,*) KG 
      IF (KG.LT.0.OR.KG.GT.MAXKG) GOTO 250
      IF (IP) WRITE (10,270) KG
270   FORMAT (' NO. OF NONLINEAR PARAMETERS: ',I2)
      IF (IQ) GO TO 770
280   PRINT 290
290   FORMAT (' NO. OF INDEPENDENT VARIABLES:  ') 
      READ (5,*) IV 
      IF (IV.LT.1.OR.IV.GT.MIV) GOTO 280
      IF (IP) WRITE (10,300) IV
300   FORMAT (' NO. OF INDEPENDENT VARIABLES: ',I2)
      IF (IQ) GO TO 770
C
C********************************************************** 
C  THIS SECTION ALLOWS OTHER FILES TO RUN THROUGH THE PROGRAM
C***********************************************************
301   DF=.FALSE.
      PRINT 302
302   FORMAT(' ARE DATA IN A FILE?')
      READ 170,IANS 
      IF(IANS.NE.'Y') GO TO 309
        DF=.TRUE.
303     PRINT 304
304     FORMAT(' WHAT IS THE FILENAME?--NO MORE THAN 26 CHARS'/
     1/'      (TYPE "QUIT" TO EXIT)')
        READ 660, FNAME
        IF (FNAME.EQ.'QUIT') GO TO 301
        INQUIRE (FILE=FNAME,EXIST=EX)
        IF (EX) GO TO 306
          PRINT 305,FNAME
305       FORMAT(' ',A26,' DOES NOT EXIST.')
          GO TO 303 
306     OPEN(UNIT=4,FILE=FNAME,STATUS='OLD')
        DO 307 I=1,M
          READ(4,*) (T(I,J),J=1,IV),Y(I)
307     CONTINUE
        REWIND 4
        CLOSE(UNIT=4,STATUS='KEEP')
       GO TO 370
C*******************************************************
C
  309 CONTINUE
C
      PRINT 310
310   FORMAT (' ARE WEIGHTS TO BE ENTERED WITH DATA?  ')
      READ 170, IANS
      IF (IANS.NE.'Y') GO TO 370
      PRINT 320
      IF (IP) WRITE (10,320)
320   FORMAT (/' DATA POINTS:  LIST INDEPENDENT VARIABLES, DEPENDENT',
     2        /' VARIABLE & WEIGHT = W SQUARED, SEPARATED BY BLANKS', 
     3        ' OR COMMAS'/)
      DO 360 J=1,M
         PRINT 330, J
330      FORMAT (' ',I3,'  ') 
         READ (5,*) (T(J,I),I=1,IV),Y(J),WEIGHT(J)
340      FORMAT (E14.8)
         IF (IP) WRITE (10,350) J,(T(J,I),I=1,IV),Y(J),WEIGHT(J)
350      FORMAT (' ',I3,3X,1P4E14.6,/6X,1P4E14.6) 
360   CONTINUE
      GO TO 460
370   PRINT 380
380   FORMAT (' ARE WEIGHTS TO BE CALCULATED IN SUBROUTINE WT?  ')
      READ 170, IANS
      IF (IANS.EQ.'Y') GO TO 410
      DO 390 J=1,M
390      WEIGHT(J)=1.
      PRINT 400
      IF (IP) WRITE (10,400)
400   FORMAT (' WEIGHTS ARE SET TO 1  '/)
      GO TO 430
410   CONTINUE
      IWT=.TRUE.
      IF (IP) WRITE (10,420)
420   FORMAT (' WEIGHTS CALCULATED IN SUBROUTINE WT '/)
C
  430 CONTINUE
      IF (DF) GO TO 799
C
      PRINT 440
440   FORMAT (/' DATA POINTS:  LIST INDEPENDENT VARIABLES AND'/
     2        ' DEPENDENT VARIABLE, SEPARATED BY BLANKS OR COMMAS'/)
      DO 450 I=1,M
         PRINT 330, I
         READ (5,*) (T(I,J),J=1,IV),Y(I)
         IF (IP) WRITE (10,350) I,(T(I,J),J=1,IV),Y(I)
450   CONTINUE
C
  799 CONTINUE
C
      PRINT 800
800   FORMAT (/' DO DEPENDENT OR INDEPENDENT VARIABLES'/
     2        ' NEED TO BE TRANSFORMED BEFORE ANALYSIS ?')
      READ 170, IANS
      IF (IANS.NE.'Y') GO TO 460
451   ITRN=.TRUE.
      CALL TRANS (T,BT,Y,BY,IV,M,MMAX)
      IF (.NOT.IQ) GOTO 460
      PRINT 459
459   FORMAT (/' IF WEIGHTS WERE CALCULATED IN SUBROUTINE WT,'/
     1        ' THEY MUST BE RECALCULATED WITH TRANSFORMED DATA.')
      GOTO 770
460   IF (IWT) CALL WT (M,WEIGHT,Y)
      IF (KG.EQ.0) THEN
                    IQ = .TRUE.
                    GO TO 770 
                    END IF
      PRINT 470
      IF (IP) WRITE (10,470)
470   FORMAT (/,' INITIAL VALUES OF NONLINEAR PARAMETERS')
      DO 490 J=1,KG 
         PRINT 330, J
         READ (5,*) ALF(J)
         IF (IP) WRITE (10,480) J,ALF(J)
480      FORMAT (' ',I3,3X,1PE13.6)
490   CONTINUE
C
 491  CONTINUE
C
C  DETERMINE INCIDENCE MATRIX 
C
      NP1=N+1
      PH(NP1)=0.
      DO 500 J=1,NP1
         DO 500 K=1,KG
500         INC(K,J)=0
      DO 530 J=1,KG 
505      PRINT 510, J
510      FORMAT (' HOW MANY TERMS OF THE MODEL DEPEND ON ', 
     2           'NONLINEAR PARAMETER ',I2,':  ') 
         READ (5,*) ITERM
         IF (ITERM.LT.1.OR.ITERM.GT.MAXN1) GOTO 505
515      PRINT 520
520      FORMAT (' WHICH TERMS ARE THEY?')
         READ (5,*) (IA(K),K=1,ITERM)
         DO 530 K=1,ITERM
            IB=IA(K)
            IF (IB.LT.1.OR.IB.GT.MAXN1) GOTO 515
            INC(J,IB)=1
530   CONTINUE
      P=0 
      DO 540 J=1,NP1
         DO 540 I=1,KG
540         IF (INC(I,J).EQ.1) P=P+1
      IF (IQ) GO TO 770
      IQ=.TRUE.
C
C  INITIALIZE NUMERICAL DERIVATIVE INCREMENTS
C
      PRINT 550
550   FORMAT (' ALL DERIVATIVES IN CLOSED FORM?  ')
      READ 170, IANS
      IF (IANS.EQ.'Y') GO TO 770
      PRINT 560
560   FORMAT (' ALL NUMERICAL DERIVATIVE INCREMENTS THE SAME?  ')
      READ 170, IANS
      IF (IANS.EQ.'Y') GO TO 600
      PRINT 570
570   FORMAT (' NUMERICAL DERIVATIVE INCREMENTS'/ 
     2        ' (ENTER 0 FOR THOSE TO BE FOUND FROM CLOSED FORM)')
      IF (IP) WRITE (10,580)
580   FORMAT (/' '/' '/' NUMERICAL DERIVATIVE INCREMENTS')
      DO 590 K=1,P
         PRINT 330, K
         READ (5,*) NDI(K)
590      IF (IP) WRITE (10,480) K,NDI(K)
      GO TO 770
600   CONTINUE
      PRINT 610
610   FORMAT (' INCREMENT VALUES:  ')
      READ (5,*) ANS
      IF (IP) WRITE (10,620) ANS
620   FORMAT (' '/' ','INCREMENT VALUES:  ',E15.7)
      DO 630 I=1,P
630      NDI(I)=ANS 
      GO TO 770
640   CONTINUE
C
C  READ IN DATA THAT HAS BEEN SAVED
C
      PRINT 650
650   FORMAT (' TYPE NAME OF FILE (MAX. OF 26 CHARACTERS):  ')
659   READ 660, FNAME
660   FORMAT (A26)
      INQUIRE (FILE=FNAME,EXIST=EX)
      IF (EX) GOTO 662
      PRINT 661,FNAME
661   FORMAT (/' ',A26,' DOES NOT EXIST. TRY AGAIN.')
      GOTO 659
662   OPEN(UNIT=3,FILE=FNAME,STATUS='OLD')
      READ (3,670) (NAM(I),I=1,11),M,(NAM(I),I=1,20),N,(NAM(I),I=1,23),
     2              KG,(NAM(I),I=1,24),IV,(NAM(I),I=1,32)
670   FORMAT (11A1,I4/20A1,I4/23A1,I4/24A1,I4/32A1)
      DO 680 I=1,M
680      READ (3,340) (T(I,J),J=1,IV),Y(I),WEIGHT(I)
      READ (3,690) (NAM(I),I=1,23),(ALF(I),I=1,KG)
690   FORMAT (23A1,8E14.8)
      NP1=N+1
      PH(NP1)=0.
      DO 700 J=1,NP1
         DO 700 K=1,KG
700         INC(K,J)=0
      DO 720 J=1,KG 
         READ (3,710) (NAM(I),I=1,38),(IA(K),K=1,NP1)
710      FORMAT (38A1,11I2)
         DO 720 I=1,NP1
            IF (IA(I).EQ.0) GO TO 720
            IB=IA(I)
            INC(J,IB)=1
720   CONTINUE
      P=0 
      DO 730 J=1,NP1
         DO 730 I=1,KG
730         IF (INC(I,J).EQ.1) P=P+1
      IQ=.TRUE.
      READ (3,740) (NAM(I),I=1,34),(NDI(I),I=1,P) 
740   FORMAT (34A1,5E14.8)
      REWIND 3
      CLOSE(UNIT=3,STATUS='KEEP')
      PRINT 750
750   FORMAT (/' WOULD YOU LIKE YOUR INITIAL',
     2        ' VALUES FILE PRINTED?  ')
      READ 170, IANS
      IF (IANS.NE.'Y') GO TO 790
      CALL PR (1,ICODE,MMAX,M,N,KG,P,IV,T,Y,WEIGHT,ALF)
      GO TO 790
C
C***********************COMMAND POINT*************************
C
770   CONTINUE
      IF (WRNBEF) GO TO 840
C
C     THE FOLLOWING QUESTIONS ARE ASKED ONLY ONCE PER EXECUTION
C
      PRINT 780
780   FORMAT (/' WOULD YOU LIKE TO SAVE THE DATA JUST INPUT?')
      READ 170, IANS
      IF (IANS.EQ.'Y') GO TO 1680
790   PRINT 820
820   FORMAT (/' DOES DATA FOLLOW AN EVENLY SPACED TIME SERIES?  ')
      READ 170, IANS
      IF (IANS.EQ.'Y') TS=.TRUE.
      PRINT 830
830   FORMAT (/' '/' '/'   NOTE:   THE COMMAND "SC" WILL SAVE',
     2        ' YOUR DATA FILE AT ANY TIME')
      WRNBEF=.TRUE. 
840   CONTINUE
      PRINT 850
850   FORMAT (/' NEXT STEP? (TYPE "H" FOR HELP)  ')
      READ 860, ICODE
860   FORMAT (A4)
      IF (ICODE.NE.'H') GO TO 950
      PRINT 870
      IF (IP) WRITE (10,870)
870   FORMAT (/' '/' THE COMMAND CODES FALL INTO FOUR MAIN CATEGORIES'/
     2/'     PRINT-  TO EXAMINE CURRENT VALUES OF INITIALIZED DATA'/
     3 '     CHANGE- TO CORRECT ANY INITIALIZED VALUE, TRANSFORM'/
     4 '             DATA, AND EITHER ADD OR DELETE DATA POINTS'/
     5 '     RUN-    TO RUN THE LEAST SQUARES PROGRAM TO FIND THE'/
     6 '             OPTIMAL PARAMETERS, OR TO RUN A STATISTICAL'/
     7 '             ANALYSIS OF THE DATA RETURNED'/
     8 '     SAVE/   TO SAVE CURRENT DATA OR '/
     9 '     STOP-   TO STOP THE PROGRAM, WITH OR WITHOUT'/ 
     * '             SAVING THE INITIALIZED DATA.'//
     1 ' TYPE "P","C","R", OR "S" FOR THE SPECIFIC COMMANDS IN'/
     2 ' THE ABOVE CATEGORIES.   TYPE "A" TO OBTAIN A LISTING OF'/
     3 ' ALL COMMAND CODES.   ')
      READ 170, IANS
      IF (IANS.EQ.'C') GO TO 890
      IF (IANS.EQ.'R') GO TO 910
      IF (IANS.EQ.'S') GO TO 930
      PRINT 880
      IF (IP) WRITE (10,880)
880   FORMAT (/' '/' '/' COMMAND CODES FOR PRINTING'//
     2 '     PM    PRINT THE NO. OF DATA POINTS'/ 
     3 '     PN    PRINT THE NO. OF LINEAR PARAMETERS'/
     4 '     PKG   PRINT THE NO. OF NONLINEAR PARAMETERS'/
     5 '     PIV   PRINT THE NO. OF INDEPENDENT VARIABLES'/ 
     6 '     PTW   PRINT THE DATA POINTS & WEIGHTS'/
     7 '     PALF  PRINT THE VALUES OF THE NONLINEAR PARAMETERS'/
     8 '     PNDI  PRINT THE NUMERICAL DERIVATIVE INCREMENTS'/
     9 '     PINC  PRINT WHICH TERMS DEPEND ON NONLIN PARAMETERS'/
     * '     LP    TURN ON THE LINE PRINTER FLAG'/
     1 '     NP    TURN OFF THE LINE PRINTER FLAG'/
     2 '     TP    TURN ON THE FLAG TO PLOT ON CRT'/
     3 '     NTP   TURN OFF THE FLAG TO PLOT ON CRT'///)
      IF (IANS.NE.'A') GO TO 770
890   CONTINUE
      PRINT 900
      IF (IP) WRITE (10,900)
900   FORMAT (//' COMMAND CODES FOR CHANGING'//
     2 '     CM    CHANGE NO. OF DATA POINTS'/
     3 '     CN    CHANGE NO. OF LINEAR PARAMETERS'/
     4 '     CKG   CHANGE NO. OF NONLINEAR PARAMETERS'/
     5 '     CIV   CHANGE NO. OF INDEPENDENT VARIABLES'/
     6 '     CT    CHANGE ONE OR MORE DATA POINTS'/
     7 '     CW    CHANGE ONE OR MORE WEIGHTS'/
     8 '     CALF  CHANGE ONE OR MORE NONLINEAR PARAMETERS'/
     9 '     CNDI  CHANGE ONE OR MORE NUMER DERIV INCREMENTS'/
     * '     CTRN  TRANSFORM THE DATA'/ 
     * '     CINC  CHANGE WHICH TERMS DEPEND ON NONLIN PARAMETERS'///)
      IF (IANS.NE.'A') GO TO 770
910   CONTINUE
      PRINT 920
      IF (IP) WRITE (10,920)
920   FORMAT (//' COMMAND CODES FOR RUNNING'//
     2 '     RV    RUN THE LEAST SQUARES PROGRAM'/
     3 '     RS    RUN A STATISTICAL ANALYSIS OF RESULTS'///)
      IF (IANS.NE.'A') GO TO 770
930   CONTINUE
      PRINT 940
      IF (IP) WRITE (10,940)
940   FORMAT (//' COMMAND CODES FOR SAVING DATA FILE AND/OR STOPPING'//
     2 '     SC    SAVE CURRENT VALUES OF DATA & CONTINUE '/
     3 '     SD    SAVE CURRENT VALUES OF DATA & STOP'/
     4 '     SR    SAVE CURRENT VALUES OF INDEP. VARIABLES ',
     5 'AND RESIDUALS & CONTINUE'/'     S0    STOP WITHOUT SAVING'///)
      GO TO 770
950   CONTINUE
      IF (ICODE.EQ.'RV') GO TO 960
      IF (ICODE.EQ.'RS') GO TO 1580
      IF ((ICODE.EQ.'PM').OR.(ICODE.EQ.'PN').OR.(ICODE.EQ.'PKG').OR.
     2    (ICODE.EQ.'PIV').OR.(ICODE.EQ.'PTW').OR.(ICODE.EQ.'PALF').OR.
     3    (ICODE.EQ.'PNDI').OR.(ICODE.EQ.'LP').OR.(ICODE.EQ.'NP').OR. 
     4    (ICODE.EQ.'PINC').OR.(ICODE.EQ.'TP').OR.(ICODE.EQ.'NTP'))
     5    GOTO 1600 
      IF (ICODE.EQ.'CM') GO TO 190
      IF (ICODE.EQ.'CN') GO TO 220
      IF (ICODE.EQ.'CKG') GO TO 250
      IF (ICODE.EQ.'CIV') GO TO 280
      IF (ICODE.EQ.'CTRN') GOTO 451
      IF (ICODE.EQ.'CINC') GO TO 491
      IF ((ICODE.EQ.'CT').OR.(ICODE.EQ.'CW').OR.(ICODE.EQ.'CALF').OR. 
     2    (ICODE.EQ.'CNDI')) GO TO 1590 
      IF (ICODE.EQ.'SC') GO TO 1680
      IF (ICODE.EQ.'SD') GO TO 1680
      IF (ICODE.EQ.'SR') GO TO 1630
      IF (ICODE.EQ.'S0') GO TO 1790
      GO TO 1610
C
C*********** RUN THE LEAST SQUARES PROGRAM ***********
C
960   CONTINUE
      IF (KG.EQ.0) GO TO 1090 
      IPRINT = 1
970   PRINT 980
980   FORMAT (/' DO YOU WISH TO CHANGE ANY DEFAULTS ',
     2        'IN THE LEAST SQUARES PROGRAM?'/
     3        '        (50: MAX. NO. OF ITERATIONS)'/
     4        '        (1.0E-5: CONVERG. TOLERANCE)'/
     5        '        (1: OUTPUT PRINTING INDEX)')
      READ 170, IANS
      IF (IANS.NE.'Y') GO TO 1110
990   PRINT 1000
1000  FORMAT (/' COMMAND CODES FOR CHANGING DEFAULTS'/
     2        '     CMNI  CHANGE MAX. NO. ITERATIONS'/
     3        '     CCT   CHANGE CONVERG. TOLERANCE'/
     4        '     COPI  CHANGE OUTPUT PRINTING INDEX')
      READ 860, ICODE
      IF (ICODE.EQ.'CMNI') GO TO 1010
      IF (ICODE.EQ.'CCT') GO TO 1050
      IF (ICODE.EQ.'COPI') GO TO 1070
      GO TO 970
1010  PRINT 1020
1020  FORMAT (/' THE NUMBER OF ITERATIONS MUST BE GREATER '/' ',
     2        'THAN 0 AND LESS THAN OR EQUAL TO 250.'/' ',
     3        'HOW MANY ITERATIONS?         ')
      READ (5,*) ITMAX
      IF (ITMAX.LE.0.OR.ITMAX.GT.250) GO TO 1010
1030  PRINT 1040
1040  FORMAT (/' DO YOU WISH TO CHANGE ANY OTHERS?')
      READ 170, IANS
      IF (IANS.NE.'Y') GO TO 1110
      GO TO 990
1050  PRINT 1060
1060  FORMAT (/' THE CONVERGENCE CRITERION MUST NOT BE LESS ',
     2        'THAN 1.5E-8. '/' ',
     3        'TYPE IN THE NEW CONVERGENCE CRITERION.    ') 
      READ (5,*) EPS1
      IF (EPS1.LT.1.5E-8) GO TO 1050
      GO TO 1030
1070  PRINT 1080
1080  FORMAT (/' TYPE OUTPUT PRINTING INDEX: ')
      READ (5,*) IPRINT
1090  IF (KG.EQ.0) IPRINT=0
      IF (IP) WRITE (10,1100) IPRINT
1100  FORMAT (' '/' ','OUTPUT PRINTING INDEX:',I3)
      IF (KG.NE.0) GO TO 1030 
1110  CONTINUE
      PRINT 1120, TITLE
      IF (IP) WRITE (10,1120) TITLE
1120  FORMAT ('1',20X,50A1//) 
      CALL VARPRO (N,KG,M,MMAX,N+P+2,IV,T,Y,W,WEIGHT,A,IPRINT,ALF,C,
     2             IERR,ITMAX,EPS1)
      IF (IERR.LE.-4) GO TO 770
      NC=.FALSE.
      NPKG=N+KG
C
C  DETERMINE WHICH TABLES ARE TO BE WRITTEN
C
      PRINT 1130
1130  FORMAT (' INDICATE WHICH OF THE FOLLOWING SHOULD BE PRINTED'/
     2        ' BY TYPING "Y": '/'   TABLE OF',
     3        ' OBSERVED & PREDICTED VALUES, WTD RESIDUAL: ')
      READ 170, IANS
      IF (IANS.EQ.'Y') ITB=.TRUE.
      PRINT 1140
1140  FORMAT ('   COVARIANCE MATRIX:  ')
      READ 170, IANS
      IF (IANS.EQ.'Y') ICV=.TRUE.
      PRINT 1150
1150  FORMAT ('   CORRELATION MATRIX:  ')
      READ 170, IANS
      IF (IANS.EQ.'Y') ICR=.TRUE.
      IF (.NOT.ICV) GO TO 1170
C
C  PRINT TABLE OF COVARIANCE MATRIX
C
      ICV=.FALSE.
      PRINT 1160
      IF (IP) WRITE (10,1160) 
1160  FORMAT ('1'/' '/' '/15X,'COVARIANCE MATRIX'/)
      CALL ARRAYD (A,NPKG,NPKG,MMAX)
1170  DO 1180 I=1,M 
         II=I
1180     CALL CALC (II,KG,N,MMAX,IV,C,ALF,T,ETA(I))
      IF (.NOT.ITB) GO TO 1250
C
C  PRINT TABLE OF OBSERVED AND PREDICTED VALUES, ETC.
C
      ITB=.FALSE.
      PRINT 1190
      IF (IP) WRITE (10,1190) 
1190  FORMAT ('1'/' '/' '/
     1         ' INDEX',2X,'IND VAR 1',4X,'OBS VAL',6X,'PRED VAL',5X, 
     2        'WTD RESID',5X,'WEIGHTS'/)
      DO 1210 I=1,M 
         PRINT 1200, I,T(I,1),Y(I),ETA(I),A(I,NPKG+1),WEIGHT(I)
         IF (IP) WRITE (10,1200) I,T(I,1),Y(I),ETA(I),A(I,NPKG+1),
     2                           WEIGHT(I)
1200     FORMAT (' ',I3,1P5E13.5)
1210  CONTINUE
C
1250  IF (.NOT.ICR) GO TO 1280
C
C  PRINT TABLE OF CORRELATION MATRIX
C
      ICR=.FALSE.
      DO 1260 I=1,NPKG
         DO 1260 J=I,NPKG
            CR(I,J)=A(I,J)/(SQRT(A(I,I))*SQRT(A(J,J)))
1260        CR(J,I)=CR(I,J)
      PRINT 1270
      IF (IP) WRITE (10,1270) 
1270  FORMAT ('1'/' '/' '/15X,' CORRELATION MATRIX')
      CALL ARRAYD (CR,NPKG,NPKG,MAXNKG) 
1280  DO 1290 I=1,M 
1290     R(I)=Y(I)-ETA(I)
      PRINT 1300
1300  FORMAT (/' PLOT ROUTINE CURRENTLY SET NOT TO DO ANY PLOTS.'/
     2        ' DO YOU WANT TO START PLOTTING? ') 
      READ 170, IANS
      IF (IANS.NE.'Y') GO TO 770
C
C  DETERMINE WHICH PLOTS ARE TO BE PRINTED
C
      IF (IV.NE.1) GO TO 1310 
      INV=1
      GO TO 1330
1310  CONTINUE
      PRINT 1320
1320  FORMAT (/' WHICH INDEPENDENT VARIABLE SHOULD BE USED?'/
     2        ' (TYPE ITS NUMBER)  ')
      READ (5,*) INV
1330  CONTINUE
1340  CONTINUE
1341  PRINT 1350
1350  FORMAT (/' IF YOU WANT TO PLOT RESIDUALS, TYPE "1"'/
     2        ' IF YOU WANT TO PLOT OBSERVED & PREDICTED VALUES, ',
     3        'TYPE "2"'/' IF YOU WANT TO PLOT BOTH, TYPE "3":  ')
      READ (5,*) IPI
      IF (IPI.NE.1.AND.IPI.NE.2.AND.IPI.NE.3) GOTO 1341
      IF (IPI.EQ.2) GO TO 1390
C
C  PLOT RESIDUALS
C
      IF (ITP) PRINT 1360
1360  FORMAT ('1',15X,'RESIDUALS')
      IF (IP) WRITE (10,1360) 
      IF (ITP) PRINT 1370, INV
1370  FORMAT (' ',15X,'IND VAR',I2)
      IF (IP) WRITE (10,1370) INV
      IF (TS.AND.ITP) CALL VPLT (R,M,IPR6,0)
      IF (TS.AND.IP) CALL VPLT (R,M,IPR10,0)
      IF (.NOT.TS) CALL SETUP1 (T,M,INV,R,ISYM,TPL,ETAPL)
      IF (IPI.EQ.1) GOTO 1550 
C
C  PLOT OBS. AND PRED. VALUES 
C
1390  IF (.NOT.ITRN.AND.ITP) PRINT 1400 
1400  FORMAT ('1',15X,'OBSERVED (*) & PREDICTED (+) VALUES')
      IF (IP.AND..NOT.ITRN) WRITE (10,1400)
      IF (ITRN.AND.ITP) PRINT 1410
1410  FORMAT ('1',15X,'OBSERVED (*) & PREDICTED ',
     2        '(+) VALUES (TRANSFORMED DATA)')
      IF (IP.AND.ITRN) WRITE (10,1410)
      IF (ITP) PRINT 1370, INV
      IF (IP) WRITE (10,1370) INV
      IF (TS.AND.ITP) CALL VPLT2 (Y,ETA,M,IPR6,1) 
      IF (TS.AND.IP) CALL VPLT2 (Y,ETA,M,IPR10,1) 
      IF (.NOT.TS) CALL SETUP2 (T,M,INV,Y,ETA,ISYM,TPL,ETAPL)
      IF (.NOT.ITRN) GOTO 1550
C
C  PLOT BACK TRANSFORMED DATA 
C
      CALL BAKTRN (ETA,BKETA,M)
      IF (ITP) PRINT 1450
1450  FORMAT ('1',15X,' OBSERVED (*) & PREDICTED ',
     2        '(+) VALUES (BACK-TRANSFORMED DATA)')
      IF (IP) WRITE (10,1450) 
      IF (ITP) PRINT 1370, INV
      IF (IP) WRITE (10,1370) INV
      IF (TS.AND.ITP) CALL VPLT2 (BY,BKETA,M,IPR6,1)
      IF (TS.AND.IP) CALL VPLT2 (BY,BKETA,M,IPR10,1)
      IF (.NOT.TS) CALL SETUP2 (BT,M,INV,BY,BKETA,ISYM,TPL,ETAPL)
1480  PRINT 1490
1490  FORMAT (/' DO YOU WISH TO PRINT A TABLE OF THE ',
     2        'BACK-TRANSFORMED VALUES?'/)
      READ 170, IANS
      IF (IANS.NE.'Y') GO TO 1550
C
C  PRINT TABLE OF BACK-TRANSFORMED DATA 
C
      PRINT 1500
1500  FORMAT ('1'/'  OBSERVED AND PREDICTED VALUES BACK-TRANSFORMED'/)
      IF (IP) WRITE (10,1500) 
      PRINT 1510
1510  FORMAT (' ','INDEX',2X,'IND VAR 1',4X,'OBS VAL',6X,'PRED VAL'/) 
      IF (IP) WRITE (10,1510) 
      DO 1530 I=1,M 
         PRINT 1520, I,BT(I,1),BY(I),BKETA(I)
1520     FORMAT (' ',I3,1P3E13.5)
         IF (IP) WRITE (10,1520) I,BT(I,1),BY(I),BKETA(I)
1530  CONTINUE
C
1550  IF (IV .EQ. 1) GOTO 770 
      PRINT 1560
1560  FORMAT (/'  CONTINUE PLOTTING USING A DIFFERENT'/
     1        ' INDEPENDENT VARIABLE ?   ')
      READ 170, IANS
      IF (IANS.NE.'Y') GO TO 770
      GO TO 1310
C
C*********** RUN THE STATISTICAL PACKAGE **********
C
1580  CONTINUE
      NKG2=N+KG+2
      CALL STAT (M,N,NKG2,KG,Y,WEIGHT,C,ALF,A,IV,T,ETA,R,ISYM,TPL,ETAPL)
      GO TO 770
C
C*********** CALL SUBROUTINE TO CHANGE ********** 
C
1590  CONTINUE
      CALL CH (ICODE,MMAX,M,KG,IV,T,Y,WEIGHT,ALF) 
      GO TO 770
C
C*********** CALL SUBROUTINE TO PRINT *********** 
C
1600  CONTINUE
      CALL PR (99999,ICODE,MMAX,M,N,KG,P,IV,T,Y,WEIGHT,ALF) 
      GO TO 770
C
C*********** PRINT COMMAND CODE ERROR MESSAGE **********
C
1610  CONTINUE
      PRINT 1620
1620  FORMAT (' COMMAND CODE WAS INVALID'//)
      GO TO 770
C
C********** SAVE RESIDUALS AND CONTINUE********** 
C
1630  CONTINUE
      PRINT 1640
1640  FORMAT (' TYPE NAME OF FILE FOR STORING RESIDUAL DATA'/
     1        ' (MAXIMUM OF 26 CHARACTERS):  ')
1641  READ 660, FNAME
      INQUIRE (FILE=FNAME,EXIST=EX)
      IF (.NOT.EX) GOTO 1643
      PRINT 1642,FNAME
1642  FORMAT (/' ',A26,' ALREADY EXISTS.  USE ANOTHER NAME.')
      GOTO 1641
1643  OPEN(UNIT=3,FILE=FNAME,STATUS='NEW')
      DO 1660 I=1,M 
         WRITE (3,1650) (T(I,J),J=1,IV),R(I)
1650     FORMAT (7(E14.8,2X)) 
1660  CONTINUE
      REWIND 3
      CLOSE(UNIT=3,STATUS='KEEP')
      IF (IP) WRITE (10,1670) FNAME
1670  FORMAT (//' RESIDUAL DATA IS STORED IN FILE:  ',A26)
      GO TO 770
C
C*********** SAVE CURRENT VALUES OF DATA ***********
C
1680  CONTINUE
      PRINT 1690
1690  FORMAT (' TYPE NAME OF FILE FOR STORING DATA'/
     2        ' (MAXIMUM OF 26 CHARACTERS):   ')
1691  READ 660, FNAME
      INQUIRE (FILE=FNAME,EXIST=EX)
      IF (.NOT.EX) GOTO 1693
      PRINT 1692,FNAME
1692  FORMAT (/' ',A26,' ALREADY EXISTS.  USE ANOTHER NAME.')
      GOTO 1691
1693  OPEN(UNIT=3,FILE=FNAME,STATUS='NEW')
      WRITE (3,1700) M,N,KG,IV
1700  FORMAT (' DATA PTS: ',I4/' LINEAR PARAMETERS: ',I4/
     2        ' NONLINEAR PARAMETERS: ',I4/' INDEPENDENT VARIABLES: ',
     3        I4/' INDEP. VAR., DEP. VAR. & WEIGHTS: ')
      DO 1710 I=1,M 
1710     WRITE (3,340) (T(I,J),J=1,IV),Y(I),WEIGHT(I)
      WRITE (3,1720) (ALF(I),I=1,KG)
1720  FORMAT (' NONLINEAR PARAMETERS: ',8E14.8)
      DO 1760 I=1,KG
         DO 1730 K=1,NP1
1730        IA(K)=0.
         K=1
         DO 1740 J=1,NP1
            IF (INC(I,J).NE.1) GO TO 1740
            IA(K)=J 
            K=K+1
1740     CONTINUE
         WRITE (3,1750) I,(IA(K),K=1,NP1)
1750     FORMAT (' TERMS WHICH DEPEND ON NONLIN PARM',I2,':  ',11I2)
1760  CONTINUE
      WRITE (3,1770) (NDI(I),I=1,P)
1770  FORMAT (' NUMERICAL DERIVATIVE INCREMENTS: ',5E14.8)
      REWIND 3
      CLOSE(UNIT=3,STATUS='KEEP')
      IF (IP) WRITE (10,1780) FNAME
1780  FORMAT (//' DATA IS STORED IN FILE:  ',A26/)
      IF (.NOT.WRNBEF) GO TO 790
      IF (ICODE.EQ.'SC') GO TO 770
C
C*********** STOP *********** 
C
1790  CONTINUE
      STOP
      END 
C
C
      SUBROUTINE CH (ICODE,MMAX,M,KG,IV,T,Y,WEIGHT,ALF)
C
C***********************************************************************
C
C   CHANGE THE VALUE OF DATA POINTS, WEIGHTS, NONLINEAR
C   PARAMETERS OR NUMERICAL DERIVATIVE INCREMENTS 
C   (DEPENDING ON THE COMMAND CODE ENTERED).
C
      PARAMETER (MAXP=21)
      DIMENSION T(MMAX,IV),Y(M),WEIGHT(M),ALF(KG) 
      REAL NDI(MAXP)
      LOGICAL IP
      CHARACTER*1 IANS
      CHARACTER*4 ICODE
      COMMON /BLKNDI/ NDI
      COMMON /BLOKIP/ IP
C
10    FORMAT (A1)
      IF (ICODE.EQ.'CT') GO TO 20
      IF (ICODE.EQ.'CW') GO TO 80
      IF (ICODE.EQ.'CALF') GO TO 120
      IF (ICODE.EQ.'CNDI') GO TO 160
C
C ********** CHANGE ONE OR MORE DATA POINTS **********
C
20    CONTINUE
      PRINT 30
30    FORMAT (' WHICH DATA POINT IS TO BE CHANGED?')
      PRINT 40
40    FORMAT (' (TYPE ITS INDEX NO. AS IT APPEARS IN THE LIST)  ')
      READ (5,*) IL 
      IF (IL.LT.1.OR.IL.GT.M) GOTO 20
      PRINT 50, IL
50    FORMAT (' TYPE NEW INDEP & DEP VARIABLES FOR DATA POINT',I3,': ')
      READ (5,*) (T(IL,J),J=1,IV),Y(IL) 
      IF (IP) WRITE (10,60) IL,(T(IL,J),J=1,IV),Y(IL)
60    FORMAT (' '/' NEW INDEP & DEPENDENT VARIABLES FOR POINT',I3,': '/
     2        1P4E15.6/15X,1P3E15.6)
      PRINT 70
70    FORMAT (' DO YOU WANT TO CHANGE ANY OTHERS?  ')
      READ 10, IANS 
      IF (IANS.EQ.'Y') GO TO 20
      GO TO 200
C
C ********** CHANGE ONE OR MORE WEIGHTS **********
C
80    CONTINUE
      PRINT 79
79    FORMAT (' ARE WEIGHTS TO BE CHANGED USING SUBROUTINE WT?')
      READ 10,IANS
      IF (IANS .NE. 'Y') GOTO 84
      CALL WT (M,WEIGHT,Y)
      GOTO 200
84    PRINT 81
81    FORMAT (' ARE ALL WEIGHTS TO BE CHANGED TO THE SAME VALUE?')
      READ 10,IANS
      IF (IANS .NE. 'Y') GOTO 85
      PRINT 82
82    FORMAT (' ENTER THE VALUE:   ')
      READ (5,*) WTVAL
      DO 83 I=1,M
         WEIGHT(I) = WTVAL
83    CONTINUE
      IF (IP) WRITE (10,101) WTVAL
101   FORMAT (' '/' NEW WEIGHT VALUE FOR ALL DATA POINTS:  ',1PE15.6) 
      GOTO 200
85    PRINT 90
90    FORMAT (' WHICH WEIGHT IS TO BE CHANGED?')
      PRINT 40
      READ (5,*) IL 
      IF (IL.LT.1.OR.IL.GT.M) GOTO 85
      PRINT 100, IL 
100   FORMAT (' ENTER NEW WEIGHT FOR DATA POINT',I4,': ')
      READ (5,*) WEIGHT(IL)
      IF (IP) WRITE (10,110) IL,WEIGHT(IL)
110   FORMAT (' '/' NEW WEIGHT FOR DATA POINT',I4,': ',1PE15.6)
      PRINT 70
      READ 10, IANS 
      IF (IANS.EQ.'Y') GO TO 80
      GO TO 200
C
C ******** CHANGE ONE OR MORE NONLINEAR PARAMETERS *********
C
120   CONTINUE
      PRINT 130
130   FORMAT (' WHICH NONLINEAR PARAMETER IS TO BE CHANGED?')
      PRINT 40
      READ (5,*) IL 
       IF (IL.LT.1.OR.IL.GT.KG) GOTO 120
      PRINT 140, IL 
140   FORMAT (' ENTER NEW VALUE FOR PARAMETER',I4,': ')
      READ (5,*) ALF(IL)
      IF (IP) WRITE (10,150) IL,ALF(IL) 
150   FORMAT (' '/' NEW VALUE FOR NONLINEAR PARAMETER',I4,': ',1PE15.6)
      PRINT 70
      READ 10, IANS 
      IF (IANS.EQ.'Y') GO TO 120
      GO TO 200
C
C ******* CHANGE ONE OR MORE NUMER DERIV INCREMENTS ********
C
160   CONTINUE
      PRINT 170
170   FORMAT (' WHICH NUMER. DERIV. INCREMENT IS TO BE CHANGED?')
      PRINT 40
      READ (5,*) IL 
      IF (IL.LT.1.OR.IL.GT.MAXP) GOTO 160
      PRINT 180, IL 
180   FORMAT (' NEW INCREMENT FOR DERIVATE',I4,': ')
      READ (5,*) NDI(IL)
      IF (IP) WRITE (10,190) IL,NDI(IL) 
190   FORMAT (' '/' NEW INCREMENT FOR DERIVATIVE',I4,': ',1PE15.6)
      PRINT 70
      READ 10, IANS 
      IF (IANS.EQ.'Y') GO TO 160
200   RETURN
C
      END 
C***********************************************************************
C
      SUBROUTINE PR (NUM,ICODE,MMAX,M,N,KG,P,IV,T,Y,WEIGHT,ALF)
C
C***********************************************************************
C
C   PRINT THE VALUES OF ANY OR ALL OF THE INITIALIZED
C   PARAMETERS.  THIS SUBROUTINE IS CALLED IN
C   RESPONSE TO ONE OF THE "PRINT" COMMANDS OR
C   WHEN THE INITIALIZED DATA IS OBTAINED FROM A
C   SAVED FILE.
C
      PARAMETER (MAXN1=11, MAXKG=8, MAXP=21)
      DIMENSION T(MMAX,IV),Y(M),WEIGHT(M),ALF(KG) 
      REAL NDI(MAXP)
      INTEGER INC(MAXKG,MAXN1),IA(MAXN1),P
      CHARACTER*1 IANS
      CHARACTER*4 ICODE,IND,VAR
      LOGICAL IP,ITP
      COMMON /BLKINC/ INC
      COMMON /BLKNDI/ NDI
      COMMON /BLOKIP/ IP
      COMMON /BLKITP/ ITP
      DATA IND /'IND '/ ,VAR /'VAR '/
210   FORMAT (A1)
220   FORMAT (' (TYPE ITS INDEX NO. AS IT APPEARS IN THE LIST)')
C
      IF (NUM.EQ.1) GO TO 640 
      IF (ICODE.EQ.'PM') GO TO 230
      IF (ICODE.EQ.'PN') GO TO 250
      IF (ICODE.EQ.'PKG') GO TO 270
      IF (ICODE.EQ.'PINC') GO TO 290
      IF (ICODE.EQ.'PIV') GO TO 340
      IF (ICODE.EQ.'PTW') GO TO 360
      IF (ICODE.EQ.'PALF') GO TO 500
      IF (ICODE.EQ.'PNDI') GO TO 590
      IF (ICODE.EQ.'LP') GO TO 620
      IF (ICODE.EQ.'NP') GO TO 630
      IF (ICODE.EQ.'TP') GO TO 632
      IF (ICODE.EQ.'NTP') GO TO 634
C
C ********** PRINT THE NO. OF DATA POINTS **********
C
230   CONTINUE
      PRINT 240, M
      IF (IP) WRITE (10,240) M
240   FORMAT (' '/' NO. OF DATA POINTS:   ',I4)
      GO TO 700
C
C ********** PRINT THE NO. OF LINEAR PARAMETERS **********
C
250   CONTINUE
      PRINT 260, N
      IF (IP) WRITE (10,260) N
260   FORMAT (' '/' NO. OF LINEAR PARAMETERS:   ',I4)
      GO TO 700
C
C ********* PRINT THE NO. OF NONLINEAR PARAMETERS **********
C
270   CONTINUE
      PRINT 280, KG 
      IF (IP) WRITE (10,280) KG
280   FORMAT (/' '/' NO. OF NONLINEAR PARAMETERS:   ',I4)
      GO TO 700
C
C ******** PRINT WHICH TERMS DEPEND ON NONLIN PARM *********
C
290   CONTINUE
      NP1=N+1
      DO 330 I=1,KG 
         K=1
         DO 300 J=1,NP1
300         IA(J)=0 
         DO 310 J=1,NP1
            IF (INC(I,J).NE.1) GO TO 310
            IA(K)=J 
            K=K+1
310      CONTINUE
         PRINT 320, I,(IA(K),K=1,NP1)
320    FORMAT (/' '/' TERMS WHICH DEPEND ON NONLIN PARM',I2,':  ',11I2)
330      IF (IP) WRITE (10,320) I,(IA(K),K=1,NP1) 
      GO TO 700
C
C ********* PRINT THE NO. OF INDEPENDENT VARIABLES *********
C
340   CONTINUE
      PRINT 350, IV 
      IF (IP) WRITE (10,350) IV
350   FORMAT (/' '/' NO. OF INDEPENDENT VARIABLES:   ',I4)
      GO TO 700
C
C ********** PRINT THE DATA POINTS & WEIGHTS **********
C
360   CONTINUE
      PRINT 370
370   FORMAT (' ALL OF THE DATA POINTS?  ')
      READ 210, IANS
      IF (IANS.EQ.'Y') GO TO 430
380   CONTINUE
      PRINT 390
390   FORMAT (' WHICH DATA POINT?  ')
      PRINT 220
      READ (5,*) IL 
      PRINT 400, IL 
      IF (IP) WRITE (10,400) IL
400   FORMAT (/' '/' FOR POINT',I4,'VALUES OF INDEP. VAR,', 
     2        ' DEP VAR & WT ARE')
      PRINT 410, (T(IL,J),J=1,IV),Y(IL),WEIGHT(IL)
      IF (IP) WRITE (10,410) (T(IL,J),J=1,IV),Y(IL),WEIGHT(IL)
410   FORMAT (' ',1P4E15.6/15X,1P4E15.6,/)
      PRINT 420
420   FORMAT (' DO YOU WANT TO EXAMINE ANY OTHERS?  ')
      READ 210, IANS
      IF (IANS.EQ.'Y') GO TO 380
      GO TO 700
430   CONTINUE
      IVP=MIN0(IV,2)
      PRINT 440, (IND,VAR,J,J=1,IVP)
      IF (IP) WRITE (10,440) (IND,VAR,J,J=1,IVP)
440   FORMAT (/' '/' INDEX',5X,'WEIGHT',9X,'DEP VAR',8X,2(2A4,I2,5X)/)
      DO 460 K=1,M
         PRINT 450, K,WEIGHT(K),Y(K),(T(K,J),J=1,IVP)
         IF (IP) WRITE (10,450) K,WEIGHT(K),Y(K),(T(K,J),J=1,IVP)
450      FORMAT (' ',I3,3X,1P4E15.6)
460   CONTINUE
      IF (IV.LE.2) GO TO 700
      PRINT 470, (IND,VAR,J,J=3,IV)
470   FORMAT (/' INDEX',3X,4(2A4,I2,5X)/)
      IF (IP) WRITE (10,470) (IND,VAR,J,J=3,IV)
      DO 490 K=1,M
         PRINT 480, K,(T(K,J),J=3,IV)
         IF (IP) WRITE (10,480) K,(T(K,J),J=3,IV) 
480      FORMAT (' ',I3,3X,1P4E15.6)
490   CONTINUE
      GO TO 700
C
C ****** PRINT THE VALUES OF THE NONLINEAR PARAMETERS ******
C
500   CONTINUE
      PRINT 510
510   FORMAT (' ALL THE NONLINEAR PARAMETERS?  ') 
      READ 210, IANS
      IF (IANS.EQ.'Y') GO TO 550
520   CONTINUE
      PRINT 530
530   FORMAT (' WHICH PARAMETER?  ')
      PRINT 220
      READ (5,*) IL 
      PRINT 540, IL,ALF(IL)
      IF (IP) WRITE (10,540) IL,ALF(IL) 
540   FORMAT (/' '/' NONLINEAR PARAMETER',I4,' IS',1PE15.6) 
      PRINT 420
      READ 210, IANS
      IF (IANS.EQ.'Y') GO TO 520
      GO TO 700
550   CONTINUE
      PRINT 560
      IF (IP) WRITE (10,560)
560   FORMAT (/' '/' NONLINEAR PARAMETERS:   ')
      DO 580 K=1,KG 
         PRINT 570, K,ALF(K)
         IF (IP) WRITE (10,570) K,ALF(K)
570      FORMAT (' ',I3,3X,1PE15.6)
580   CONTINUE
      GO TO 700
C
C ********** PRINT THE NUMER DERIV INCREMENTS **********
C
590   CONTINUE
      PRINT 600
      IF (IP) WRITE (10,600)
600   FORMAT (/' '/' NUMERICAL DERIVATIVE INCREMENTS:')
      DO 610 J=1,P
         PRINT 570, J,NDI(J)
         IF (IP) WRITE (10,570) J,NDI(J)
610   CONTINUE
      GO TO 700
C
C ********** TURN ON THE LINE PRINTER FLAG **********
C
620   CONTINUE
      IP=.TRUE.
      GO TO 700
C
C ********** TURN OFF THE LINE PRINTER FLAG **********
C
630   CONTINUE
      IP=.FALSE.
      GO TO 700
C
C ********** TURN ON THE CRT PLOTTING FLAG **********
C
632   CONTINUE
      ITP = .TRUE.
      GOTO 700
C
C ********** TURN OFF THE CRT PLOTTING FLAG **********
C
634   CONTINUE
      ITP = .FALSE. 
      GOTO 700
C
C ********** PRINT ALL INITIAL VALUES **********
C
640   CONTINUE
      PRINT 240, M
      IF (IP) WRITE (10,240) M
      PRINT 260, N
      IF (IP) WRITE (10,260) N
      PRINT 280, KG 
      IF (IP) WRITE (10,280) KG
      PRINT 350, IV 
      IF (IP) WRITE (10,350) IV
      IVP=MIN0(IV,2)
      PRINT 440, (IND,VAR,J,J=1,IVP)
      IF (IP) WRITE (10,440) (IND,VAR,J,J=1,IVP)
      DO 650 K=1,M
         PRINT 450, K,WEIGHT(K),Y(K),(T(K,J),J=1,IVP)
650      IF (IP) WRITE (10,450) K,WEIGHT(K),Y(K),(T(K,J),J=1,IVP)
      IF (IV.LE.2) GO TO 670
      PRINT 470, (IND,VAR,J,J=3,IV)
      IF (IP) WRITE (10,470) (IND,VAR,J,J=3,IV)
      DO 660 K=1,M
         PRINT 480, K,(T(K,J),J=3,IV)
         IF (IP) WRITE (10,480) K,(T(K,J),J=3,IV) 
660   CONTINUE
670   IF (KG.EQ.0) GO TO 700
      PRINT 560
      IF (IP) WRITE (10,560)
      DO 680 K=1,KG 
         PRINT 570, K,ALF(K)
         IF (IP) WRITE (10,570) K,ALF(K)
680   CONTINUE
      PRINT 600
      IF (IP) WRITE (10,600)
      DO 690 K=1,P
         PRINT 570, K,NDI(K)
         IF (IP) WRITE (10,570) K,NDI(K)
690   CONTINUE
      GO TO 290
700   RETURN
C
      END 
C***********************************************************************
C
      SUBROUTINE ARRAYD (BA,NRR,NCC,LVV)
C
C***********************************************************************
C
C   PROGRAM AUTHORS  L. FELDMAN AND J DUNPHEY FORMERLY
C   WITH COMPUTING SCIENCES DIVISION, UNION CARBIDE,
C   NUCLEAR DIVISION, OAK RIDGE
C
C
C  ARRAYD IS AN OUTPUT SUBROUTINE WHICH PRINTS A
C  SINGLE PRECISION MATRIX.  THE CALLING SEQUENCE 
C  IS
C      CALL ARRAYD (A,NR,NC,ID)
C  WHERE
C   A      IS A SINGLE PRECISION ARRAY TO BE
C          LISTED.
C   NR     IS THE INTEGER NUMBER OF ROWS OF THE
C          MATRIX A DESIRED IN THE OUTPUT.
C   NC     IS THE INTEGER NUMBER OF COLUMNS OF THE
C          MATRIX A DESIRED IN THE OUTPUT.
C   ID     IS THE INTEGER FIRST DIMENSION OF A IN 
C          THE CALLING PROGRAM, I.E. A(ID,*) OR A(ID)
C
      CHARACTER*8 COLUMN
      DIMENSION BA(1)
      LOGICAL IP
      COMMON /BLOKIP/ IP
      DATA COLUMN /'COLUMN  '/
C
      DO 1801 ML=1,NCC,4
         MH=ML+3
         IF (NCC.GE.MH) GO TO 1760
         MH=NCC
1760     MLM1=(ML-1)*LVV
         MHM1=(MH-1)*LVV
         DO 1801 LI=1,NRR,50
            NL=LI
            NH=NL+49
            IF (NRR.GE.NH) GO TO 1770
            NH=NRR
1770        PRINT 1780, (COLUMN,J,J=ML,MH)
            IF (IP) WRITE (10,1780) (COLUMN,J,J=ML,MH)
1780        FORMAT (/' '/' ',4H ROW,6X,4(A6,I5,4X)/)
            DO 1790 L=NL,NH
               MML=MLM1+L
               MMH=MHM1+L
               PRINT 1800, L,(BA(J),J=MML,MMH,LVV)
1790           IF (IP) WRITE (10,1800) L,(BA(J),J=MML,MMH,LVV)
1800        FORMAT (' ',I3,3X,1P4E15.6) 
C
 1801 CONTINUE
C
      RETURN
      END 
C
C
      SUBROUTINE SETUP1 (T,M,INV,Y,ISYM,TPL,YPL)
C
C***********************************************************************
C
C  THIS SUBROUTINE SETS UP ARRAYS CONTAINING THE VALUES TO BE
C  PLOTTED ON A GRAPH PLOTTING ONLY ONE SET OF POINTS
C
      LOGICAL IP,ITP
      COMMON /BLOKIP/ IP
      COMMON /BLKITP/ ITP
      DIMENSION T(M,INV),Y(M),ISYM(M),TPL(M),YPL(M)
      DATA IPR6 /6/ ,IPR10 /10/
C
      DO 10 K=1,M
         TPL(K)=T(K,INV)
         YPL(K)=Y(K)
         ISYM(K)=3
10    CONTINUE
C
C
      IF (ITP) CALL SPLTH (YPL,TPL,M,ISYM,IPR6)
C
      IF (IP) CALL SPLT (YPL,TPL,M,ISYM,IPR10)
C
      RETURN
C
C
C
      END 
C***********************************************************************
C
      SUBROUTINE SETUP2 (T,M,INV,Y,ETA,ISYM,TPL,ETAPL)
C
C***********************************************************************
C
C  THIS SUBROUTINE SETS UP ARRAYS CONTAINING THE VALUES TO BE
C  PLOTTED ON A GRAPH PLOTTING TWO SETS OF POINTS 
C
      LOGICAL IP,ITP
      COMMON /BLOKIP/ IP
      COMMON /BLKITP/ ITP
      DIMENSION T(M,INV),Y(M),ETA(M),ISYM(2*M),TPL(2*M),ETAPL(2*M)
      DATA IPR6 /6/ ,IPR10 /10/
C
      DO 20 K=1,M
         TPL(K)=T(K,INV)
         ETAPL(K)=Y(K)
         ISYM(K)=3
20    CONTINUE
C
      MP1=M+1
      ITWOM=2*M
      DO 30 K=MP1,ITWOM
         KL=K-M
         TPL(K)=T(KL,INV)
         ETAPL(K)=ETA(KL)
         ISYM(K)=1
30    CONTINUE
C
      IF (ITP) CALL SPLTH (ETAPL,TPL,ITWOM,ISYM,IPR6)
C
      IF (IP) CALL SPLT (ETAPL,TPL,ITWOM,ISYM,IPR10)
C
      RETURN
C
C
C
      END 
C***********************************************************************
C
      SUBROUTINE SPLT (YARAY,XARAY,M,ISYM,IPR)
C
C***********************************************************************
C
C  THIS SUBROUTINE SETS UP THE CALL TO PLTMAIN FOR A WHOLE PAGE PLOT
C  OF X VS Y IN WHICH THE COMPUTER SETS UP THE GRAPH LIMITS AND THE
C  USER SPECIFIES THE N VECTOR ISYM IN WHICH THE PLOTTING SYMBOLS ARE 
C  PROVIDED
C
C  WRITTEN BY --
C              JANET R. DONALDSON
C              STATISTICAL ENGINEERING LABORATORY/BOULDER
C              NATIONAL BUREAU OF STANDARDS
C
C  LATEST REVISION  -  MAY 30, 1978  (JRD)
C
C
C  INPUT ARGUMENTS - (IN ORDER OF USAGE)
C
C  XARAY THE N VECTOR OF OBSERVATIONS FOR THE X (HORIZONTAL)
C        COORDINATES
C  YARAY THE N VECTOR OF OBSERVATIONS FOR THE Y (VERTICAL)
C        COORDINATES
C  M     THE INTEGER NUMBER OF OBSERVATIONS TO BE PLOTTED
C  ISYM  THE N VECTOR OF INTEGERS DETERMINING THE PLOTTING SYMBOLS TO 
C        USED FOR EACH POINT INDIVIDUALLY, WHERE
C            1 = +    2 = .    3 = *    4 = -
C            5 = A    6 = B    7 = C    8 = D
C            9 = E   10 = F   11 = G   12 = H
C           13 = I   14 = J   15 = K   16 = L
C           17 = M   18 = N   19 = O   20 = P
C           21 = Q   22 = R   23 = S   24 = T
C           25 = U   26 = V   27 = W   28 = Y
C           29 = Z   30 (AND ABOVE) = $ 
C        THE NUMBERS 1 TO 9 INDICATE MULTIPLE POINTS ON A GIVEN
C        PLOT LOCATION, WHERE THE NUMBER INDICATES HOW MANY POINTS
C        ARE REPRESENTED
C        NOTE THAT  X  IS NOT USED AS A PLOTTING SYMBOL EXCEPT TO
C        INDICATE THAT MORE THAN 9 POINTS FELL ON THE SAME PLOT
C        LOCATION
C
C
C  ADDITIONAL VARIABLES USED
C
C  IPCODE   THE INTEGER CODE, USED IN ERROR CHECKING, WHICH DETERMINES
C           WHICH PLOT SUBROUTINE HAS BEEN CALLED 
C  NSIZE    THE INTEGER CODE FOR THE PLOT SIZE, WHERE
C           0 INDICATES A PLOT 100 COL BY 50 ROWS 
C
C
C  ADDITIONAL SUBROUTINES REQUIRED
C     PLTERR,PLTMN
C
C  MACHINE DEPENDENT VARIABLES - NONE
C
      REAL XARAY(M),YARAY(M)
      INTEGER ISYM(M)
C
      COMMON /ERRCHK/ IERROR
C
      NSIZE=0
      IERROR=0
      IPCODE=1
      CALL PLTERR (M,IPCODE,IPR)
      IF (IERROR.NE.0) RETURN 
      CALL PLTMN (YARAY,XARAY,M,ISYM,NSIZE,IPR)
      RETURN
C
C
C
      END 
C***********************************************************************
C
      SUBROUTINE SPLTH (YARAY,XARAY,M,ISYM,IPR)
C
C***********************************************************************
C
C  THIS SUBROUTINE SETS UP THE CALL TO PLTMAIN FOR A HALF PAGE PLOT
C  OF X VS Y IN WHICH THE COMPUTER SETS UP THE GRAPH LIMITS AND THE
C  USER SPECIFIES THE N VECTOR ISYM IN WHICH THE PLOTTING SYMBOLS ARE 
C  PROVIDED
C
C  WRITTEN BY --
C              JANET R. DONALDSON
C              STATISTICAL ENGINEERING LABORATORY/BOULDER
C              NATIONAL BUREAU OF STANDARDS
C
C  LATEST REVISION  -  MAY 30, 1978  (JRD)
C
C
C  INPUT ARGUMENTS - (IN ORDER OF USAGE)
C
C  XARAY THE N VECTOR OF OBSERVATIONS FOR THE X (HORIZONTAL)
C        COORDINATES
C  YARAY THE N VECTOR OF OBSERVATIONS FOR THE Y (VERTICAL)
C        COORDINATES
C  M     THE INTEGER NUMBER OF OBSERVATIONS TO BE PLOTTED
C  ISYM  THE N VECTOR OF INTEGERS DETERMINING THE PLOTTING SYMBOLS TO 
C        USED FOR EACH POINT INDIVIDUALLY, WHERE
C            1 = +    2 = .    3 = *    4 = -
C            5 = A    6 = B    7 = C    8 = D
C            9 = E   10 = F   11 = G   12 = H
C           13 = I   14 = J   15 = K   16 = L
C           17 = M   18 = N   19 = O   20 = P
C           21 = Q   22 = R   23 = S   24 = T
C           25 = U   26 = V   27 = W   28 = Y
C           29 = Z   30 (AND ABOVE) = $ 
C        THE NUMBERS 1 TO 9 INDICATE MULTIPLE POINTS ON A GIVEN
C        PLOT LOCATION, WHERE THE NUMBER INDICATES HOW MANY POINTS
C        ARE REPRESENTED
C        NOTE THAT  X  IS NOT USED AS A PLOTTING SYMBOL EXCEPT TO
C        INDICATE THAT MORE THAN 9 POINTS FELL ON THE SAME PLOT
C        LOCATION
C
C
C  ADDITIONAL VARIABLES USED
C
C  IPCODE   THE INTEGER CODE, USED IN ERROR CHECKING, WHICH DETERMINES
C           WHICH PLOT SUBROUTINE HAS BEEN CALLED 
C  NSIZE    THE INTEGER CODE FOR THE PLOT SIZE, WHERE
C           1 INDICATES A PLOT  50 COL BY 50 ROWS 
C
C
C  ADDITIONAL SUBROUTINES REQUIRED
C     PLTERR,PLTMN
C
C  MACHINE DEPENDENT VARIABLES - NONE
C
      REAL XARAY(M),YARAY(M)
      INTEGER ISYM(M)
C
      COMMON /ERRCHK/ IERROR
C
      NSIZE=1
      IERROR=0
      IPCODE=2
      CALL PLTERR (M,IPCODE,IPR)
      IF (IERROR.NE.0) RETURN 
      CALL PLTMN (YARAY,XARAY,M,ISYM,NSIZE,IPR)
      RETURN
C
C
C
      END 
C***********************************************************************
C
      SUBROUTINE PLTMN (YARAY,XARAY,M,ISYM,NSIZE,IPR)
C
C***********************************************************************
C
C       THIS IS THE MAIN SUBROUTINE FOR ONE PAGE PLOTS
C
C  THE FOLLOWING PLOTS ARE SET UP FROM THIS ROUTINE
C
C  ONE PAGE PLOTS WITH ONLY THE PLOTTING SYMBOL + USED
C     - WHOLE PAGE PLOTS, ROUTINE COMPUTES GRAPH LIMITS (PLOT)
C     - HALF PAGE PLOT, ROUTINE COMPUTES GRAPH LIMITS (PLOTH)
C
C  ONE PAGE PLOT OF MULTI COLUMN MATRIX Y VERSUS VECTOR X
C     - WHOLE PAGE PLOTS, ROUTINE COMPUTES GRAPH LIMITS (MPLOT)
C     - HALF PAGE PLOT, ROUTINE COMPUTES GRAPH LIMITS (MPLOTH)
C
C  WRITTEN BY  --
C              JANET R. DONALDSON
C              STATISTICAL ENGINEERING LABORATORY/BOULDER
C              NATIONAL BUREAU OF STANDARDS
C
C  EDITED BY --
C            CLAIRE WOLFE
C            SCIENTIFIC COMPUTING DIVISION
C            NATIONAL BUREAU OF STANDARDS
C
C  LATEST REVISION  -  OCTOBER 5,1982  (CMW)
C
C
C  INPUT ARGUMENTS - (IN ORDER OF USAGE)
C
C  XARAY THE N VECTOR OF OBSERVATIONS FOR THE X (HORIZONTAL)
C        COORDINATES
C  YARAY THE N BY M MATRIX OF OBSERVATIONS FOR THE Y (VERTICAL)
C        COORDINATES
C        THE ITH COLUMN IS PLOTTED USING THE FOLLOWING SYMBOL
C            1 = A    2 = B    3 = C    4 = D
C            5 = E    6 = F    7 = G    8 = H
C            9 = I   10 = J   11 = K   12 = L
C           13 = M   14 = N   15 = O   16 = P
C           17 = M   18 = N   19 = O   20 = P
C           21 = U   22 = V   23 = W   24 = Y
C           25 = Z   26 (AND ABOVE) = $ 
C        THE NUMBERS 1 TO 9 INDICATE MULTIPLE POINTS ON A GIVEN
C        PLOT LOCATION, WHERE THE NUMBER INDICATES HOW MANY POINTS
C        ARE REPRESENTED
C        NOTE THAT  X  IS NOT USED AS A PLOTTING SYMBOL EXCEPT TO
C        INDICATE THAT MORE THAN 9 POINTS FELL ON THE SAME PLOT
C        LOCATION
C  N     THE INTEGER NUMBER OF OBSERVATIONS TO BE PLOTTED (IN EACH
C        COLUMN)
C  ISYM  THE N VECTOR OF INTEGERS DETERMINING THE PLOTTING SYMBOLS TO 
C        USED FOR THE S SERIES OF PLOTS, WHERE
C            1 = +    2 = .    3 = *    4 = -
C            5 = A    6 = B    7 = C    8 = D
C            9 = E   10 = F   11 = G   12 = H
C           13 = I   14 = J   15 = K   16 = L
C           17 = M   18 = N   19 = O   20 = P
C           21 = Q   22 = R   23 = S   24 = T
C           25 = U   26 = V   27 = W   28 = Y
C           29 = Z   30 (AND ABOVE) = $ 
C        THE NUMBERS 1 TO 9 INDICATE MULTIPLE POINTS ON A GIVEN
C        PLOT LOCATION, WHERE THE NUMBER INDICATES HOW MANY POINTS
C        ARE REPRESENTED
C        NOTE THAT  X  IS NOT USED AS A PLOTTING SYMBOL EXCEPT TO
C        INDICATE THAT MORE THAN 9 POINTS FELL ON THE SAME PLOT
C        LOCATION
C  NSIZE THE INTEGER CODE FOR THE PLOT SIZE, WHERE
C        0 INDICATES A PLOT 100 COL BY 50 ROWS
C        1 INDICATES A PLOT 50 COL BY  50 ROWS
C
C
C  ADDITIONAL VARIABLES USED - (IN ALPHABETICAL ORDER)
C
C  ALINE    THE VECTOR OF THE CURRENT PLOT LINE
C  ALPHAI   THE PLOT AXIS SYMBOL I
C  AXISCH   THE Y A AXIS SYMBOL TO BE USED FOR THE CURRENT LINE,
C           EITHER I OR -
C  ALPHAX   THE PLOTTING SYMBOL X DESIGNATING MORE THAN 9 POINTS FELL ON
C           A SINGLE PLOTTING LOCATION
C  BLANK    THE PLOTTING SYMBOL BLANK
C  DELX     THE RANGE OF THE X AXIS
C  DELY     THE RANGE OF THE Y AXIS
C  HYPHEN   THE PLOT AXIS SYMBOL -
C  IC       THE COUNT OF THE NUMBER OF VALUES FALLING OUTSIDE OF THE
C           GRAPH BOUNDS
C  ICOL     THE COLUMN LOCATION FOR THE PLOT LINE 
C  IPCODE   THE INTEGER CODE, USED IN ERROR CHECKING, WHICH DETERMINES
C           WHICH PLOT SUBROUTINE HAS BEEN CALLED 
C  IPR      THE LOGICAL UNIT NUMBER OF THE PRINTER
C  ITEST    THE INDICATOR VARIABLE FOR WHETHER THE X AXIS LABELS ARE
C           PRINTED IN E OR F FORMAT
C  ISS      THE INTEGER VALUE USED TO DETERMINE THE SPECIFIC SYMBOL TO
C           BE PLOTTED
C  KSS      AN INTEGER VECTOR USED IN DETERMINING THE PLOT SYMBOL
C           NEEDED
C  NN       THE NUMBER OF Y LABELS TO BE LISTED ON THE LEFT AXIS,
C           (DEPENDENT ON THE GRAPH SIZE)
C  NUMCOL   THE INTEGER VALUE OF THE NUMBER OF COLUMNS IN THE GRAPH
C  NUMLAB   THE INTEGER NUMBER OF X LABELS TO BE LISTED AT THE BOTTOM 
C           OF THE GRAPH (DEPENDENT ON GRAPH SIZE)
C  NUMROW   THE INTEGER NUMBER OF ROWS IN THE GRAPH
C  SYM      THE VECTOR OF PLOT SYMBOL ASSIGNMENTS (SEE ISS ABOVE)
C  SYM1     THE VECTOR OF INTEGER VALUES USED TO INDICATE MULTIPLE
C           POINTS ON THE SAME PLOT LOCATION
C  TEMP     THE ARRAY OF VALUES TO BE PRINTED WHICH FALL OUTSIDE THE
C           GRAPH LIMITS
C  XLAB     THE VECTOR OF X AXIS LABELS 
C  XN       THE MINIMUM X VALUE TO BE PLOTTED, COMPUTED FROM DATA OR
C           ASSIGNED BY XMIN
C  XWIDTH   THE VALUE OF AN INDIVIDUAL X AXIS GRAPH INTERVAL
C  XX       THE MAXIMUM X VALUE TO BE PLOTTED, COMPUTED FROM DATA OR
C           ASSIGNED BY XMAX
C  YLAB     THE VALUE OF THE Y AXIS LABEL TO BE PRINTED
C  YLOWER   THE LOWER BOUND FOR Y VALUES TO BE PLOTTED ON THE CURRENT 
C           LINE
C  YN       THE MINIMUM Y VALUE TO BE PLOTTED, COMPUTED FROM DATA OR
C           ASSIGNED BY YMIN
C  YUPPER   THE UPPER BOUND FOR Y VALUES TO BE PLOTTED ON THE CURRENT 
C           LINE
C  YWIDTH   THE VALUE OF AN INDIVIDUAL Y AXIS GRAPH INTERVAL
C  YX       THE MAXIMUM Y VALUE TO BE PLOTTED, COMPUTED FROM DATA OR
C           ASSIGNED BY YMAX
C
C
C  MACHINE DEPENDENT VARIABLES DEFINED  -  NONE
C
      REAL XARAY(M),YARAY(M),XLAB(15),YLAB(20)
      INTEGER ISYM(M),KSS(101)
      CHARACTER*1 SYM(30),SYM1(9),ALINE(103),BLANK,HYPHEN,
     2        ALPHAI,ALPHAX,AXISCH
      DATA BLANK /' '/ ,HYPHEN /'-'/ ,ALPHAI /'I'/ ,ALPHAX /'X'/
      DATA SYM(1) /'+'/ ,SYM(2) /'.'/ ,SYM(3) /'*'/ ,SYM(4) /'-'/ ,
     2     SYM(5) /'A'/ ,SYM(6) /'B'/ ,SYM(7) /'C'/ ,SYM(8) /'D'/ ,
     3     SYM(9) /'E'/ ,SYM(10) /'F'/ ,SYM(11) /'G'/ ,SYM(12) /'H'/ ,
     4     SYM(13) /'I'/ ,SYM(14) /'J'/ ,SYM(15) /'K'/ ,SYM(16) /'L'/ ,
     5     SYM(17) /'M'/ ,SYM(18) /'N'/ ,SYM(19) /'O'/ ,SYM(20) /'P'/ ,
     6     SYM(21) /'Q'/ ,SYM(22) /'R'/ ,SYM(23) /'S'/ ,SYM(24) /'T'/ ,
     7     SYM(25) /'U'/ ,SYM(26) /'V'/ ,SYM(27) /'W'/ ,SYM(28) /'Y'/ ,
     8     SYM(29) /'Z'/ ,SYM(30) /'$'/ 
      DATA SYM1(1) /'1'/ ,SYM1(2) /'2'/ ,SYM1(3) /'3'/ ,SYM1(4) /'4'/ ,
     2     SYM1(5) /'5'/ ,SYM1(6) /'6'/ ,SYM1(7) /'7'/ ,SYM1(8) /'8'/ ,
     3     SYM1(9) /'9'/
C
C  DEFINE THE NUMBER OF ROWS AND COLUMNS WITHIN THE PLOT
C
      NUMCOL=101
      IF (NSIZE.NE.0) NUMCOL=51
      NUMCP2=NUMCOL+2
      NUMROW=51
      NN=(NUMROW-1)/5
C
C  ROUTINE DETERMINES LIMITS FOR BOTH AXES
C
      YN=YARAY(1)
      YX=YARAY(1)
      XN=XARAY(1)
      XX=XARAY(1)
      DO 70 I=1,M
         IF (YARAY(I).GE.YN) GO TO 40
         YN=YARAY(I)
40       IF (YARAY(I).LE.YX) GO TO 50
         YX=YARAY(I)
50       IF (XARAY(I).GE.XN) GO TO 60
         XN=XARAY(I)
60       IF (XARAY(I).LE.XX) GO TO 70
         XX=XARAY(I)
70    CONTINUE
C
C  CHECK FOR ERRORS IN THE PLOT LIMITS
C
      IF (XN.LT.XX.AND.YN.LT.YX) GO TO 80
      IF (XN.EQ.XX) WRITE (IPR,240) XN
      IF (YN.EQ.YX) WRITE (IPR,330) YN
      RETURN
C
C  DETERMINE THE VALUE OF A DIVISION FOR X AND Y
C
80    DELY=YX-YN
      DELX=XX-XN
      YWIDTH=DELY/FLOAT(NUMROW-1)
      XWIDTH=DELX/FLOAT(NUMCOL-1)
C
C  TEST FOR FORMAT FOR Y AXIS LABELS
C
      ITEST=0
      IF ((YN.GE.1.0E6.OR.YN.LE.(-1.0E5)).OR.(ABS(YN).GT.0.0.AND.
     2    ABS(YN).LT.0.001)) ITEST=1
      IF ((YX.GE.1.0E6.OR.YX.LE.(-1.0E5)).OR.(ABS(YX).GT.0.0.AND.
     2    ABS(YX).LT.0.001)) ITEST=1
C
C  BLANK OUT THE PLOT PRINT LINE ALINE
C
      DO 90 ICOL=1,103
         ALINE(ICOL)=BLANK
90    CONTINUE
C
C  WRITE OUT THE TOP HORIZONTAL AXIS OF THE PLOT, AND
C  THE VALUE OF A DIVISION FOR THE X AND Y AXIS LISTED
C
      WRITE (IPR,260) YWIDTH,XWIDTH
      DO 100 ICOL=1,NUMCOL
         ALINE(ICOL)=HYPHEN
100   CONTINUE
      ALINE(1)=ALPHAI
      DO 110 ICOL=1,NUMCOL,10 
         ALINE(ICOL)=ALPHAI
110   CONTINUE
      WRITE (IPR,250) HYPHEN,(ALINE(I),I=1,NUMCOL),HYPHEN
C
C  DETERMINE AND WRITE OUT THE PLOT POSITIONS ONE LINE AT A TIME.
C
      IC=0
      LN=-1
      IK1=1
      DO 190 IROW=1,NUMROW
         LN=LN+1
         IF (LN.GT.NN) LN=1
         DO 120 ICOL=1,NUMCOL 
            ALINE(ICOL)='0'
            KSS(ICOL)=0
120      CONTINUE
         AIROW=IROW 
         YUPPER=YX+(1.5-AIROW)*YWIDTH
         YLOWER=YX+(0.5-AIROW)*YWIDTH
         DO 130 I=1,M
            IF (YLOWER.GT.YARAY(I).OR.YARAY(I).GE.YUPPER) GO TO 130
            IF (YN.GT.YARAY(I).OR.YARAY(I).GT.YX) GO TO 130 
            IF (XARAY(I).LT.XN.OR.XARAY(I).GT.XX) GO TO 130 
            ICOL=((XARAY(I)-XN)/XWIDTH)+1.5
            ALINE(ICOL)=CHAR(ICHAR(ALINE(ICOL))+1)
C       SYMBOL PLOTS (SPLT)
            KSS(ICOL)=ISYM(I) 
C       MULTIPLE PLOTS (MPLT) 
            IF (KSS(ICOL).GT.30) KSS(ICOL)=30
            IC=IC+1 
130      CONTINUE
         ISTOP = 1
         DO 170 ICOL=1,NUMCOL 
            ISS=ICHAR(ALINE(ICOL))-ICHAR('0')
            IF (ISS.NE.0) GO TO 140
            ALINE(ICOL)=BLANK 
            GO TO 170
140         IF (ISS.GE.2) GO TO 150
            K=KSS(ICOL)
            ALINE(ICOL)=SYM(K)
            ISTOP = ICOL
            GO TO 170
150         IF (ISS.GT.9) GO TO 160
            ALINE(ICOL)=SYM1(ISS)
            ISTOP = ICOL
            GO TO 170
160         ALINE(ICOL)=ALPHAX
            ISTOP = ICOL
170      CONTINUE
         ALINE(NUMCP2)=HYPHEN 
         AXISCH=HYPHEN
         IF (IROW.NE.1.AND.IROW.NE.NUMROW.AND.LN.NE.NN) GO TO 180
         IK1=1
         YLAB(IK1)=YX+(1.0-AIROW)*YWIDTH
         IF (ITEST.EQ.0) WRITE (IPR,270) YLAB(IK1),AXISCH,(ALINE(ICOL),
     2                                   ICOL=1,NUMCP2)
         IF (ITEST.EQ.1) WRITE (IPR,280) YLAB(IK1),AXISCH,(ALINE(ICOL),
     2                                   ICOL=1,NUMCP2)
         IK1=IK1+1
         GO TO 190
180      ALINE(NUMCP2)=ALPHAI 
         AXISCH=ALPHAI
         IF (NUMCP2.LE.53) WRITE (IPR,290) AXISCH,(ALINE(ICOL),
     2                                     ICOL=1,ISTOP)
         IF (NUMCP2.GT.53) WRITE (IPR,290) AXISCH,(ALINE(ICOL),
     2                                     ICOL= 1,NUMCP2)
190   CONTINUE
C
C  SKIP 1 LINE FOR A BOTTOM MARGIN WITHIN THE PLOT, WRITE OUT THE
C  BOTTOM HORIZONTAL AXIS, AND WRITE OUT THE X AXIS LABELS. 
C
      DO 200 ICOL=1,103
         ALINE(ICOL)=BLANK
200   CONTINUE
      ALINE(NUMCOL+2)=ALPHAI
      IF (NUMCP2.LE.53) WRITE (IPR,290) ALPHAI
      IF (NUMCP2.GT.53) WRITE (IPR,290) ALPHAI,(ALINE(ICOL),
     2                                  ICOL=1,NUMCP2)
      DO 210 ICOL=1,NUMCOL
         ALINE(ICOL)=HYPHEN
210   CONTINUE
      ALINE(1)=ALPHAI
      DO 220 ICOL=11,NUMCOL,10
         ALINE(ICOL)=ALPHAI
220   CONTINUE
      ALINE(NUMCOL+1)=HYPHEN
      ALINE(NUMCOL+2)=BLANK
      WRITE (IPR,250) HYPHEN,(ALINE(ICOL),ICOL=1,NUMCP2)
      NUMLAB=NUMCOL/10+1
      ANUMLM=NUMLAB-1
      ITEST=0
      DO 230 I=1,NUMLAB
         AIM1=I-1
         XLAB(I)=XN+(AIM1/ANUMLM)*DELX
         IF (ABS(XLAB(I)).GT.0.0.AND.ABS(XLAB(I)).LT.0.01) ITEST=1
         IF (XLAB(I).GE.1.0E4.OR.XLAB(I).LE.(-1.0E3)) ITEST=1
230   CONTINUE
      IF (ITEST.EQ.0) WRITE (IPR,300) (XLAB(I),I=1,NUMLAB)
      IF (ITEST.EQ.1) WRITE (IPR,320) (XLAB(I),I=1,NUMLAB)
C
C  CHECK FOR POINTS OUTSIDE OF GRAPH LIMITS
C  THE TOTAL NUMBER OF POINTS TO BE PLOTTED IS M
C
      IC=M-IC
      IF (IC.EQ.0) RETURN
      WRITE (IPR,310) IC
      RETURN
C
C  FORMAT STATEMENTS
C
240   FORMAT(/52H THE VALUES IN THE X ARRAY ARE ALL IDENTICALLY EQUAL/
     2        4H TO ,E13.7)
250   FORMAT (1H ,16X,105A1)
260   FORMAT (1H ,15X,'   1 Y DIV=',G12.5,' AND 1 X DIV=',G12.5)
270   FORMAT (4X,F11.4,1X,A1,1X,103A1)
280   FORMAT (1H ,E14.7,1X,A1,1X,103A1) 
290   FORMAT (1H ,15X,A1,1X,103A1)
300   FORMAT (13X,11(1X,F9.4))
310   FORMAT (16X,6H**NOTE,I4,
     2        43H VALUES FELL OUTSIDE THE SPECIFIED LIMITS**)
320   FORMAT (14X,11E10.3)
330   FORMAT(/52H THE VALUES IN THE Y ARRAY ARE ALL IDENTICALLY EQUAL/
     2        4H TO ,E13.6)
C
C
C
      END 
C***********************************************************************
C
      SUBROUTINE PLTERR (M,IPCODE,IPR)
C
C***********************************************************************
C
C  LATEST REVISION  -  OCTOBER 5, 1982  (CMW)
C
      COMMON /ERRCHK/ IERROR
C
      IERROR=0
      IF (M.GE.1) RETURN
      IERROR=1
      WRITE (IPR,420)
      WRITE (IPR,370)
      IF (M.GE.1) GO TO 340
      WRITE (IPR,380) M
340   WRITE (IPR,410)
      GO TO (350,360), IPCODE 
350   WRITE (IPR,390)
      RETURN
360   WRITE (IPR,400)
      RETURN
C
C  FORMAT STATEMENTS
C
370   FORMAT (51H AN ERROR HAS BEEN DETECTED IN THE PLOTTING ROUTINE, 
     2        15H ARGUMENT LIST.)
380   FORMAT (/5X,' THE NUMBER OF OBSERVATIONS TO BE PLOTTED, M, MUST',
     2        /5X,' BE GREATER THAN OR EQUAL TO 1.  THE INPUT VALUE ',
     3        'OF M IS ',I6,1H.)
390   FORMAT (10X,' CALL SPLT(YARAY, XARAY, M, ISYM, IPR)') 
400   FORMAT (10X,' CALL SPLTH(YARAY, XARAY, M, ISYM, IPR)')
410   FORMAT (/5X,34H THE FORM OF THE CALL STATEMENT IS)
420   FORMAT (1H1)
C
C
C
      END 
C***********************************************************************
C
      SUBROUTINE ORPLOT (T,M,IOR)
C
C***********************************************************************
C
C  ARRANGE THE INDEPENDENT VARIABLES TO BE USED IN ASCENDING
C  ORDER.  THE INDEXING IS STORED IN ARRAY IOR.
C
      PARAMETER ( MMAX=300 )
      DIMENSION IOR(M)
      REAL T(M),SAV(MMAX),HN,H
C
      DO 430 I=1,M
         SAV(I)=T(I)
430   CONTINUE
      H=SAV(1)
      IX=1
      DO 440 I=2,M
         HN=AMAX1(SAV(I),H)
         IF (HN.NE.H) IX=I
         H=HN
440   CONTINUE
      IOR(M)=IX
      MAX=IX
      MM1=M-1
      DO 460 J=1,MM1
         H=SAV(1)
         IX=1
         DO 450 I=2,M
            HN=AMIN1(SAV(I),H)
            IF (HN.NE.H) IX=I 
            H=HN
450      CONTINUE
         IOR(J)=IX
         SAV(IX)=SAV(MAX)
460   CONTINUE
      RETURN
C
C
C
      END 
C***********************************************************************
C
      SUBROUTINE VPLT (X1,M,IPR,IMEAN)
C
C***********************************************************************
C
C  THIS SUBROUTINE SETS UP THE CALL TO VPMAIN TO PLOT A SINGLE
C  VECTOR VERSUS THE INPUT ORDER ON A VERTICAL SCALE
C
C
C  WRITTEN BY - JANET R. DONALDSON
C               STATISTICAL ENGINEERING LABORATORY/BOULDER
C               NATIONAL BUREAU OF STANDARDS
C
C  LATEST REVISION  -  MAY 30, 1978  (JRD)
C
C
C  INPUT ARGUMENTS - (IN ORDER OF USAGE)
C
C  X1      THE N VECTOR CONTAINING THE 1ST OBSERVED SERIES
C  M       THE INTEGER NUMBER OF OBSERVATIONS IN EACH SERIES
C
C
C  ADDITIONAL VARIABLES USED - (IN ALPHABETICAL ORDER)
C
C  IMEAN   THE INTEGER INDICATOR AS TO WHETHER THE MEAN OR ZERO VALUE 
C          IS TO BE PLOTTED.
C          IF IMEAN.NE.0 GRAPH WILL PLOT MEAN VALUE
C  IPCODE  THE INTEGER CODE, USED IN ERROR CHECKING, WHICH DETERMINES 
C          WHICH PLOT SUBROUTINE HAS BEEN CALLED
C  LMSPEC  THE INTEGER VALUE USED TO DESIGNATE WHICH SET OF SPECIFIED 
C          LIMITS ARE TO BE USED WHEN NSER.EQ.2
C          IF NSER.EQ.1, THIS IS A DUMMY VARIABLE 
C  NSER    THE INTEGER WHICH INDICATES THE NUMBER OF SERIES (VECTORS) 
C          TO BE PLOTTED
C          IF NSER.EQ.1, ONE SERIES IS PLOTTED
C  X2      IN THIS CASE, A DUMMY VARIABLE
C
C
C  ADDITIONAL SUBROUTINES REQUIRED
C
C     VPLTER,VPMAIN(VPLTER)
C
C
C  MACHINE DEPENDENT VARIABLES - NONE
C
      REAL X1(M),X2(1)
      COMMON /ERRCHK/ IERROR
C
      NSER=1
      LMSPEC=1
      IERROR=0
      IPCODE=1
      CALL VPLTER (M,NSER,LMSPEC,IPCODE,IPR)
      IF (IERROR.NE.0) RETURN 
      IF (IPR.EQ.6) CALL VMAINH (X1,X2,M,NSER,LMSPEC,IMEAN,IPR)
      IF (IPR.EQ.10) CALL VPMAIN (X1,X2,M,NSER,LMSPEC,IMEAN,IPR)
      RETURN
C
C
C
      END 
C***********************************************************************
C
      SUBROUTINE VPLT2 (X1,X2,M,IPR,IMEAN)
C
C***********************************************************************
C
C  THIS SUBROUTINE SETS UP A CALL TO VPMAIN TO PLOT TWO VECTORS
C  VERSUS THE INPUT ORDER ON A VERTICAL SCALE, IN WHICH THE X AXIS
C  LIMITS ARE COMPUTED FROM THE MINIMUM AND MAXIMUM VALUES OF THE TWO 
C  VECTORS
C
C  WRITTEN BY - JANET R. DONALDSON
C               STATISTICAL ENGINEERING LABORATORY/BOULDER
C               NATIONAL BUREAU OF STANDARDS
C
C  LATEST REVISION  -  MAY 31, 1978  (JRD)
C
C
C  INPUT ARGUMENTS - (IN ORDER OF USAGE)
C
C  X1    THE N VECTOR CONTAINING THE 1ST OBSERVED SERIES
C  X2    THE N VECTOR CONTAINING THE 2ND OBSERVED SERIES
C  M     THE INTEGER NUMBER OF OBSERVATIONS IN EACH SERIES
C
C
C  ADDITIONAL VARIABLES USED - (IN ALPHABETICAL ORDER)
C
C  IMEAN    THE INTEGER INDICATOR AS TO WHETHER THE MEAN OR ZERO VALUE
C           IS TO BE PLOTTED. 
C           IF IMEAN.NE.0 GRAPH WILL PLOT MEAN VALUE
C  IPCODE   THE INTEGER CODE, USED IN ERROR CHECKING, WHICH DETERMINES
C           WHICH PLOT SUBROUTINE HAS BEEN CALLED 
C  LMSPEC   THE INTEGER VALUE USED TO DESIGNATE WHICH SET OF SPECIFIED
C           LIMITS ARE TO BE USED WHEN NSER.EQ.2
C           FOR THIS CASE, LMSPEC SHOULD BE SET EQUAL TO 0, 
C           DESIGNATING THAT THE LIMITS ARE TO BE SPECIFIED FROM THE
C           MINIMUM AND MAXIMUM OF BOTH SERIES
C  NSER     THE INTEGER WHICH INDICATES THE NUMBER OF SERIES (VECTORS)
C           TO BE PLOTTED
C           IF NSER.EQ.2, TWO SERIES ARE PLOTTED
C
C  ADDITIONAL SUBROUTINES REQUIRED  -
C     VPLTER,VPMAIN 
C
C  MACHINE DEPENDENT VARIABLES DEFINED  -  NONE
C
      REAL X1(M),X2(M)
      COMMON /ERRCHK/ IERROR
C
      LMSPEC=0
      NSER=2
      IERROR=0
      IPCODE=2
      CALL VPLTER (M,NSER,LMSPEC,IPCODE,IPR)
      IF (IERROR.NE.0) RETURN 
      IF (IPR.EQ.6) CALL VMAINH (X1,X2,M,NSER,LMSPEC,IMEAN,IPR)
      IF (IPR.EQ.10) CALL VPMAIN (X1,X2,M,NSER,LMSPEC,IMEAN,IPR)
      RETURN
C
C
C
      END 
C***********************************************************************
C
      SUBROUTINE VPLTER (M,NSER,LMSPEC,IPCODE,IPR)
C
C***********************************************************************
C
C  LATEST REVISION  -  JUNE 30, 1978  (JRD)
C
      COMMON /ERRCHK/ IERROR
C
      IERROR=0
      IF (M.GE.1.AND.(NSER.EQ.1.OR.NSER.EQ.2).AND.(LMSPEC.EQ.0.OR.
     2    LMSPEC.EQ.1)) RETURN
      IERROR=1
      WRITE (IPR,590)
      WRITE (IPR,520)
      IF (M.GE.1) GO TO 470
      WRITE (IPR,530) M
470   IF (NSER.EQ.1.OR.NSER.EQ.2) GO TO 480
      WRITE (IPR,540) NSER
480   IF (LMSPEC.EQ.0.OR.LMSPEC.EQ.1) GO TO 490
      WRITE (IPR,550) LMSPEC
490   WRITE (IPR,580)
      GO TO (500,510), IPCODE 
500   WRITE (IPR,560)
      RETURN
510   WRITE (IPR,570)
      RETURN
C
C  FORMAT STATEMENTS
C
520   FORMAT (51H AN ERROR HAS BEEN DETECTED IN THE PLOTTING ROUTINE, 
     2        15H ARGUMENT LIST.)
530   FORMAT (5X/' ','THE NUMBER OF OBSERVATIONS TO BE PLOTTED, N,'
     2        /5X,'MUST BE GREATER THAN OR EQUAL TO 1.  THE INPUT ',
     3        'VALUE OF N IS ',I6,1H.)
540   FORMAT (5X/' ','THE NUMBER OF SERIES TO BE PLOTTED, NSER, MUST' 
     2        /5X,' BE EQUAL TO 1 OR 2. THE INPUT VALUE OF NSER IS ', 
     3        I6,1H.)
550   FORMAT (5X/' ','THE INDICATOR VARIABLE USED TO DETERMINE HOW '/5X,
     2        ' THE SCALES ARE TO BE COMPUTED, LMSPEC, MUST BE '/5X,
     3        ' EQUAL TO A 0, 1 OR 2.  THE INPUT VALUE OF LMSPEC IS ',
     4        I6,1H.)
560   FORMAT (10X,' CALL VPLT(X1, M, IPR, IMEAN)')
570   FORMAT (10X,' CALL VPLT2(X1, X2, M, IPR, IMEAN)')
580   FORMAT (5X/34H THE FORM OF THE CALL STATEMENT IS)
590   FORMAT (1H1)
C
C
C
      END 
C***********************************************************************
C
      SUBROUTINE VPMAIN (X1,X2,M,NSER,LMSPEC,IMEAN,IPR)
C
C***********************************************************************
C
C  THIS IS THE MAIN SUBROUTINE FOR PRODUCING VERTICAL MULTIPAGE
C  PLOTS
C  THE FOLLOWING PLOTS ARE PRODUCED BY THIS ROUTINE -
C
C     VERTICAL PLOT OF ONE VECTOR, COMPUTED LIMITS (VPLOT)
C     VERTICAL PLOT OF ONE VECTOR, COMPUTED LIMITS, CONNECTED POINTS
C        (VPLOTB)
C     VERTICAL PLOT OF ONE VECTOR, INPUT LIMITS (VPLOTL)
C     VERTICAL PLOT OF ONE VECTOR, INPUT LIMITS, CONNECTED POINTS
C        (VPLOTBL)
C     VERTICAL PLOT OF TWO VECTORS, COMPUTED LIMITS (VPLOT2)
C     VERTICAL PLOT OF TWO VECTORS, INPUT LIMITS (VPLOT2L)
C
C
C  WRITTEN BY - JANET R. DONALDSON
C               STATISTICAL ENGINEERING LABORATORY/BOULDER
C               NATIONAL BUREAU OF STANDARDS
C
C  LATEST REVISION  -  14 JULY 1978   JEK
C
C
C  INPUT ARGUMENTS - (IN ORDER OF USAGE)
C
C  X1     THE N VECTOR CONTAINING THE 1ST OBSERVED SERIES
C  X2     THE N VECTOR CONTAINING THE 2ND OBSERVED SERIES
C  M      THE INTEGER NUMBER OF OBSERVATIONS IN EACH SERIES 
C  NSER   THE INTEGER WHICH INDICATES THE NUMBER OF SERIES (VECTORS)
C         TO BE PLOTTED
C         IF NSER.EQ.1, ONE SERIES IS PLOTTED
C         IF NSER.EQ.2, TWO SERIES ARE PLOTTED
C  LMSPEC THE INTEGER VALUE USED TO DESIGNATE WHICH SET OF SPECIFIED
C         LIMITS ARE TO BE USED WHEN NSER.EQ.2
C         IF NSER.EQ.1, THIS IS A DUMMY VARIABLE
C         IF LMSPEC = 0, THE LIMITS ARE SPECIFIED BY THE MINIMUM AND
C                        MAXIMUM OF THE TWO SERIES
C         IF LMSPEC = 1, THE LIMITS ARE SPECIFIED BY THOSE DESIGNATED 
C                        FOR THE X1 VECTOR
C  IMEAN  THE INTEGER INDICATOR AS TO WHETHER THE MEAN OR ZERO VALUE
C         IS TO BE PLOTTED.
C         IF IMEAN.EQ.0 GRAPH WILL PLOT ZERO VALUE
C         IF IMEAN.NE.0 GRAPH WILL PLOT MEAN VALUE
C  IPR    THE LOGICAL UNIT WHERE OUTPUT WILL BE DIRECTED
C
C
C  ADDITIONAL VARIABLES USED
C
C  BLANK  THE PLOTTING SYMBOL  (BLANK)
C  D      THE PLOTTING SYMBOL  "+"
C  ALPHAI THE PLOTTING SYMBOL  "I"
C  IPR    THE LOGICAL UNIT NUMBER OF THE PRINTER.  ITS VALUE IS
C         INSTALLATION DEPENDENT
C  IX1    THE PLOT LINE LOCATION FOR THE ITH VALUE OF THE X1 VECTOR
C  IX2    THE PLOT LINE LOCATION FOR THE ITH VALUE OF THE X2 VECTOR
C  MEAN1  THE PLOT LINE LOCATION OF THE MEAN(OR ZERO VALUE) FOR THE X1
C         VECTOR
C  ONE    THE SYMBOL  "1"  USED FOR PLOT LABELS
C  STAR   THE PLOTTING SYMBOL  "*"
C  PL     THE VECTOR CONTAINING HORIZONTAL PRINT POSITIONS FOR EACH
C         LINE OF THE GRAPH
C  SCALE1 THE SCALE OF THE X1 VECTOR PLOT
C  SCALE2 THE SCALE OF THE X2 VECTOR PLOT
C  TWO    THE SYMBOL  "2"  USED FOR PLOT LABLES AND TO INDICATE
C         MULTIPLE POINTS ON THE GRAPH
C  XLIM   A VECTOR OF THE X AXIS (HORIZONTAL) GRAPH LABELS
C  XMN1   THE MINIMUM VALUE TO BE PLOTTED FOR X1
C  XMN2   THE MINIMUM VALUE TO BE PLOTTED FOR X2
C  XMX1   THE MAXIMUM VALUE TO BE PLOTTED FOR X1
C  XMX2   THE MAXIMUM VALUE TO BE PLOTTED FOR X2
C  X1MEAN THE MEAN VALUE OF THE X1 VECTOR
C  X2MEAN THE MEAN VALUE OF THE X2 VECTOR
C
C
C  MACHINE DEPENDENT VARIABLES DEFINED  -  NONE
C
      DIMENSION X1(M),X2(M),PL(103),XLIM(9)
      DOUBLE PRECISION SM
      CHARACTER*1 BLANK,D,ONE,TWO,STAR,PL,ALPHAI
      DATA BLANK /' '/ ,D /'+'/ ,ONE /'1'/ ,TWO /'2'/ ,ALPHAI /'I'/ , 
     2     STAR /'*'/
C
C  COMPUTE SCALE FOR X1
C
      XMN1=X1(1)
      XMX1=X1(1)
      DO 600 I=2,M
         IF (X1(I).LT.XMN1) XMN1=X1(I)
         IF (X1(I).GT.XMX1) XMX1=X1(I)
600   CONTINUE
      SCALE1=(XMX1-XMN1)/100. 
      SCALE2=SCALE1 
C
C  CHECK FOR SECOND VECTOR
C
      IF (NSER.EQ.1) GO TO 620
C
C  COMPUTE SCALE FOR X2
C
      XMN2=X2(1)
      XMX2=X2(1)
      DO 610 I=2,M
         IF (X2(I).LT.XMN2) XMN2=X2(I)
         IF (X2(I).GT.XMX2) XMX2=X2(I)
610   CONTINUE
      SCALE2=(XMX2-XMN2)/100. 
C
C
C  ONE SCALE SPECIFIED BY MINIMUM AND MAXIMUM OF BOTH SERIES
C
      IF (XMN1.GT.XMN2) XMN1=XMN2
      XMN2=XMN1
      IF (XMX1.LT.XMX2) XMX1=XMX2
      XMX2=XMX1
      SCALE1=(XMX1-XMN1)/100. 
      SCALE2=SCALE1 
C
C  COMPUTE MEAN OF X1
C
620   IF (SCALE1.GT.0.0.AND.SCALE2.GT.0.0) GO TO 640
      IF (NSER.EQ.2) GO TO 630
      WRITE (IPR,810) X1(1)
      RETURN
630   IF (LMSPEC.EQ.0) WRITE (IPR,820) X1(1)
      IF (LMSPEC.EQ.1) WRITE (IPR,810) X1(1)
      RETURN
640   SM=0.0D0
      DO 650 I=1,M
         SM=SM+DBLE(X1(I))
650   CONTINUE
      X1MEAN=SNGL(SM/DBLE(FLOAT(M)))
C
C  COMPUTE LOCATION OF ZERO IN PLOT STRING
C
      IF (IMEAN.EQ.0) MEAN1=INT(((-XMN1)/SCALE1)+2.5)
      IF (NSER.EQ.1) GO TO 670
C
C  COMPUTE MEAN OF X2
C
      SM=0.0D0
      DO 660 I=1,M
         SM=SM+DBLE(X2(I))
660   CONTINUE
      X2MEAN=SNGL(SM/DBLE(FLOAT(M)))
C
C  PRINT HEADINGS FOR X1 (AND X2)
C
670   IF (NSER.EQ.1) WRITE (IPR,800) X1MEAN
      IF (NSER.EQ.2) WRITE (IPR,830) ONE,X1MEAN,TWO,X2MEAN
      WRITE (IPR,740) SCALE1
      WRITE (IPR,750)
      DO 680 I=1,9
         XLIM(I)=XMN1+SCALE1*FLOAT(I)*10.
680   CONTINUE
      WRITE (IPR,760) XMN1,(XLIM(I),I=2,8,2),XMX1,(XLIM(I),I=1,9,2)
C
C  BEGIN PLOTTING
C
      IF (NSER.EQ.2) GO TO 710
C
C  ONE VECTOR
C
      WRITE (IPR,770)
      LABEL=1
      DO 700 I=1,M
         DO 690 J=1,103
            PL(J)=BLANK
690      CONTINUE
         IX1=INT(((X1(I)-XMN1)/SCALE1)+2.5)
         IF (IX1.LT.2) IX1=1
         IF (IX1.GT.102) IX1=103
         IF (IMEAN.EQ.0) PL(MEAN1)=ALPHAI
         PL(IX1)=STAR
         WRITE (IPR,840) LABEL,(PL(J),J=1,103),X1(I)
         LABEL=LABEL+1
700   CONTINUE
      WRITE (IPR,850)
      RETURN
C
C  TWO VECTORS
C
710   WRITE (IPR,790)
      LABEL=1
      DO 730 I=1,M
         DO 720 J=1,103
            PL(J)=BLANK
720      CONTINUE
         IX1=INT(((X1(I)-XMN1)/SCALE1)+2.5)
         IF (IX1.LT.2.OR.IX1.GT.102) IX1=103
         PL(IX1)=STAR
         IX2=INT(((X2(I)-XMN2)/SCALE2)+2.5)
         IF (IX2.LT.2.OR.IX2.GT.102) IX2=103
         PL(IX2)=D
         IF (IX1.EQ.IX2) PL(IX1)=TWO
         WRITE (IPR,780) LABEL,(PL(J),J=1,103),X1(I),X2(I)
         LABEL=LABEL+1
730   CONTINUE
      WRITE (IPR,850)
      RETURN
C
C  FORMAT STATEMENTS
C
740   FORMAT (//19H GRAPH INTERVAL IS ,G15.8)
750   FORMAT (10X,41H-COMPUTED FROM MAXIMUM AND MINIMUM VALUES)
760   FORMAT (2X,6(G15.8,5X)/12X,5(G15.8,5X)/9X,10(10HI---------),1HI)
770   FORMAT (8X,1HI,101X,1HI)
780   FORMAT (1X,I5,1X,A1,1HI,101A1,1HI,A1,2(1X,E9.3))
790   FORMAT (8X,1HI,101X,1HI,5X,6HX1 = *,4X,6HX2 = +)
800   FORMAT (' '/9H MEAN IS ,G15.8)
810   FORMAT(/53H THE VALUES IN THE X1 ARRAY ARE ALL IDENTICALLY EQUAL/
     2        4H TO ,E13.7)
820   FORMAT(/54H THE VALUES IN THE X1 AND X2 ARRAY ARE ALL IDENTICALLY/
     2        10H EQUAL TO ,E13.7)
830   FORMAT (' '/2H X,A1,7H VECTOR/9H MEAN IS ,G15.8)
840   FORMAT (1X,I5,1X,A1,1HI,101A1,1HI,A1,1X,E11.5)
850   FORMAT (1H ,8X,10(10HI---------),1HI)
C
C
C
      END 
C***********************************************************************
C
      SUBROUTINE VMAINH (X1,X2,M,NSER,LMSPEC,IMEAN,IPR)
C
C***********************************************************************
C
C  THIS SUBROUTINE IS IDENTICAL TO VPMAIN EXCEPT THAT IT PRODUCES
C  PLOTS ONE-HALF THE SIZE OF VPMAIN.
C
C  FOR DETAILED COMMENTS AND A LIST OF THE INPUT ARGUMENTS AND
C  VARIABLES, SEE THE COMMENTS AT THE BEGINNING OF VPMAIN BY
C  JANET R. DONALDSON.
C
C
      DIMENSION X1(M),X2(M),PL(53),XLIM(9)
      DOUBLE PRECISION SM
      CHARACTER*1 BLANK,D,ONE,TWO,STAR,PL,ALPHAI
      DATA BLANK /' '/ ,D /'+'/ ,ONE /'1'/ ,TWO /'2'/ ,ALPHAI /'I'/ , 
     2     STAR /'*'/
C
C  COMPUTE SCALE FOR X1
C
      XMN1=X1(1)
      XMX1=X1(1)
      DO 860 I=2,M
         IF (X1(I).LT.XMN1) XMN1=X1(I)
         IF (X1(I).GT.XMX1) XMX1=X1(I)
860   CONTINUE
      SCALE1=(XMX1-XMN1)/50.
      SCALE2=SCALE1 
C
C  CHECK FOR SECOND VECTOR
C
      IF (NSER.EQ.1) GO TO 880
C
C  COMPUTE SCALE FOR X2
C
      XMN2=X2(1)
      XMX2=X2(1)
      DO 870 I=2,M
         IF (X2(I).LT.XMN2) XMN2=X2(I)
         IF (X2(I).GT.XMX2) XMX2=X2(I)
870   CONTINUE
      SCALE2=(XMX2-XMN2)/50.
C
C  ONE SCALE SPECIFIED BY MINIMUM AND MAXIMUM OF BOTH SERIES
C
      IF (XMN1.GT.XMN2) XMN1=XMN2
      XMN2=XMN1
      IF (XMX1.LT.XMX2) XMX1=XMX2
      XMX2=XMX1
      SCALE1=(XMX1-XMN1)/50.
      SCALE2=SCALE1 
      GO TO 880
C
C  COMPUTE MEAN OF X1
C
880   IF (SCALE1.GT.0.0.AND.SCALE2.GT.0.0) GO TO 900
      IF (NSER.EQ.2) GO TO 890
      WRITE (IPR,1070) X1(1)
      RETURN
890   IF (LMSPEC.EQ.0) WRITE (IPR,1080) X1(1)
      IF (LMSPEC.EQ.1) WRITE (IPR,1070) X1(1)
      RETURN
900   SM=0.0D0
      DO 910 I=1,M
         SM=SM+DBLE(X1(I))
910   CONTINUE
      X1MEAN=SNGL(SM/DBLE(FLOAT(M)))
C
C  COMPUTE LOCATION OF ZERO IN PLOT STRING
C
      IF (IMEAN.EQ.0) MEAN1=INT(((-XMN1)/SCALE1)+2.5)
      IF (NSER.EQ.1) GO TO 930
C
C  COMPUTE MEAN OF X2
C
      SM=0.0D0
      DO 920 I=1,M
         SM=SM+DBLE(X2(I))
920   CONTINUE
      X2MEAN=SNGL(SM/DBLE(FLOAT(M)))
C
C  PRINT HEADINGS FOR X1 (AND X2)
C
930   IF (NSER.EQ.1) WRITE (IPR,1060) X1MEAN
      IF (NSER.EQ.2) WRITE (IPR,1090) ONE,X1MEAN,TWO,X2MEAN 
      WRITE (IPR,1000) SCALE1 
      WRITE (IPR,1010)
      DO 940 I=1,9
         XLIM(I)=XMN1+SCALE1*FLOAT(I)*5.
940   CONTINUE
      WRITE (IPR,1020) XMN1,(XLIM(I),I=2,8,2),XMX1,(XLIM(I),I=1,9,2)
C
C  BEGIN PLOTTING
C
      IF (NSER.EQ.2) GO TO 970
C
C  ONE VECTOR
C
      WRITE (IPR,1030)
      LABEL=1
      DO 960 I=1,M
         DO 950 J=1,53
            PL(J)=BLANK
950      CONTINUE
         IX1=INT(((X1(I)-XMN1)/SCALE1)+2.5)
         IF (IX1.LT.2) IX1=1
         IF (IX1.GT.52) IX1=53
         IF (IMEAN.EQ.0) PL(MEAN1)=ALPHAI
         PL(IX1)=STAR
         WRITE (IPR,1100) LABEL,(PL(J),J=1,53),X1(I)
         LABEL=LABEL+1
960   CONTINUE
      WRITE (IPR,1110)
      RETURN
C
C  TWO VECTORS
C
970   WRITE (IPR,1050)
      LABEL=1
      DO 990 I=1,M
         DO 980 J=1,53
            PL(J)=BLANK
980      CONTINUE
         IX1=INT(((X1(I)-XMN1)/SCALE1)+2.5)
         IF (IX1.LT.2.OR.IX1.GT.52) IX1=53
         PL(IX1)=STAR
         IX2=INT(((X2(I)-XMN2)/SCALE2)+2.5)
         IF (IX2.LT.2.OR.IX2.GT.52) IX2=53
         PL(IX2)=D
         IF (IX1.EQ.IX2) PL(IX1)=TWO
         WRITE (IPR,1040) LABEL,(PL(J),J=1,53),X1(I),X2(I)
         LABEL=LABEL+1
990   CONTINUE
      WRITE (IPR,1110)
      RETURN
C
C  FORMAT STATEMENTS
C
1000  FORMAT (//19H GRAPH INTERVAL IS ,G15.8)
1010  FORMAT (10X,41H-COMPUTED FROM MAXIMUM AND MINIMUM VALUES)
1020  FORMAT (1X,6(F9.2,1X)/6X,5(F9.2,1X)/7X,10(5HI----),1HI)
1030  FORMAT (6X,1HI,51X,1HI) 
1040  FORMAT (1X,I3,1X,A1,1HI,51A1,1HI,A1,E9.3,1X,E9.3)
1050  FORMAT (6X,1HI,51X,1HI,3X,6HX1 = *,4X,6HX2 = +)
1060  FORMAT (' '/9H MEAN IS ,G15.8)
1070  FORMAT(/53H THE VALUES IN THE X1 ARRAY ARE ALL IDENTICALLY EQUAL/
     2        4H TO ,E13.7)
1080  FORMAT(/54H THE VALUES IN THE X1 AND X2 ARRAY ARE ALL IDENTICALLY/
     2        10H EQUAL TO ,E13.7)
1090  FORMAT (' ',2H0X,A1,7H VECTOR/9H MEAN IS ,G15.8)
1100  FORMAT (1X,I3,1X,A1,1HI,51A1,1HI,A1,1X,E11.5)
1110  FORMAT (1H ,6X,10(5HI----),1HI)
      END 
C
C
      SUBROUTINE STAT (M,N,NKG2,KG,Y,WEIGHT,C,ALF,A,IV,T,ETA,R,ISYM,XPL,
C
     2                 YPL)
C
C***********************************************************************
C
C   PROVIDE THE INDICATED INFORMATION IN RESPONSE TO
C   THE RS COMMAND: 
C  (1)  STANDARD DEVIATION OF THE PARAMETERS
C  (2)  T-RATIO
C  (3)  RESIDUAL SUM OF SQUARES
C  (4)  RESIDUAL STANDARD ERROR
C  (5)  COEFFICIENT OF DETERMINATION
C  (6)  STANDARDIZED RESIDUALS
C  (7)  PLOT OF STANDARDIZED RESIDUALS VS INDENDENT
C       VARIABLES
C  (8)  PLOT OF STANDARDIZED RESIDUALS VS PREDICTED
C       VALUES
C  (9)  NORMAL PROBABILITY PLOT.
C
      PARAMETER ( MMAX=300, MAXN1=11, MAXNKG=18 ) 
      CHARACTER*1 IANS
      CHARACTER*4 IND,VAR
      REAL Y(M),C(N),ALF(KG),A(MMAX,NKG2)
      REAL T(MMAX,IV),WEIGHT(M),ETA(M),R(M)
      REAL BM(MMAX,MAXNKG),XX1(MMAX,MAXNKG),PH(MAXN1)
      REAL SYH(MMAX),SR(MMAX),SRV(MMAX) 
      LOGICAL IP,IPT,IR,ITB,IPIV,IPPV,IPP,TS,ITP
      DIMENSION IOR(MMAX)
      DIMENSION ISYM(1),XPL(1),YPL(1)
      COMMON /IPIBLK/ IPI
      COMMON /BLOKIP/ IP
      COMMON /BLKPH/ PH
      COMMON /BLKBM/ BM
      COMMON /BLKTS/ TS
      COMMON /BLKITP/ ITP
      DATA PI /3.141592654/
      DATA IND /'IND '/ ,VAR /'VAR '/
      DATA IPR6 /6/ ,IPR10 /10/
C
C  INITIALIZE VALUES
C
      IPI=1
      NP1=N+1
      NPKG=N+KG
      RMS=A(2,NKG2) 
C
C  DETERMINE WHICH STATISTICS ARE TO BE PRINTED
C
      PRINT 10
10    FORMAT (' INDICATE WHICH OF THE FOLLOWING SHOULD',' BE PRINTED:'/
     2        '   PARAMETERS, THEIR STD DEV',' & T-RATIO:  ')
      READ (5,20) IANS
20    FORMAT (A1)
      IF (IANS.EQ.'Y') IPT=.TRUE.
      PRINT 30
30    FORMAT ('   RES SUM OF SQS, RES MEAN SQ, RES ','STD DEV &'/
     2        '       COEFF OF DETERMINATION:  ') 
      READ (5,20) IANS
      IF (IANS.EQ.'Y') IR=.TRUE.
      PRINT 40
40    FORMAT ('   TABLE OF OBS, PRED VAL, RESID & STD',' RESID:  ')
      READ (5,20) IANS
      IF (IANS.EQ.'Y') ITB=.TRUE.
      PRINT 50
50    FORMAT ('   PLOT OF STD RES VS IND VAR:  ') 
      READ (5,20) IANS
      IF (IANS.EQ.'Y') IPIV=.TRUE.
      PRINT 60
60    FORMAT ('   PLOT OF STD RES VS PRED VAL:  ')
      READ (5,20) IANS
      IF (IANS.EQ.'Y') IPPV=.TRUE.
      PRINT 70
70    FORMAT ('   NORMAL PROBABILITY PLOT:  ')
      READ (5,20) IANS
      IF (IANS.EQ.'Y') IPP=.TRUE.
      IF (.NOT.IPT) GO TO 130 
C
C  PRINT PARAMETERS, THEIR STD DEV & T-RATIOS
C
      IPT=.FALSE.
      PRINT 80
      IF (IP) WRITE (10,80)
80    FORMAT ('1'/' '/
     1          4X,'PARAMETER',6X,'S.D. OF PARM.',5X,'T-RATIO'/)
      IF (N.EQ.0) GO TO 110
      DO 100 I=1,N
         AI=ABS(A(I,I))
         BB=SQRT(AI)
         TR=C(I)/BB 
         PRINT 90, C(I),BB,TR 
         IF (IP) WRITE (10,90) C(I),BB,TR
90       FORMAT (' ',1P3E16.6)
100   CONTINUE
110   IF (KG.EQ.0) GO TO 130
      DO 120 I=1,KG 
         J=N+I
         AJ=ABS(A(J,J))
         BB=SQRT(AJ)
         TR=ALF(I)/BB
         PRINT 90, ALF(I),BB,TR
120      IF (IP) WRITE (10,90) ALF(I),BB,TR
130   IF (.NOT.IR) GO TO 200
C
C  PRINT RESIDUAL SUM OF SQUARES, RES MEAN SQUARE &
C     RES STD ERROR.
C  PRINT COEFFICIENT OF DETERMINATION.
C
      IR=.FALSE.
      RSS=(A(1,NKG2))**2
      RSD=SQRT(ABS(RMS))
      PRINT 140, RSS
      IF (IP) WRITE (10,140) RSS
140   FORMAT (' '/' '/' WTD RESIDUAL SUM OF SQUARES:  ',1PE16.6)
      PRINT 150, RMS
      IF (IP) WRITE (10,150) RMS
150   FORMAT (' '/' ','WTD RESIDUAL MEAN SQUARE:  ',3X,1PE16.6)
      PRINT 160, RSD
      IF (IP) WRITE (10,160) RSD
160   FORMAT (' '/' ','WTD RESIDUAL STANDARD ERROR:  ',1PE16.6//)
      YB=0.
      SWS=0.
      DO 170 I=1,M
         SWS=SWS+WEIGHT(I)
170      YB=YB+Y(I)*WEIGHT(I) 
      YB=YB/SWS
      CTSS=0.
      DO 180 I=1,M
180      CTSS=CTSS+WEIGHT(I)*(Y(I)-YB)**2
      R2=1.-RSS/CTSS
      PRINT 190, R2 
      IF (IP) WRITE (10,190) R2
190   FORMAT (' ','COEFFICIENT OF DETERMINATION (R-SQUARE): ',1PE16.6//)
C
C  CALCULATE STANDARDIZED RESIDUALS
C
200   IF (N.EQ.0) GO TO 220
      DO 210 I=1,M
         II=I
         CALL CALC (II,KG,N,MMAX,IV,C,ALF,T,DUMMY)
         DO 210 J=1,N
210         BM(I,J)=PH(J)
220   DO 240 J=1,NPKG
         DO 240 I=1,M
            SUM=0.
            DO 230 K=1,NPKG
230            SUM=SUM+BM(I,K)*A(K,J)
240         XX1(I,J)=SUM
      DO 260 I=1,M
         SUM=0.
         DO 250 K=1,NPKG
250         SUM=SUM+XX1(I,K)*BM(I,K)
         SYH(I)=SQRT(ABS(SUM))
260      SR(I)=R(I)/SQRT(ABS(RMS/WEIGHT(I)-SUM))
      IF (.NOT.ITB) GO TO 330 
C
C  PRINT TABLE OF OBS & PRED VALUES ALONG WITH STD RES
C
      ITB=.FALSE.
      PRINT 270
      IF (IP) WRITE (10,270)
270   FORMAT ('1'/' INDEX',2X,'IND VAR 1',4X,'OBS VAL',5X,'PRED VAL',6X,
     2        'RESID',8X,'STD RES'/)
      DO 290 I=1,M
         PRINT 280, I,T(I,1),Y(I),ETA(I),R(I),SR(I)
         IF (IP) WRITE (10,280) I,T(I,1),Y(I),ETA(I),R(I),SR(I)
280      FORMAT (' ',I3,1P5E13.5)
290   CONTINUE
C
330   IF (.NOT.IPIV) GO TO 370
C
C  PLOT STD RES VS IND VAR
C
      IPIV=.FALSE.
      IPI=1
      DO 360 I=1,IV 
         IF (ITP) PRINT 340, I
         IF (IP) WRITE (10,340) I
340      FORMAT ('1',15X,'STANDARDIZED RESIDUAL', 
     2           ' VS INDEPENDENT VARIABLE',I2/)
         IF (TS.AND.ITP) CALL VPLT (SR,M,IPR6,0)
         IF (TS.AND.IP) CALL VPLT (SR,M,IPR10,0)
         IF (.NOT.TS) CALL SETUP1 (T,M,I,SR,ISYM,XPL,YPL)
360   CONTINUE
370   IF (.NOT.IPPV) GO TO 390
C
C  PLOT STD RES VS PRED VAL
C
      IPPV=.FALSE.
      IF (ITP) PRINT 380
      IF (IP) WRITE (10,380)
380   FORMAT ('1',15X,'STANDARDIZED RESIDUAL VS. PREDICTED VALUES'/)
      CALL SETUP1 (ETA,M,1,SR,ISYM,XPL,YPL)
390   IF (.NOT.IPP) GO TO 460 
C
C  DO NORMAL PROBABILITY PLOT 
C
      IPP=.FALSE.
      IF (M.GT.10) GO TO 410
      DO 400 I=1,M
         PP=(3.*I-1.)/(3.*M+1.)
400      SYH(I)=4.91*(PP**(.14)-(1.-PP)**(.14))
      GO TO 430
410   PP=2./(3.*M+1.)
      SYH(1)=4.91*(PP**(.14)-(1.-PP)**(.14))
      PI8=PI/8
      PI4=PI/4
      DO 420 I=2,M
         PP=(I-PI8)/(M+1.-PI4)
         SYH(I)=4.91*(PP**(.14)-(1.-PP)**(.14))
420   CONTINUE
430   CALL ORPLOT (SR,M,IOR)
      DO 440 I=1,M
440      SRV(I)=SR(IOR(I))
      IF (ITP) PRINT 450
      IF (IP) WRITE (10,450)
450   FORMAT ('1',15X,'NORMAL PROBABILITY PLOT'/) 
      CALL SETUP1 (SYH,M,1,SRV,ISYM,XPL,YPL)
460   RETURN
      END 
C
C
C***********************************************************************
C
      SUBROUTINE VARPRO (N,KG,M,MMAX,NPP2,IV,T,Y,W,WEIGHT,A,IPRINT,ALF,C
     2,IERR,ITMAX,EPS1)
C
C***********************************************************************
C
C  THE FOLLOWING COMMENTS MOSTLY REFER TO
C  THE BATCH VERSION OF VARPRO.
C
C  GIVEN A SET OF M OBSERVATIONS, CONSISTING OF
C  VALUES Y(1), Y(2),...,Y(M) OF A DEPENDENT VARIABLE
C  Y, WHERE Y(I) CORRESPONDS TO THE IV INDEPENDENT
C  VARIABLE(S) T(I,1), T(I,2),..., T(I,IV), VARPRO
C  ATTEMPTS TO COMPUTE A WEIGHTED LEAST SQUARES FIT TO
C  A FUNCTION ETA (THE 'MODEL') WHICH IS A LINEAR 
C  COMBINATION
C                 N 
C    ETA(C,ALF;T)  =  SUM C * PHI (ALF;T) + PHI  (ALF;T)
C                J=1  J     J             N+1
C
C  OF NONLINEAR FUNCTIONS PHI(J) (E.G. A SUM OF
C  EXPONENTIALS AND/OR GAUSSIANS),  THAT IS, DETERMINE
C  THE LINEAR PARAMETERS C(J) AND THE VECTOR OF
C  NONLINEAR PARAMETERS ALF BY MINIMIZING
C
C                2     M                           2
C  NORM(RESIDUAL)  =  SUM WEIGHT * (Y - ETA(C,ALF;T )) .
C                     I=1       I    I             I
C
C  THE (N+1)-ST TERM IS OPTIONAL, AND IS USED WHEN IT
C  IS DESIRED TO FIX ONE OR MORE OF THE C'S (RATHER
C  THAN LET THEM BE DETERMINED).  VARPRO REQUIRES 
C  FIRST DERIVATIVES OF THE PHI'S.
C
C                  NOTES: 
C
C  A)  THE ABOVE PROBLEM IS ALSO REFERRED TO AS
C  'MULTIPLE NONLINEAR REGRESSION'.  FOR USE IN
C  STATISTICAL ESTIMATION, VARPRO RETURNS THE
C  RESIDUALS, THE COVARIANCE MATRIX OF THE LINEAR AND
C  NONLINEAR PARAMETERS, AND THE ESTIMATED VARIANCE
C  OF THE OBSERVATIONS.  OTHER OPTIONS ARE AVAILABLE
C  IN THE INTERACTIVE VERSION.
C
C  B)  AN ETA OF THE ABOVE FORM IS CALLED 'SEPARABLE'.
C  THE CASE OF A NONSEPARABLE ETA CAN BE HANDLED BY
C  SETTING N = 0 AND USING PHI(N+1).
C
C  C)  VARPRO MAY ALSO BE USED TO SOLVE LINEAR LEAST
C  SQUARES PROBLEMS (IN THAT CASE NO ITERATIONS ARE
C  PERFORMED).  SET KG = 0.
C
C  D)  THE MAIN ADVANTAGE OF VARPRO OVER OTHER LEAST
C  SQUARES PROGRAMS IS THAT NO INITIAL GUESSES ARE
C  NEEDED FOR THE LINEAR PARAMETERS.  NOT ONLY DOES
C  THIS MAKE IT EASIER TO USE, BUT IT OFTEN LEADS 
C  TO FASTER CONVERGENCE.
C
C
C   DESCRIPTION OF PARAMETERS 
C
C  N       NUMBER OF LINEAR PARAMETERS C (MUST BE 
C          .GE. 0). 
C  KG      NUMBER OF NONLINEAR PARAMETERS ALF (MUST
C          BE .GE. 0).
C  M       NUMBER OF OBSERVATIONS, M MUST BE GREATER
C          THAN KG + N (I.E., THE NUMBER OF OBSERVA-
C          TIONS MUST EXCEED THE NUMBER OF PARAMETERS).
C  IV      NUMBER OF INDEPENDENT VARIABLES T.
C  T       REAL M BY IV ARRAY OF INDEPENDENT VARIABLES.
C          T(I,J) CONTAINS THE VALUE OF THE I-TH
C          OBSERVATION OF THE J-TH INDEPENDENT
C          VARIABLE.
C  Y       M-VECTOR OF OBSERVATIONS, ONE FOR EACH 
C          ROW OF T.
C  W       M-VECTOR OF SQRT OF WEIGHTS. 
C  WEIGHT  M-VECTOR OF NONNEGATIVE WEIGHTS.  SHOULD
C          BE SET TO 1'S IF WEIGHTS ARE NOT DESIRED.
C  INC     KG X (N+1) INTEGER INCIDENCE MATRIX.
C          INC(K,J) = 1 IF NONLINEAR PARAMETER ALF(K)
C          APPEARS IN THE J-TH FUNCTION PHI(J),
C          OTHERWISE, INC(K,J) = 0.  IF PHI(N+1) IS
C          NOT INCLUDED IN THE MODEL, THE (N+1)-ST
C          COLUMN OF INC SHOULD BE SET TO ZERO.  IF
C          PHI(N+1) IS INCLUDED, THE APPROPRIATE
C          ELEMENTS OF THE (N+1)-ST COLUMN SHOULD BE
C          SET TO 1'S.  INC IS NOT NEEDED WHEN N = 0 OR
C          KG = 0.  CAUTION:  THE DECLARED ROW DIMENSION
C          OF INC (IN ADA) MUST CURRENTLY BE SET TO 12.
C          SEE 'RESTRICTIONS' BELOW.
C  MMAX    THE DECLARED ROW DIMENSION OF THE MATRICES
C          A AND T.  IT MUST BE AT LEAST MAX(M,2*KG+3).
C  NPP2    N+P+2, WHERE P IS THE NUMBER OF ONES IN
C          MATRIX INC..  THE DECLARED COLUMN
C          DIMENSION OF A MUST BE AT LEAST NPP2.
C          (IF N = 0, SET NPP2 = KG+2.  IF KG = 0,
C          SET NPP2 = N+2.)
C  A       REAL MATRIX OF SIZE MAX(M,2*KG+3) BY N+P+2.
C          ON INPUT IT CONTAINS THE PHI(J)'S AND
C          THEIR DERIVATIVES (SEE BELOW).  ON OUTPUT
C          THE FIRST N+KG ROWS AND COLUMNS OF A WILL
C          CONTAIN AN APPROXIMATION TO THE (WEIGHTED)
C          COVARIANCE MATRIX AT THE SOLUTION (THE 
C          FIRST N ROWS CORRESPOND TO THE LINEAR
C          PARAMETERS, THE LAST KG TO THE NONLINEAR
C          ONES), COLUMN N+KG+1 WILL CONTAIN THE
C          WEIGHTED RESIDUALS (Y - ETA), A(1,N+KG+2)
C          WILL CONTAIN THE (EUCLIDEAN) NORM OF THE
C          WEIGHTED RESIDUAL, AND A(2,N+KG+2) WILL
C          CONTAIN AN ESTIMATE OF THE (WEIGHTED)
C          VARIANCE OF THE OBSERVATIONS,
C          NORM(RESIDUAL)**2/(M-N-KG).
C  IPRINT  INPUT INTEGER CONTROLLING PRINTED OUTPUT.
C          IF IPRINT IS POSITIVE, THE NONLINEAR
C          PARAMETERS, THE WEIGHTED NORM OF THE RESIDUAL, AND
C          THE MARQUARDT PARAMETER WILL BE OUTPUT 
C          EVERY IPRINT-TH ITERATION (AND INITIALLY,
C          AND AT THE FINAL ITERATION).  THE LINEAR
C          PARAMETERS WILL BE PRINTED AT THE FINAL
C          ITERATION.  ANY ERROR MESSAGES WILL ALSO
C          BE PRINTED.  (IPRINT = 1 IS RECOMMENDED
C          AT FIRST.)  IF IPRINT = 0,  ONLY THE FINAL
C          QUANTITIES WILL BE PRINTED, AS WELL AS ANY
C          ERROR MESSAGES.  IF IPRINT = -1, NO
C          PRINTING WILL BE DONE.  THE USER IS THEN
C          RESPONSIBLE FOR CHECKING THE PARAMETER 
C          IERR FOR ERRORS.
C  ALF     KG-VECTOR OF ESTIMATES OF NONLINEAR
C          PARAMETERS (INPUT).  ON OUTPUT IT WILL 
C          CONTAIN OPTIMAL VALUES OF THE NONLINEAR
C          PARAMETERS.
C  C       N-VECTOR OF LINEAR PARAMETERS (OUTPUT ONLY).
C  IERR    INTEGER ERROR FLAG (OUTPUT): 
C          .GT. 0 - SUCCESSFUL CONVERGENCE. IERR IS
C          THE NUMBER OF ITERATIONS TAKEN.
C          -1  TERMINATED FOR TOO MANY ITERATIONS.
C          -2  TERMINATED FOR ILL-CONDITIONING
C              (MARQUARDT PARAMETER TOO LARGE.) ALSO
C              SEE IERR = -8 BELOW.
C          -4  INPUT ERROR IN PARAMETER N, M, KG, 
C              NPP2 OR MMAX.
C          -5  INC MATRIX IMPROPERLY SPECIFIED, OR
C              P DISAGREES WITH NPP2.
C          -6  WEIGHTS NOT NONNEGATVE.
C          -7  'CONSTANT' COLUMN WAS COMPUTED MORE
C              THAN ONCE.
C          -8  CATASTROPHIC FAILURE - A COLUMN OF THE
C              A MATRIX HAS BECOME ZERO.  SEE
C              'CONVERGENCE FAILURES' BELOW.
C
C          (IF IERR .LE. -4, THE LINEAR PARAMETERS,
C          COVARIANCE MATRIX, ETC. ARE NOT RETURNED.)
C
C   SUBROUTINES REQUIRED
C
C  SEVEN SUBROUTINES, DPA, ORFAC1, ORFAC2, BACSUB,
C  POSTPR, COV, VARERR AND ADA ARE PROVIDED.  IN ADDITION,
C  THE USER MUST PROVIDE SUBROUTINES CALC, DERIV AND WT,
C  WHICH, GIVEN ALF, WILL EVALUATE THE FUNCTIONS PHI(J) AND 
C  THEIR PARTIAL DERIVATIVES D PHI(J)/D ALF(K), AT THE SAMPLE
C  POINTS T(I).
C
C  RESTRICTIONS
C
C        THE SUBROUTINES DPA AND ADA CONTAIN THE LOCALLY
C     DIMENSIONED MATRIX INC, WHOSE DIMENSIONS ARE
C     CURRENTLY SET FOR MAXIMA OF N+1 = 15, KG = 12.
C     THEY MUST BE CHANGED FOR LARGER PROBLEMS.
C     DATA PLACED IN ARRAY A IS OVERWRITTEN ('DESTROYED').
C     DATA PLACED IN ARRAYS T, Y AND INC IS LEFT INTACT.
C     THE PROGRAM RUNS IN WATFIV, EXCEPT WHEN N = 0
C     OR KG = 0.
C
C        IT IS ASSUMED THAT THE MATRIX PHI(J, ALF; T(I))
C     HAS FULL COLUMN RANK.  THIS MEANS THAT THE FIRST N
C     COLUMNS OF THE MATRIX A MUST BE LINEARLY
C     INDEPENDENT.
C
C  CONVERGENCE FAILURES
C
C        IF CONVERGENCE FAILURES OCCUR, FIRST CHECK
C     FOR INCORRECT REPLIES TO WHICH TERMS OF THE 
C     MODEL DEPEND ON WHICH NONLINEAR PARAMETER.
C     IF THESE ARE CORRECT,
C     TRY SEVERAL STARTING GUESSES FOR ALF.  IF ADA
C     IS CODED CORRECTLY, AND IF ERROR RETURNS IERR = -2
C     OR -8 PERSISTENTLY OCCUR, THIS IS A SIGN OF ILL-
C     CONDITIONING, WHICH MAY BE CAUSED BY SEVERAL THINGS.
C     ONE IS POOR SCALING OF THE PARAMETERS; ANOTHER IS
C     AN UNFORTUNATE INITIAL GUESS FOR THE PARAMETERS,
C     STILL ANOTHER IS A POOR CHOICE OF THE MODEL.
C
C  ALGORITHM
C
C        THE RESIDUAL R IS MODIFIED TO INCORPORATE, FOR
C     ANY FIXED ALF, THE OPTIMAL LINEAR PARAMETERS FOR
C     THAT ALF.  IT IS THEN POSSIBLE TO MINIMIZE ONLY ON
C     THE NONLINEAR PARAMETERS.  AFTER THE OPTIMAL VALUES
C     OF THE NONLINEAR PARAMETERS HAVE BEEN DETERMINED,
C     THE LINEAR PARAMETERS CAN BE RECOVERED BY LINEAR
C     LEAST SQUARES TECHNIQUES (SEE REF. 1).
C
C        THE MINIMIZATION IS BY A MODIFICATION OF 
C     OSBORNE'S (REF. 3) MODIFICATION OF THE LEVENBERG-
C     MARQUARDT ALGORITHM.  INSTEAD OF SOLVING THE NORMAL
C     EQUATIONS WITH MATRIX
C
C           T      2
C         (J J + NU  * D),    WHERE  J = D(ETA)/D(ALF),
C
C     STABLE ORTHOGONAL (HOUSEHOLDER) REFLECTIONS ARE
C     USED ON A MODIFICATION OF THE MATRIX
C
C                           (   J  )
C                           (------) ,
C                           ( NU*D )
C
C     WHERE D IS A DIAGONAL MATRIX CONSISTING OF THE
C     LENGTHS OF THE COLUMNS OF J.  THIS MARQUARDT
C     STABILIZATION ALLOWS THE ROUTINE TO RECOVER FROM
C     SOME RANK DEFICIENCIES IN THE JACOBIAN.  OSBORNE'S
C     EMPIRICAL STRATEGY FOR CHOOSING THE MARQUARDT
C     PARAMETER HAS PROVEN REASONABLY SUCCESSFUL IN
C     PRACTICE.  (GAUSS-NEWTON WITH STEP CONTROL CAN BE
C     OBTAINED BY MAKING THE CHANGE INDICATED BEFORE THE
C     INSTRUCTION LABELED 1).  A DESCRIPTION CAN BE FOUND
C     IN REF. (3), AND A FLOW CHART IN (2), P. 22.
C
C     FOR REFERENCE, SEE
C
C     1.  GENE H. GOLUB AND V. PEREYRA, 'THE DIFFEREN-
C         TIATION OF PSEUDO-INVERSES AND NONLINEAR LEAST
C         SQUARES PROBLEMS WHOSE VARIABLES SEPARATE,'
C         SIAM J. NUMER. ANAL. 10, 413-432 (1973).
C     2.  ------, SAME TITLE, STANFORD C.S. REPORT 72-261,
C         FEB. 1972.
C     3.  OSBORNE, MICHAEL R., 'SOME ASPECTS OF NON-
C         LINEAR LEAST SQUARES CALCULATIONS,' IN LOOTSMA,
C         ED., 'NUMERICAL METHODS FOR NON-LINEAR
C         OPTIMIZATION,' ACADEMIC PRESS, LONDON, 1972.
C     4.  KROGH, FRED, 'EFFICIENT IMPLEMENTATION OF A
C         VARIABLE PROJECTION ALGORITHM FOR NONLINEAR
C         LEAST SQUARES PROBLEMS,' COMM. ACM 17, PP. 167-
C         169 (MARCH, 1974).
C     5.  KAUFMAN, LINDA, 'A VARIABLE PROJECTION METHOD
C         FOR SOLVING SEPARABLE NONLINEAR LEAST SQUARES
C         PROBLEMS', B.I.T. 15,49-57 (1975).
C
C               JOHN BOLSTAD
C               COMPUTER SCIENCE DEPT., SERRA HOUSE
C               STANFORD UNIVERSITY
C               SEPTEMBER, 1976
C
C  ......................................................
C
      REAL A(MMAX,NPP2),C(N),ALF(KG),T(MMAX,IV),WEIGHT(M),W(M),Y(M),NU
      INTEGER B1,OUTPUT
      LOGICAL SKIP,IP
      COMMON /BLOKIP/ IP
      DATA OUTPUT /6/
C
C  THE FOLLOWING THREE PARAMETERS ARE USED IN
C  THE CONVERGENCE TEST:  EPS1 IS AN ABSOLUTE
C  TOLERANCE FOR THE NORM OF THE PROJECTION OF
C  THE RESIDUAL ONTO THE RANGE OF THE JACOBIAN
C  OF THE VARIABLE PROJECTION FUNCTIONAL.  EPS2
C  IS A RELATIVE AND ABSOLUTE TOLERANCE FOR THE
C  CHANGE IN THE NORM OF THE RESIDUAL.  ITMAX IS
C  THE MAXIMUM NUMBER OF FUNCTION AND DERIVATIVE
C  EVALUATIONS ALLOWED.  CAUTION:  EPS1 AND
C  EPS2 MUST NOT BE SET SMALLER THAN 10 TIMES
C  THE UNIT ROUND-OFF OF THE MACHINE.
C
      EPS2=EPS1
      IERR=1
      ITER=0
      NP1=N+1
      B1=N+2
      NKG2=N+KG+2
      KGP1=KG+1
      KG23=2*KG+3
      SKIP=.FALSE.
      MODIT=IPRINT
      IF (IPRINT.LE.0) MODIT=ITMAX+2
      NU=0.
C
C  IF GAUSS-NEWTON IS DESIRED, REMOVE THE
C  NEXT STATEMENT.
C
      NU=1.
C
C  CHECK FOR VALID INPUT
C
      IF (N.GE.0.AND.KG.GE.0.AND.N+KG.LT.M.AND.NKG2.LE.NPP2.AND.
     2    KG23.LE.MMAX.AND.M.LE.MMAX.AND.IV.GT.0.AND..NOT.(KG.EQ.0.AND.
     3    N.EQ.0)) GOTO 710
      IERR=-4
      CALL VARERR (IPRINT,IERR,1)
      GO TO 880
C
C           BEGIN OUTER ITERATION LOOP TO UPDATE ALF.
C           CALCULATE THE NORM OF THE RESIDUAL AND THE
C           DERIVATIVE OF THE MODIFIED RESIDUAL THE
C           FIRST TIME, BUT ONLY THE DERIVATIVE IN
C           SUBSEQUENT ITERATIONS.
C
710   CALL DPA (N,KG,M,MMAX,NPP2,IV,T,Y,W,WEIGHT,ALF,IERR,IPRINT,A,C, 
     2          A(1,NP1),RNORM)
      GNSTEP=1.0
      ITERIN=0
      IF (ITER.GT.0) GO TO 720
      IF (KG.EQ.0) GO TO 860
      IF (IERR.NE.1) GO TO 880
C
      IF (IPRINT.LE.0) GO TO 720
      WRITE (OUTPUT,920) ITERIN,RNORM
      IF (IP) WRITE (10,920) ITERIN,RNORM
      WRITE (OUTPUT,890) NU
      IF (IP) WRITE (10,890) NU
C       BEGIN TWO-STAGE ORTHOGONAL FACTORIZATION
720   CALL ORFAC1 (KGP1,MMAX,M,N,IPRINT,A(1,B1),PRJRES,IERR)
      IF (IERR.LT.0) GO TO 880
      IERR=2
      IF (NU.EQ.0.) GO TO 740 
C
C       BEGIN INNER ITERATION LOOP FOR GENERATING 
C       NEW ALF AND TESTING IT FOR ACCEPTANCE.
C
730   CALL ORFAC2 (KGP1,MMAX,NU,A(1,B1))
C
C       SOLVE A KG X KG UPPER TRIANGULAR SYSTEM
C        FOR DELTA-ALF.  THE TRANSFORMED RESIDUAL 
C       (IN COL. NKG2 OF A) IS OVERWRITTEN BY
C       THE RESULT DELTA-ALF. 
C
740   CALL BACSUB (MMAX,KG,A(1,B1),A(1,NKG2))
      DO 750 K=1,KG 
750      A(K,B1)=ALF(K)+A(K,NKG2)
C
C       NEW ALF(K) = ALF(K) + DELTA-ALF(K)
C
C            STEP TO THE NEW POINT NEW ALF, AND COMPUTE
C            THE NEW NORM OF THE RESIDUAL.  NEW ALF
C            IS STORED IN COLUMN B1 OF A.
C
760   CALL DPA (N,KG,M,MMAX,NPP2,IV,T,Y,W,WEIGHT,A(1,B1),IERR,IPRINT,A,
     2          C,A(1,NP1),RNEW)
      IF (IERR.NE.2) GO TO 880
      ITER=ITER+1
      ITERIN=ITERIN+1
      SKIP=MOD(ITER,MODIT).NE.0
      IF (SKIP) GO TO 770
      WRITE (OUTPUT,900) ITER 
      IF (IP) WRITE (10,900) ITER
      WRITE (OUTPUT,940) (A(K,B1),K=1,KG)
      IF (IP) WRITE (10,940) (A(K,B1),K=1,KG)
      WRITE (OUTPUT,920) ITERIN,RNEW
      IF (IP) WRITE (10,920) ITERIN,RNEW
C
770   IF (ITER.LT.ITMAX) GO TO 780
      IERR=-1
      CALL VARERR (IPRINT,IERR,1)
      GO TO 870
780   IF (RNEW.LT.RNORM) GO TO 830
C
C             RETRACT THE STEP JUST TAKEN
C
      IF (NU.NE.0) GO TO 800
C             GAUSS-NEWTON OPTION ONLY
      GNSTEP=0.5*GNSTEP
      IF (GNSTEP.LT.EPS1) GO TO 870
      DO 790 K=1,KG 
790      A(K,B1)=ALF(K)+GNSTEP*A(K,NKG2)
      GO TO 760
C        ENLARGE THE MARQUARDT PARAMETER
800   NU=1.5*NU
      IF (.NOT.SKIP) WRITE (OUTPUT,910) NU
      IF ((.NOT.SKIP).AND.IP) WRITE (10,910) NU
      IF (NU.LE.100) GO TO 810
      IERR=-2
      CALL VARERR (IPRINT,IERR,1)
      GO TO 870
C        RETRIEVE UPPER TRIANGULAR FORM 
C        AND RESIDUAL OF FIRST STAGE.
810   DO 820 K=1,KG 
         KSUB=NP1+K 
         DO 820 J=K,KGP1
            JSUB=NP1+J
            ISUB=KGP1+J
820         A(K,JSUB)=A(ISUB,KSUB)
      GO TO 730
C
C       END OF INNER ITERATION LOOP
C
C     ACCEPT THE STEP JUST TAKEN
C
830   DELR=RNORM-RNEW
      RNORM=RNEW
      ALFNRM=0.
      DALFNM=0.
      DO 840 K=1,KG 
         ALF(K)=A(K,B1)
         ALFNRM=ALFNRM+ALF(K)**2
840      DALFNM=DALFNM+A(K,NKG2)**2
      ALFNRM=SQRT(ALFNRM)
      DALFNM=GNSTEP*SQRT(DALFNM)
      ACUM=DALFNM/ALFNRM
C
C        IF ITERIN IS GREATER THAN 1, A STEP WAS
C        RETRACTED DURING THIS OUTER ITERATION.
C
      IF (ITERIN.EQ.1) NU=0.5*NU
      IF (SKIP) GO TO 850
      WRITE (OUTPUT,890) NU
      IF (IP) WRITE (10,890) NU
      WRITE (OUTPUT,930) ACUM 
      IF (IP) WRITE (10,930) ACUM
850   IERR=3
      IF (PRJRES.GT.EPS1.OR.DELR.GT.(RNORM+1.D0)*EPS2) GO TO 710
C        END OF OUTER ITERATION LOOP
C
C        CALCULATE FINAL QUANTITIES -- LINEAR PARAMETERS,
C        RESIDUALS, COVARIANCE MATRIX, ETC.
C
860   IERR=ITER
870   IF (KG.GT.0) CALL DPA (N,KG,M,MMAX,NPP2,IV,T,Y,W,WEIGHT,ALF,4,
     2                       IPRINT,A,C,A(1,NP1),RNORM)
      CALL POSTPR (N,KG,M,NKG2,RNORM,IPRINT,ALF,W,A,A(1,NP1),C)
880   RETURN
C
890   FORMAT (1H ,4X,'NU =',E15.7)
900   FORMAT (' '/' '/3X,'ITERATION',I4,4X,'NONLINEAR ',10HPARAMETERS)
910   FORMAT (1H ,4X,'STEP RETRACTED, NU =',1PE15.7)
920   FORMAT (' '/' ',I5,29H  WEIGHTED NORM OF RESIDUAL =,1PE16.7)
930   FORMAT (1H ,4X,'NORM(DELTA-ALF) / NORM(ALF) =',1PE12.3)
940   FORMAT (' '/' ',1PE16.7,7X,1PE16.7)
C
      END 
C***********************************************************************
C
      SUBROUTINE ORFAC1 (KGP1,MMAX,M,N,IPRINT,B,PRJRES,IERR)
C
C***********************************************************************
C
C         STAGE 1:  HOUSEHOLDER REDUCTION OF
C
C                (    .    )      ( DR'. R3 )    KG
C                ( DR . R2 )  TO  (----. -- ),
C                (    .    )      (  0 . R4 )  M-N-KG
C
C                  KG    1          KG   1
C
C      WHERE DR = -D(Q2)*Y IS THE DERIVATIVE OF THE MODI-
C      FIED RESIDUAL PRODUCED BY DPA, R2 IS THE TRANS-
C      FORMED RESIDUAL FROM DPA, AND DR' IS IN UPPER
C      TRIANGULAR FORM (AS IN REF. (2), P. 18).
C      DR IS STORED IN ROWS N+1 TO M AND COLUMNS N+2 TO
C      N + KG + 1 OF THE MATRIX A (I.E.,COLUMNS 1 TO KG
C      OF THE MATRIX B).  R2 IS STORED IN COLUMN
C      N + KG + 2 OF THE MATRIX A (COLUMN KG + 1 OF B).
C      FOR K = 1, 2, ..., KG, FIND REFLECTION
C      I - U * U' / BETA WHICH ZEROES B(I, K),
C      I = N+K+1, ..., M.
C
C  ......................................................
C
      REAL B(MMAX,KGP1)
C
      KG=KGP1-1
      KG23=2*KG+3
      NP1=N+1
      NPKG=N+KG
C
      DO 990 K=1,KG 
         ALPHA=0.0
         NPK=N+K
         DO 950 I=NPK,M
950         ALPHA=ALPHA+B(I,K)**2
         ALPHA=SIGN(SQRT(ALPHA),B(NPK,K))
         U1=B(NPK,K)+ALPHA
         B(NPK,K)=U1
         BETA=ALPHA*U1
         IF (ALPHA.NE.0.0) GO TO 960
C                          COLUMN WAS ZERO
         IERR=-8
         CALL VARERR (IPRINT,IERR,NP1+K)
         GO TO 1040 
C                    APPLY REFLECTIONS TO REMAINING COLUMNS 
C                    OF B AND TO RESIDUAL VECTOR. 
960      KP1=K+1
         DO 980 J=KP1,KGP1
            ACUM=0.0
            DO 970 I=NPK,M
970            ACUM=ACUM+B(I,K)*B(I,J)
            ACUM=ACUM/BETA
            DO 980 I=NPK,M
980            B(I,J)=B(I,J)-B(I,K)*ACUM
990      B(NPK,K)=-ALPHA
C
      PRJRES=0.
      DO 1000 J=NP1,NPKG
1000     PRJRES=PRJRES+B(J,KGP1)**2
      PRJRES=SQRT(PRJRES)
C
C        SAVE UPPER TRIANGULAR FORM AND TRANSFORMED
C        RESIDUAL, FOR USE IN CASE A STEP IS RETRACTED.
C        ALSO COMPUTE COLUMN LENGTHS.
C
      IF (IERR.EQ.4) GO TO 1040
      DO 1030 K=1,KG
         NPK=N+K
         ACUM=0.
         DO 1010 J=NP1,NPK
1010        ACUM=ACUM+B(J,K)**2
         DO 1020 J=K,KGP1
            JSUB=KGP1+J
            B(K,J)=B(NPK,J)
1020        B(JSUB,K)=B(NPK,J)
         B(KG23,K)=SQRT(ACUM) 
1030  CONTINUE
C
1040  RETURN
C
      END 
C***********************************************************************
C
      SUBROUTINE ORFAC2 (KGP1,MMAX,NU,B)
C
C***********************************************************************
C
C     STAGE 2:  SPECIAL HOUSEHOLDER REDUCTION OF
C
C         KG     ( DR' . R3 )      (DR'' . R5 )
C                (-----. -- )      (-----. -- )
C     M-N-KG     (  0  . R4 )  TO  (  0  . R4 )
C                (-----. -- )      (-----. -- )
C         KG     (NU*D . 0  )      (  0  . R6 )
C
C                   KG    1          KG    1
C
C     WHERE DR', R3, AND R4 ARE AS IN ORFAC1, NU IS THE
C     MARQUARDT PARAMETER, D IS A DIAGONAL MATRIX 
C     CONSISTING OF THE LENGTHS OF THE COLUMNS OF DR',
C     AND DR'' IS IN UPPER TRIANGULAR FORM.  DETAILS IN
C     (1), PP. 423-424.  NOTE THAT THE (M-N-KG) BAND OF
C     ZEROES, AND R4, ARE OMITTED IN STORAGE.
C
C  ......................................................
C
      REAL B(MMAX,KGP1),NU
C
      KG=KGP1-1
      KG2=2*KG
      KG23=KG2+3
      DO 1070 K=1,KG
         KP1=K+1
         KGPK=KG+K
         KGPKM1=KGPK-1
         B(KGPK,K)=NU*B(KG23,K)
         ALPHA=B(K,K)**2
         DO 1050 I=KGP1,KGPK
1050        ALPHA=ALPHA+B(I,K)**2
         ALPHA=SIGN(SQRT(ALPHA),B(K,K)) 
         U1=B(K,K)+ALPHA
         BETA=ALPHA*U1
         B(K,K)=-ALPHA
C              THE K-TH REFLECTION MODIFIES ONLY ROWS K,
C              KG+1, KG+2, ..., KG+K, AND COLUMNS K TO KG+1.
         DO 1070 J=KP1,KGP1
            B(KGPK,J)=0.
            ACUM=U1*B(K,J)
            DO 1060 I=KGP1,KGPKM1
1060           ACUM=ACUM+B(I,K)*B(I,J)
            ACUM=ACUM/BETA
            B(K,J)=B(K,J)-U1*ACUM
            DO 1070 I=KGP1,KGPK
1070           B(I,J)=B(I,J)-B(I,K)*ACUM
C
      RETURN
C
      END 
C***********************************************************************
C
      SUBROUTINE DPA (N,KG,M,MMAX,NPP2,IV,T,Y,W,WEIGHT,ALF,ISEL,IPRINT,A
     2,U,R,RNORM)
C
C***********************************************************************
C
C     COMPUTE THE NORM OF THE RESIDUAL (IF ISEL = 1 OR 2),
C     OR THE (M-N) X KG DERIVATIVE OF THE MODIFIED
C     RESIDUAL (M-N) VECTOR Q2*Y (IF ISEL = 1 OR 3).
C     HERE Q * PHI = S, I.E., 
C
C    N     ( Q1 ) (     .   .        )   ( S  . R1 .  F1  ) 
C          (----) ( PHI . Y . D(PHI) ) = (--- . -- . ---- ) 
C    M-N   ( Q2 ) (     .   .        )   ( 0  . R2 .  F2  ) 
C
C            M       N    1      P         N     1     P
C
C     WHERE Q IS M X M ORTHOGONAL, AND S IS N X N UPPER
C     TRIANGULAR.  THE NORM OF THE RESIDUAL = NORM(R2),
C     AND THE DESIRED DERIVATIVE ACCORDING TO REF.(5), IS
C
C                                         -1
C              D(Q2 * Y) = -Q2 * D(PHI)* S  * Q1* Y.
C
C  ......................................................
C
      PARAMETER (MAXN1=11, MAXKG=8)
      COMMON /BLKINC/ INC
      REAL A(MMAX,NPP2),ALF(KG),T(MMAX,IV),W(M),WEIGHT(M),Y(M),R(M),U(N)
      INTEGER FIRSTC,FIRSTR,P,INC(MAXKG,MAXN1)
      LOGICAL PHINP1
C
      NP1=N+1
      CALL ADA (NP1,N,KG,M,MMAX,NPP2,IV,A,T,ALF,MIN0(ISEL,3),PHINP1,
     2          NCONP1)
      GO TO (1100,1090,1080,1080), ISEL 
C                                      ISEL = 3 OR 4
1080  FIRSTC=NP2
      LASTC=NPP1
      FIRSTR=(4-ISEL)*N+1
      GO TO 1200
C                                      ISEL = 2
1090  FIRSTC=NCONP1 
      LASTC=NP1
      IF (NCON.EQ.0) GO TO 1150
      IF (A(1,NCON).EQ.SAVE) GO TO 1150 
      ISEL=-7
      CALL VARERR (IPRINT,ISEL,NCON)
      GO TO 1320
C                             ISEL = 1 (FIRST TIME ONLY)
1100  NP2=N+2
      NKG2=N+2+KG
      NPP1=NPP2-1
      FIRSTC=1
      LASTC=NPP1
      FIRSTR=NP1
      DO 1110 I=1,M 
         IF (WEIGHT(I).GE.0.) GO TO 1110
C                                      ERROR IN WEIGHTS
         ISEL=-6
         II=I
         CALL VARERR (IPRINT,ISEL,II)
         GO TO 1320 
1110     W(I)=SQRT(WEIGHT(I)) 
C
      NCON=N
      NCONP1=NP1
      IF (KG.EQ.0) GO TO 1160 
C
      IF (N.NE.0) GO TO 1120
      PHINP1=.TRUE. 
      GO TO 1180
C                      DETERMINE NUMBER OF CONSTANT FCNS.
1120  P=0 
      DO 1130 J=1,NP1
         IF (P.EQ.0) NCONP1=J 
         DO 1130 K=1,KG
            INCKJ=INC(K,J)
            IF (INCKJ.EQ.1) P=P+1
1130  CONTINUE
C
      NCON=NCONP1-1 
C                DETERMINE IF PHI(N+1) IS IN THE MODEL.
      PHINP1=.FALSE.
      DO 1140 K=1,KG
         IF (INC(K,NP1).EQ.1) PHINP1=.TRUE.
1140  CONTINUE
C                                     ISEL = 1 OR 2
1150  IF (PHINP1) GO TO 1180
1160  DO 1170 I=1,M 
1170     R(I)=Y(I)
      GO TO 1200
1180  DO 1190 I=1,M 
1190     R(I)=Y(I)-R(I)
C                             WEIGHT APPROPRIATE COLUMNS
1200  DO 1210 I=1,M 
         ACUM=W(I)
         DO 1210 J=FIRSTC,LASTC
1210        A(I,J)=A(I,J)*ACUM
C
C        COMPUTE ORTHOGONAL FACTORIZATIONS BY HOUSEHOLDER
C        REFLECTIONS.  IF ISEL = 1 OR 2, REDUCE PHI
C        (STORED IN THE FIRST N COLUMNS OF THE MATRIX A)
C        TO UPPER TRIANGULAR FORM, (Q*PHI = S), AND
C        TRANSFORM Y (STORED IN COLUMN N+1), GETTING
C        Q*Y = R.  IF ISEL = 1, ALSO TRANSFORM J = D PHI
C        (STORED IN COLUMNS N+2 THROUGH N+P+1 OF THE
C        MATRIX A), GETTING Q*J = F.  IF ISEL = 3 OR 4,
C        PHI HAS ALREADY BEEN REDUCED, TRANSFORM ONLY J.
C        S, R, AND F OVERWRITE PHI, Y, AND J, RESPEC-
C        TIVELY AND A FACTORED FORM OF Q IS SAVED IN U
C        AND THE LOWER TRIANGLE OF PHI. 
C
      IF (N.EQ.0) GO TO 1260
      DO 1250 K=1,N 
         KP1=K+1
         IF (ISEL.GE.3.OR.(ISEL.EQ.2.AND.K.LT.NCONP1)) GO TO 1230
         ALPHA=0.
         DO 1220 I=K,M
1220        ALPHA=ALPHA+A(I,K)**2
         ALPHA=SIGN(SQRT(ALPHA),A(K,K)) 
         U(K)=A(K,K)+ALPHA
         A(K,K)=-ALPHA
         FIRSTC=KP1 
         IF (ALPHA.NE.0.0) GO TO 1230
         ISEL=-8
         KK=K
         CALL VARERR (IPRINT,ISEL,KK)
         GO TO 1320 
C                              APPLY REFLECTIONS TO COLUMNS 
C                              FIRSTC TO LASTC.
1230     BETA=-A(K,K)*U(K)
         DO 1250 J=FIRSTC,LASTC
            ACUM=U(K)*A(K,J)
            DO 1240 I=KP1,M
1240           ACUM=ACUM+A(I,K)*A(I,J)
            ACUM=ACUM/BETA
            A(K,J)=A(K,J)-U(K)*ACUM
            DO 1250 I=KP1,M
1250           A(I,J)=A(I,J)-A(I,K)*ACUM
C
1260  IF (ISEL.GE.3) GO TO 1280
      RNORM=0.
      DO 1270 I=NP1,M
1270     RNORM=RNORM+R(I)**2
      RNORM=SQRT(RNORM)
      IF (ISEL.EQ.2) GO TO 1320
      IF (NCON.GT.0) SAVE=A(1,NCON)
C
C        F2 IS NOW CONTAINED IN ROWS N+1 TO M AND COLUMNS
C        N+2 TO N+P+1 OF THE MATRIX A.  NOW SOLVE THE
C        N X N UPPER TRIANGULAR SYSTEM S*C = R1 FOR THE
C        LINEAR PARAMETERS C.  C OVERWRITES R1.
C
1280  IF (N.GT.0) CALL BACSUB (MMAX,N,A,R)
C
C        MAJOR PART OF KAUFMAN'S SIMPLIFICATION OCCURS
C        HERE.  COMPUTE THE DERIVATIVE OF ETA WITH
C        RESPECT TO THE NONLINEAR PARAMETERS
C
C    T  D ETA      T   N       D PHI(J)   D PHI(N+1)
C   Q *-------- = Q *(SUM C(J) -------- + ----------) =  F2*C
C      D ALF(K)        J=1      D ALF(K)   D ALF(K)
C
C        AND STORE THE RESULT IN COLUMNS N+2 TO N+KG+1.
C        IF ISEL NOT = 4, THE FIRST N ROWS ARE OMITTED.
C        THIS IS -D(Q2)*Y.  IF ISEL NOT = 4 THE RESIDUAL
C        R2 = Q2*Y (IN COL. N+1) IS COPIED TO COLUMN
C        N+KG+2.  OTHERWISE ALL OF COLUMN N+1 IS COPIED.
C
      DO 1310 I=FIRSTR,M
         IF (N.EQ.NCON) GO TO 1310
         L=NP1
         DO 1300 K=1,KG
            ACUM=0. 
            DO 1290 J=NCONP1,N
               IF (INC(K,J).EQ.0) GO TO 1290
               L=L+1
               ACUM=ACUM+A(I,L)*R(J)
1290        CONTINUE
            KSUB=NP1+K
            IF (INC(K,NP1).EQ.0) GO TO 1300
            L=L+1
            ACUM=ACUM+A(I,L)
1300        A(I,KSUB)=ACUM
1310     A(I,NKG2)=R(I)
C
1320  RETURN
C
      END 
C***********************************************************************
C
      SUBROUTINE BACSUB (MMAX,N1,A1,Z)
C
C***********************************************************************
C
C     BACKSOLVE THE N1 X N1 UPPER TRIANGULAR SYSTEM
C     A1*Z = B.  THE SOLUTION Z OVERWRITES THE RIGHT
C     SIDE B.
C
      REAL A1(MMAX,N1),Z(N1)
C
      Z(N1)=Z(N1)/A1(N1,N1)
      IF (N1.EQ.1) GO TO 1350 
      NP1=N1+1
      DO 1340 IBACK=2,N1
         I=NP1-IBACK
C           I = N-1, N-2, ..., 2, 1
         IP1=I+1
         ACUM=Z(I)
         DO 1330 J=IP1,N1
1330        ACUM=ACUM-A1(I,J)*Z(J)
1340     Z(I)=ACUM/A1(I,I)
C
1350  RETURN
C
      END 
C***********************************************************************
C
      SUBROUTINE POSTPR (N,KG,M,NKG2,RNORM,IPRINT,ALF,W,A,R,U)
C
C***********************************************************************
C
C     CALCULATE RESIDUALS, SAMPLE VARIANCE, AND
C     COVARIANCE MATRIX.  ON INPUT, U CONTAINS
C     INFORMATION ABOUT HOUSEHOLDER REFLECTIONS FROM DPA.
C     ON OUTPUT, IT CONTAINS THE LINEAR PARAMETERS.
C
      PARAMETER (MMAX=300, MAXNKG=18)
      REAL A(MMAX,NKG2),ALF(KG),R(M),U(N),W(M)
      INTEGER OUTPUT
      LOGICAL IP
      REAL BM(MMAX,MAXNKG)
      COMMON /BLKBM/ BM
      COMMON /BLOKIP/ IP
      DATA OUTPUT /6/
C
      NP1=N+1
      NPKG=NKG2-2
      NKG1=NPKG+1
C
C           UNWIND HOUSEHOLDER TRANSFORMATIONS TO GET
C           RESIDUALS AND MOVE THE LINEAR PARAMETERS
C           FROM R TO U.
C
      IF (N.EQ.0) GO TO 1390
      DO 1370 KBACK=1,N
         K=NP1-KBACK
         KP1=K+1
         ACUM=0.
         DO 1360 I=KP1,M
1360        ACUM=ACUM+A(I,K)*R(I)
         SAVE1=R(K) 
         R(K)=ACUM/A(K,K)
         ACUM=-ACUM/(U(K)*A(K,K))
         U(K)=SAVE1 
         DO 1370 I=KP1,M
1370        R(I)=R(I)-A(I,K)*ACUM
      IF (KG.EQ.0) GO TO 1410 
      DO 1380 J=NP1,NPKG
         JP1=J+1
         DO 1380 I=1,M
1380        BM(I,J)=A(I,JP1)/W(I)
C
C           THE FIRST N COLUMNS OF THE MATRIX HAVE BEEN
C           REDUCED TO UPPER TRIANGULAR FORM IN DPA.
C           FINISH BY REDUCING ROWS N+1 TO M AND COLUMNS
C           N+2 THROUGH N+KG+1 TO TRIANGULAR FORM.  THEN
C           SHIFT COLUMNS OF DERIVATIVE MATRIX OVER ONE
C           TO THE LEFT TO BE ADJACENT TO THE FIRST N
C           COLUMNS.
C
1390  CALL ORFAC1 (KG+1,MMAX,M,N,IPRINT,A(1,N+2),PRJRES,4)
      DO 1400 I=1,M 
         A(I,NKG2)=R(I)
         DO 1400 K=NP1,NKG1
1400        A(I,K)=A(I,K+1)
C                           COMPUTE COVARIANCE MATRIX
1410  A(1,NKG2)=RNORM
      ACUM=RNORM*RNORM/(M-N-KG)
      A(2,NKG2)=ACUM
      CALL COV (MMAX,NPKG,ACUM,A)
C
      IF (IPRINT.LT.0) GO TO 1420
      WRITE (OUTPUT,1430)
      IF (IP) WRITE (10,1430) 
      IF (N.GT.0) WRITE (OUTPUT,1440) (U(J),J=1,N)
      IF ((N.GT.0).AND.(IP)) WRITE (10,1440) (U(J),J=1,N)
      IF (KG.GT.0) WRITE (OUTPUT,1450) (ALF(K),K=1,KG)
      IF ((KG.GT.0).AND.(IP)) WRITE (10,1450) (ALF(K),K=1,KG)
      WRITE (OUTPUT,1460) RNORM,ACUM
      IF (IP) WRITE (10,1460) RNORM,ACUM
      WRITE (OUTPUT,1430)
      IF (IP) WRITE (10,1430) 
1420  RETURN
C
1430  FORMAT (1H /1H ,59(1H'))
1440  FORMAT (1H /20H   LINEAR PARAMETERS//(4E15.7))
1450  FORMAT (1H /23H   NONLINEAR PARAMETERS//(4E15.7))
1460  FORMAT (1H /30H   WEIGHTED NORM OF RESIDUAL =,E15.7/1H /
     2        32H   WEIGHTED ESTIMATED VARIANCE =,E15.7)
C
      END 
C***********************************************************************
C
      SUBROUTINE COV (MMAX,NPKG,SIGMA2,A)
C
C***********************************************************************
C
C        COMPUTE THE UNSCALED COVARIANCE MATRIX OF THE
C     N + KG PARAMETERS.  THIS INVOLVES COMPUTING 
C
C                       2     -1    -T
C                  SIGMA  *  T   * T
C
C     WHERE THE (N+KG) X (N+KG) UPPER TRIANGULAR MATRIX T
C     IS DESCRIBED IN SUBROUTINE POSTPR.  THE RESULT
C     OVERWRITES THE FIRST N+KG ROWS AND COLUMNS OF THE
C     MATRIX A.  THE RESULTING MATRIX IS SYMMETRIC.
C
C        REFERENCE:  LAWSON AND HANSON, SOLVING LEAST
C     SQUARES PROBLEMS, PRENTICE-HALL, 1974, PP. 67-
C     70, 281.
C
C  ......................................................
C
      REAL A(MMAX,NPKG)
C
      DO 1470 J=1,NPKG
1470     A(J,J)=1./A(J,J)
C
C              INVERT T UPON ITSELF
C
      IF (NPKG.EQ.1) GO TO 1500
      NM1=NPKG-1
      DO 1490 I=1,NM1
         IP1=I+1
         DO 1490 J=IP1,NPKG
            JM1=J-1 
            SUM=0.
            DO 1480 L=I,JM1
1480           SUM=SUM+A(I,L)*A(L,J)
1490        A(I,J)=-SUM*A(J,J)
C
C              NOW FORM THE MATRIX PRODUCT
C
1500  DO 1520 I=1,NPKG
         DO 1520 J=I,NPKG
            SUM=0.
            DO 1510 L=J,NPKG
1510           SUM=SUM+A(I,L)*A(J,L)
            SUM=SUM*SIGMA2
            A(I,J)=SUM
1520        A(J,I)=SUM
C
      RETURN
C
      END 
C***********************************************************************
C
      SUBROUTINE VARERR (IPRINT,IERR,K) 
C
C***********************************************************************
C
C     PRINT ERROR MESSAGES
C
      LOGICAL IP
      COMMON /BLOKIP/ IP
      INTEGER ERRNO,OUTPUT
      DATA OUTPUT /6/
C
      IF (IPRINT.LT.0) GO TO 1600
      ERRNO=IABS(IERR)
      GO TO (1530,1540,1600,1550,1560,  1570,1580,1590), ERRNO
C
1530  WRITE (OUTPUT,1610)
      IF (IP) WRITE (10,1610) 
      GO TO 1600
1540  WRITE (OUTPUT,1620)
      IF (IP) WRITE (10,1620) 
      GO TO 1600
1550  WRITE (OUTPUT,1630)
      IF (IP) WRITE (10,1630) 
      GO TO 1600
1560  WRITE (OUTPUT,1640)
      IF (IP) WRITE (10,1640) 
      GO TO 1600
1570  WRITE (OUTPUT,1650) K
      IF (IP) WRITE (10,1650) K
      GO TO 1600
1580  WRITE (OUTPUT,1660) K
      IF (IP) WRITE (10,1660) K
      GO TO 1600
1590  WRITE (OUTPUT,1670) K
      IF (IP) WRITE (10,1670) K
C
1600  RETURN
1610  FORMAT (' '/2X,'PROBLEM TERMINATED FOR EXCESSIVE ITERATIONS'//) 
1620  FORMAT (' '/2X,'PROBLEM TERMINATED BECAUSE OF ILL-CONDITIONING'//)
1630  FORMAT (' '/2X,'INPUT ERROR IN PARAMETER N,KG,M,NPP2, OR MMAX.'//)
1640  FORMAT (' '/2X,'ERROR -- INC MATRIX IMPROPERLY SPECIFIED',
     2        ' OR DISAGREES WITH NPP2.'//)
1650  FORMAT (' '/2X,'ERROR -- WEIGHT(',I4,') NOT NONNEGATIVE.'//)
1660  FORMAT (' '/2X,'ERROR -- CONSTANT COLUMN',I3,
     2        ' MUST BE COMPUTED ONLY WHEN ISEL = 1.'//)
1670  FORMAT (' '/2X,'CATASTROPHIC FAILURE -- COLUMN',I4,
     2        ' IS ZERO. SEE RALL AND FUNDERLIC''S DOCUMENTATION.'//) 
C
      END 
C***********************************************************************
C
      SUBROUTINE ADA (NP1,N,KG,M,MMAX,NPP2,IV,A,T,ALF,ISEL,PHINP1,
     2 NCONP1)
C
C***********************************************************************
C
C  EVALUATE THE FUNCTIONS PHI AND THEIR PARTIAL DERIVATIVES 
C  D PHI(J)/D ALF(K) AT THE SAMPLE POINTS T(I).
C    ISEL=1 MEANS COMPUTE BOTH THE FUNCTIONS AND
C           DERIVATIVES
C    ISEL=2 MEANS COMPUTE ONLY THE FUNCTIONS PHI
C    ISEL=3 MEANS COMPUTE ONLY THE DERIVATIVES
C
      PARAMETER (MAXN1=11, MAXN=10, MAXKG=8, MAXP=21)
      REAL A(MMAX,NPP2),ALF(KG),T(MMAX,IV)
      REAL C(MAXN),NDI(MAXP),DPH(MAXN1,MAXKG)
      INTEGER INC(MAXKG,MAXN1),P
      LOGICAL PHINP1,NC
      COMMON /NCBLK/ NC
      COMMON /BLKINC/ INC
      COMMON /BLKPH/ PH
      COMMON /BLKNDI/ NDI
      COMMON /BLKDPH/ DPH
      COMMON /BLOCKP/ P
      REAL PH(MAXN1)
C
      GO TO (1680,1690,1710), ISEL
C
C  CALCULATE FUNCTION VALUES & STORE THEM IN FIRST (N+1)
C  COLUMNS OF A.
C
1680  PHINP1=.TRUE. 
      NPH=NP1
      NCONP1=1
1690  IF (.NOT.PHINP1) NPH=N
      DO 1700 I=1,M 
         II=I
         CALL CALC (II,KG,N,MMAX,IV,C,ALF,T,ETA)
         DO 1700 J=NCONP1,NPH 
            A(I,J)=PH(J)
1700  CONTINUE
      NC=.TRUE.
      IF (ISEL.EQ.2) GO TO 1730
      IF (KG.EQ.0) GO TO 1730 
C
C  CALCULATE NONVANISHING PARTIAL DERIVATIVES & STORE THEM
C  IN COLUMNS (N+2) THROUGH (N+P+2) OF A.
C
C  CHECK TO SEE IF ANY DERIVATIVES ARE TO BE CALCULATED NUMERICALLY
C
1710  DO 1712 I=1,P 
         IF (NDI(I).NE.0.0) GOTO 1715
1712  CONTINUE
C
C  ALL DERIVATIVES ARE TO BE CALCULATED FROM CLOSED FORM
C
      DO 1714 K=1,M 
         KK = K
         CALL DERIV (KK,T,MMAX,IV,ALF,KG,NP1)
         IK = NP1
         DO 1714 I=1,KG
            DO 1714 J=1,NP1
               IF (INC(I,J).NE.1) GOTO 1714
               IK = IK+1
               A(K,IK) = DPH(J,I)
1714  CONTINUE
      RETURN
1715  CONTINUE
C
C  ONE OR MORE DERIVATIVES ARE TO BE CALCULATED NUMERICALLY 
C
      IK=NP1
      DO 1720 I=1,KG
         II=I
         DO 1720 J=1,NP1
            JJ=J
            IF (INC(I,J).NE.1) GO TO 1720
            IK=IK+1 
            IN=IK-NP1
            DO 1719 K=1,M
               KK=K 
               CALL NUMDER (II,JJ,KK,IN,N,KG,IV,MMAX,T,ALF,A(K,IK))
1719        CONTINUE
1720  CONTINUE
1730  RETURN
C
      END 
C***********************************************************************
C
      SUBROUTINE NUMDER (I,J,K,IN,N,KG,IV,MMAX,T,ALF,DV)
C
C***********************************************************************
C
C  GIVEN AN ARRAY OF NUMERICAL DERIVATIVE INCREMENTS,
C  CALCULATE DERIVATES NUMERICALLY IF INCREMENT IS NOT ZERO,
C  USE CLOSED FORM IN SUBROUTINE DERIV IF INCREMENT IS ZERO.
C
      PARAMETER (MAXN=10, MAXN1=11, MAXKG=8, MAXP=21)
      REAL T(MMAX,IV),ALF(KG),NDI(MAXP) 
      LOGICAL NC
      COMMON /NCBLK/ NC
      COMMON /BLKNDI/ NDI
      COMMON /BLKDPH/ DPH
      REAL DPH(MAXN1,MAXKG),CC(MAXN)
C
      DK=NDI(IN)
      IF (DK.EQ.0.) GO TO 1750
C
C  DERIVATIVE TO BE CALCULATED NUMERICALLY USING FORWARD
C  DIFFERENCE QUOTIENT.
C
      NP1=N+1
      DO 1740 JJ=1,N
1740     CC(JJ)=0.
      IF (J.NE.NP1) CC(J)=1.
      NC=.FALSE.
      CALL CALC (K,KG,N,MMAX,IV,CC,ALF,T,YC)
      ASAVE=ALF(I)
      ALF(I)=ASAVE+DK
      CALL CALC (K,KG,N,MMAX,IV,CC,ALF,T,YD)
      NC=.TRUE.
      DV=(YD-YC)/DK 
      ALF(I)=ASAVE
      RETURN
C
C  DERIVATIVE TO BE CALCULATED FROM CLOSED FORM IN
C  SUBROUTINE DERIV.
C
1750  CALL DERIV (K,T,MMAX,IV,ALF,KG,NP1)
      DV=DPH(J,I)
      RETURN
C
      END 
