C--------------------- DOUBLE PRECISION VERSION 2B ---------------------MAI00100
C  DISCRETE - VERSION 2B (December 1990).  FOR THE AUTOMATIC ANALYSIS OFMAI00200
C      MULTICOMPONENT EXPONENTIAL DECAY DATA FOR UP TO 9 COMPONENTS.    MAI00300
C  MUST BE RUN ON A COMPUTER WITH AT LEAST 5 SIGNIFICANT FIGURES IN     MAI00400
C      ITS SINGLE PRECISION FORTRAN ARITHMETIC AND 9 SIG. FIGS. IN ITS  MAI00500
C      DOUBLE PRECISION PRECISION FORTRAN ARITHMETIC.                   MAI00600
C  REFERENCES - S. W. PROVENCHER (1976) WRITE UP OF DISCRETE,           MAI00700
C                                (1980) Math. Biosci., Vol. 50, 251-262,MAI00800
C                                (1976) J. CHEM. PHYS., VOL. 64, 2772-7,MAI00900
C                                (1976) BIOPHYS. J., VOL. 16, 27-41.    MAI01000
C-----------------------------------------------------------------------MAI01100
C  CALLS SUBPROGRAMS - BLOCK DATA, DATAIO, WEIGHT, FANLYZ, YANLYZ       MAI01200
C  WHICH IN TURN CALL - LSTSQR, EVAR, VARF, ETHEOR, PIVOT, PIVOT1,      MAI01300
C                       ANLERR, FISHER, RESIDU, PLPRES                  MAI01400
C-----------------------------------------------------------------------MAI01500
      LOGICAL LAST, REGINT, LBLOKA                                      MAI01600
      DOUBLE PRECISION GSE, YLYFIT, DELTAT, E, TSTART, DBLOKA, T        MAI01700
      COMMON /BLOK/ DBLOKA(5933), RBLOKA(471), IBLOKA(301), NMAX,       MAI01800
     1 NINTMX, IBLOKB(12), REGINT, LBLOKA(45)                           MAI01900
C                                                                       MAI02000
C***********************************************************************MAI02100
C***********************************************************************MAI02200
C                                                                       MAI02300
C  THE FOLLOWING INSTRUCTIONS DESCRIBE ALL POSSIBLE NECESSARY CHANGES   MAI02400
C  THAT YOU MAY HAVE TO MAKE IN THE MAIN PROGRAM. (SEE ALSO THE         MAI02500
C  BLOCK DATA SUBPROGRAM.)                                              MAI02600
C                                                                       MAI02700
C***********************************************************************MAI02800
C  IF YOU CHANGE NMAX IN THE BLOCK DATA SUBPROGRAM, YOU MUST READJUST   MAI02900
C  THE DIMENSIONS IN THE FOLLOWING STATEMENT TO NMAX -                  MAI03000
C                                                                       MAI03100
      DIMENSION T(200), Y(200), SQRTW(200), YLYFIT(200)                 MAI03200
C                                                                       MAI03300
C***********************************************************************MAI03400
C  ORDINARILY E SHOULD BE DIMENSIONED E(NMAX,10) IN THE FOLLOWING       MAI03500
C  STATEMENT.  HOWEVER, IF YOU WILL ALWAYS INPUT REGINT=.TRUE., AND     MAI03600
C  IF NMAX IS GREATER THAN 200, YOU CAN SAVE STORAGE BY DIMENSIONING    MAI03700
C  IT AS E(1,10) -                                                      MAI03800
C                                                                       MAI03900
      DIMENSION E(200,10)                                               MAI04000
C                                                                       MAI04100
C***********************************************************************MAI04200
C  ORDINARILY GSE SHOULD BE DIMENSIONED GSE(500,4).  HOWEVER, IF YOU    MAI04300
C  WILL ALWAYS INPUT REGINT=.FALSE., AND IF NMAX IS LESS THAN 200,      MAI04400
C  YOU CAN SAVE A LITTLE STORAGE BY DIMENSIONING IT AS GSE(1,4) IN      MAI04500
C  THE FOLLOWING STATEMENT -                                            MAI04600
C                                                                       MAI04700
      DIMENSION GSE(500,4)                                              MAI04800
C                                                                       MAI04900
C***********************************************************************MAI05000
C  IF YOU CHANGE NINTMX IN THE BLOCK DATA SUBPROGRAM, YOU MUST          MAI05100
C  READJUST THE DIMENSIONS IN THE FOLLOWING STATEMENT TO NINTMX -       MAI05200
C                                                                       MAI05300
      DIMENSION TSTART(10), NT(10), DELTAT(10)                          MAI05400
C                                                                       MAI05500
C  THIS IS THE END OF ALL POSSIBLE CHANGES YOU MAY HAVE TO MAKE IN      MAI05600
C  THE MAIN PROGRAM.                                                    MAI05700
C                                                                       MAI05800
C***********************************************************************MAI05900
C***********************************************************************MAI06000
C                                                                       MAI06100
      EQUIVALENCE (E(1,1),GSE(1,1))                                     MAI06200
C-----------------------------------------------------------------------MAI06300
C  DATAIO IS FOR INPUT AND OUTPUT OF INPUT DATA                         MAI06400
  100 CALL DATAIO (T,Y,SQRTW,NMAX,TSTART,DELTAT,NT,NINTMX,LAST)         MAI06500
      IF (REGINT) GO TO 110                                             MAI06600
      NDIME=NMAX                                                        MAI06700
      NDIMG=1                                                           MAI06800
      GO TO 200                                                         MAI06900
  110 NDIME=1                                                           MAI07000
      NDIMG=500                                                         MAI07100
C-----------------------------------------------------------------------MAI07200
C  WEIGHT GENERATES CRUDE STARTING VALUES FOR CONSTRAINED LEAST SQUARES MAI07300
C      FITS TO RAW DATA TO GENERATE WEIGHTS AND DOES INITIAL            MAI07400
C      COMPUTATIONS OF QUANTITIES FOR LATER USE.                        MAI07500
  200 CALL WEIGHT (.FALSE.,Y,T,SQRTW,YLYFIT,E,NMAX,NDIME,               MAI07600
     1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG)                               MAI07700
C-----------------------------------------------------------------------MAI07800
C  FANLYZ IS FOR THE CONSTRAINED STEPWISE LEAST SQUARES ANALYSIS OF     MAI07900
C      THE RAW DATA AND THE TRANSFORMS TO GET STARTING ESTIMATES        MAI08000
C      FOR THE FINAL LEAST SQUARES ANALYSES OF THE RAW DATA.            MAI08100
      CALL FANLYZ (T,Y,SQRTW,YLYFIT,E,NMAX,NDIME,                       MAI08200
     1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG)                               MAI08300
C-----------------------------------------------------------------------MAI08400
C  WEIGHT USES STARTING VALUES FROM FANLYZ FOR CONSTRAINED LEAST        MAI08500
C      SQUARES FITS TO RAW DATA TO GENERATE WEIGHTS FOR THE FINAL       MAI08600
C      LEAST SQUARES FITS TO THE RAW DATA.                              MAI08700
      CALL WEIGHT (.TRUE.,Y,T,SQRTW,YLYFIT,E,NMAX,NDIME,                MAI08800
     1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG)                               MAI08900
C-----------------------------------------------------------------------MAI09000
C  YANLYZ IS FOR THE CONSTRAINED STEPWISE LEAST SQUARES ANALYSIS OF     MAI09100
C      THE RAW DATA (USING THE STARTING VALUES OBTAINED FROM FANLYZ     MAI09200
C      AND THE WEIGHTS FROM WEIGHT).                                    MAI09300
      CALL YANLYZ (T,Y,SQRTW,YLYFIT,E,NMAX,NDIME,                       MAI09400
     1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG)                               MAI09500
C-----------------------------------------------------------------------MAI09600
      IF (.NOT.LAST) GO TO 100                                          MAI09700
      STOP                                                              MAI09800
      END                                                               MAI09900
C--------------------- DOUBLE PRECISION VERSION 2B ---------------------BLK00100
      BLOCK DATA                                                        BLK00200
      LOGICAL REGINT, FAILED, PRPREL, NOBASE, WTED, REPEAT, PLOTRS,     BLK00300
     1 NONNEG, FFAIL, PIVALP, PRFINL                                    BLK00400
      REAL LAMST, LAMF, LAMNMX                                          BLK00500
      DOUBLE PRECISION TRBASE, GF, DGFL, DGRIDF, GFSTL, FHAT, VARMIN,   BLK00600
     1 SUMWTY, SUMWT, ZERO, DGRIDE, GESTL, VAR, PALPHA, PLAM,           BLK00700
     2 SE, PLMTRY, DELTAP, DFCAPL, SFDE, ADUB, DDSE, DTOLER, F, DF      BLK00800
      COMMON /BLOK/ TRBASE(19), GF(19,126), DGFL(19,126), DGRIDF,       BLK00900
     1 GFSTL, FHAT(19), VARMIN, SUMWTY, SUMWT, ZERO(10),                BLK01000
     2 DGRIDE, GESTL, VAR, PALPHA(10), PLAM(10), SE(10,10),             BLK01100
     3 PLMTRY(10), DELTAP(19), DFCAPL(10,9), SFDE(9), ADUB(19,19),      BLK01200
     4 DDSE(9,9), DTOLER(19), F(19,10), DF(19,10), DELTA, FHATMX,       BLK01300
     5 CONVRG, LAMST(9), LAMF(9,9,2), PRECIS, LAMNMX(2,2), SIGYY        BLK01400
      COMMON /BLOK/ ENPHI, SIGMAP(19), PCERR(19), VGRID2(252), JGEADD,  BLK01500
     1 NBINOM(20,9), NVGRID(9), NPERG2(9), MDIMA, IBASE, IWTSGN,        BLK01600
     2 MINTER, MLAMMX, NGRID2(9), NLAMMX, IWT, LABEL(80), NINT, NF,     BLK01700
     3 N, IWRITE, IREAD, LINEPG, NMAX, NINTMX, MTRY, ITER, JLAMP1,      BLK01800
     4 NFTRY(9), REGINT, FAILED, PRPREL, NOBASE, WTED, REPEAT,          BLK01900
     5 PLOTRS, PRFINL, NONNEG, FFAIL(9,2), PIVALP(19)                   BLK02000
C                                                                       BLK02100
C***********************************************************************BLK02200
C***********************************************************************BLK02300
C                                                                       BLK02400
C  THE FOLLOWING INSTRUCTIONS DESCRIBE ALL POSSIBLE NECESSARY CHANGES   BLK02500
C  THAT YOU MAY HAVE TO MAKE IN THE BLOCKDATA SUBPROGRAM. (SEE THE      BLK02600
C  MAIN PROGRAM ALSO.)                                                  BLK02700
C                                                                       BLK02800
C***********************************************************************BLK02900
C  ALL STATEMENTS FOR READING CARDS AND FOR PRINTING OUTPUT ARE OF      BLK03000
C  THE FORM READ (IREAD,... AND WRITE(IWRITE,....  THUS, IREAD AND      BLK03100
C  IWRITE MAY HAVE TO BE CHANGED TO THE VALUES APPROPRIATE FOR YOUR     BLK03200
C  COMPUTER CENTER IN THE FOLLOWING STATEMENT -                         BLK03300
C                                                                       BLK03400
      DATA IREAD/5/,IWRITE/6/                                           BLK03500
C                                                                       BLK03600
C***********************************************************************BLK03700
C  IN THE FOLLOWING STATEMENT, LINEPG MUST BE SET TO THE NUMBER OF      BLK03800
C  LINES PER PAGE AVAILABLE FOR PRINTED OUTPUT ON YOUR PRINTER -        BLK03900
C                                                                       BLK04000
      DATA LINEPG/62/                                                   BLK04100
C                                                                       BLK04200
C***********************************************************************BLK04300
C  NMAX MUST BE GREATER THAN OR EQUAL TO THE MAXIMUM NUMBER OF DATA     BLK04400
C  POINTS YOU WILL EVER USE.  CORE STORAGE CAN BE SAVED BY KEEPING      BLK04500
C  NMAX AS LITTLE ABOVE THIS MAXIMUM AS POSSIBLE.  IF YOU CHANGE NMAX   BLK04600
C  IN THE FOLLOWING STATEMENT, YOU MUST ALSO RESET THE DIMENSIONS AS    BLK04700
C  DIRECTED IN THE MAIN PROGRAM -                                       BLK04800
C                                                                       BLK04900
      DATA NMAX/200/                                                    BLK05000
C                                                                       BLK05100
C***********************************************************************BLK05200
C  NINTMX MUST BE GREATER THAN OR EQUAL TO THE MAXIMUM VALUE OF NINT    BLK05300
C  YOU WILL EVER USE WITH REGINT=.TRUE.. IF YOU CHANGE                  BLK05400
C  NINTMX IN THE FOLLOWING STATEMENT, YOU MUST ALSO RESET THE           BLK05500
C  DIMENSIONS AS DIRECTED IN THE MAIN PROGRAM -                         BLK05600
C                                                                       BLK05700
      DATA NINTMX /10/                                                  BLK05800
C                                                                       BLK05900
C  THIS IS THE END OF ALL POSSIBLE NECESSARY CHANGES YOU MAY HAVE TO    BLK06000
C  MAKE IN THE BLOCK DATA SUBPROGRAM                                    BLK06100
C                                                                       BLK06200
C***********************************************************************BLK06300
C***********************************************************************BLK06400
C                                                                       BLK06500
      DATA NGRID2/0, 20, 12, 5*10, 12/,MINTER/7/,                       BLK06600
     1 CONVRG/5.E-5/,MLAMMX/9/, SIGMAP/19*0./, PRECIS/1.E-7/            BLK06700
      END                                                               BLK06800
C--------------------- DOUBLE PRECISION VERSION 2B ---------------------DAT00100
C  SUBROUTINE DATAIO.  FOR INPUT AND OUTPUT OF INITIAL INPUT DATA       DAT00200
C-----------------------------------------------------------------------DAT00300
C  CALLS NO EXTERNAL SUBPROGRAMS.                                       DAT00400
C-----------------------------------------------------------------------DAT00500
      SUBROUTINE DATAIO (T,Y,SQRTW,NMAX,TSTART,DELTAT,NT,NINTMX,LAST)   DAT00600
      LOGICAL REGINT, LAST, NONNEG, NOBASE, PRY, PRPREL, PLOTRS, REPEAT,DAT00700
     1 LBLOKA, LBLOKB, LBLOKC, ABORT, PRFINL                            DAT00800
      DOUBLE PRECISION TEND, DUB, T, TSTART, DELTAT, DBLOKA             DAT00900
      COMMON /BLOK/ DBLOKA(5933), RBLOKA(471), IBLOKA(201), IWTSGN,     DAT01000
     1 IBLOKB, MLAMMX, IBLOKC(9), NLAMMX, IWT, LABEL(80), NINT,         DAT01100
     2 IBLOKD, N, IWRITE, IREAD, IBLOKE(3), MTRY, IBLOKF(11),           DAT01200
     3 REGINT, LBLOKA, PRPREL, NOBASE, LBLOKB, REPEAT, PLOTRS,          DAT01300
     4 PRFINL, NONNEG, LBLOKC(37)                                       DAT01400
      DIMENSION T(NMAX), Y(NMAX), SQRTW(NMAX), TSTART(NINTMX),          DAT01500
     1 NT(NINTMX), DELTAT(NINTMX), LA(9), LATF(2)                       DAT01600
      DATA LATF/1HT, 1HF/                                               DAT01700
 5010 FORMAT (80A1)                                                     DAT01800
      READ (IREAD,5010) LABEL                                           DAT01900
 5020 FORMAT (38H1DISCRETE - Version 2B (December 1990),10X,80A1)       DAT02000
      WRITE (IWRITE,5020) LABEL                                         DAT02100
      ABORT=.FALSE.                                                     DAT02200
 5100 FORMAT (9A1,3I3)                                                  DAT02300
      READ (IREAD,5100) LA,NLAMMX,IWT,MTRY                              DAT02400
      DO 101 J=1,9                                                      DAT02500
      IF (LA(J).EQ.LATF(1) .OR. LA(J).EQ.LATF(2)) GO TO 101             DAT02600
 5101 FORMAT (/43H THE FOLLOWING CARD HAS NO T OR F IN COLUMN,I2)       DAT02700
      WRITE (IWRITE,5101) J                                             DAT02800
      ABORT=.TRUE.                                                      DAT02900
  101 CONTINUE                                                          DAT03000
      IF (NLAMMX.GE.1 .AND. NLAMMX.LE.MLAMMX) GO TO 102                 DAT03100
 5102 FORMAT (/9H NLAMMX =,I4,21H IS NOT BETWEEN 1 AND,I2)              DAT03200
      WRITE (IWRITE,5102) NLAMMX, MLAMMX                                DAT03300
      ABORT=.TRUE.                                                      DAT03400
  102 IF (IWT.EQ.1 .OR. IABS(IWT).EQ.2 .OR. IABS(IWT).EQ.3              DAT03500
     1 .OR. IWT.EQ.4) GO TO 103                                         DAT03600
 5103 FORMAT (/1X,I4,29H IS AN IMPROPER VALUE FOR IWT)                  DAT03700
      WRITE (IWRITE,5103) IWT                                           DAT03800
      ABORT=.TRUE.                                                      DAT03900
  103 IF (MTRY.GE.1 .AND. MTRY.LE.45) GO TO 104                         DAT04000
 5104 FORMAT (/7H MTRY =,I4,24H IS NOT BETWEEN 1 AND 45)                DAT04100
      WRITE (IWRITE,5104) MTRY                                          DAT04200
      ABORT=.TRUE.                                                      DAT04300
  104 IF (.NOT.ABORT) GO TO 105                                         DAT04400
 5105 FORMAT (/21H ERRONEOUS CARD IMAGE,10X,9A1,3I3)                    DAT04500
      WRITE (IWRITE,5105) LA,NLAMMX,IWT,MTRY                            DAT04600
      STOP                                                              DAT04700
  105 LAST=LA(1) .EQ. LATF(1)                                           DAT04800
      REGINT=LA(2) .EQ. LATF(1)                                         DAT04900
      NOBASE=LA(3) .EQ. LATF(1)                                         DAT05000
      NONNEG=LA(4) .EQ. LATF(1)                                         DAT05100
      PRY=LA(5) .EQ. LATF(1)                                            DAT05200
      PRPREL=LA(6) .EQ. LATF(1)                                         DAT05300
      PRFINL=LA(7) .EQ. LATF(1)                                         DAT05400
      PLOTRS=LA(8) .EQ. LATF(1)                                         DAT05500
      REPEAT=LA(9) .EQ. LATF(1)                                         DAT05600
 5106 FORMAT (/9H   LAST =,L2,4X,8HREGINT =,L2,5X,8HNOBASE =,L2,5X,     DAT05700
     1 8HNONNEG =,L2,5X,5HPRY =,L2,8X,8HPRPREL =,L2,5X,8HPRFINL =,L2,5X,DAT05800
     2 8HPLOTRS =,L2,5X,8HREPEAT =,L2/9H NLAMMX =,I2,4X,8HIWT    =,     DAT05900
     3 I2,5X,8HMTRY   =,I3)                                             DAT06000
      WRITE (IWRITE,5106) LAST,REGINT,NOBASE,NONNEG,PRY,PRPREL,PRFINL,  DAT06100
     1 PLOTRS,REPEAT,NLAMMX,IWT,MTRY                                    DAT06200
      IWTSGN=IWT                                                        DAT06300
      IWT=IABS(IWT)                                                     DAT06400
      IF (REGINT) GO TO 150                                             DAT06500
C-----------------------------------------------------------------------DAT06600
C  READS IN T VALUES (REGINT=.FALSE.).                                  DAT06700
C-----------------------------------------------------------------------DAT06800
      READ (IREAD,5110) N                                               DAT06900
 5110 FORMAT (I5)                                                       DAT07000
      READ (IREAD,5120) (T(J),J=1,N)                                    DAT07100
 5120 FORMAT (5D15.6)                                                   DAT07200
      NINT=1                                                            DAT07300
      TSTART(1)=0.D0                                                    DAT07400
      DELTAT(1)=0.D0                                                    DAT07500
      NT(1)=N                                                           DAT07600
      GO TO 200                                                         DAT07700
C-----------------------------------------------------------------------DAT07800
C  CALCULATES T VALUES FROM TSTART, TEND, AND NT (REGINT=.TRUE.).       DAT07900
C-----------------------------------------------------------------------DAT08000
  150 READ (IREAD,5110) NINT                                            DAT08100
      IF (NINT.LE.NINTMX .AND. NINT.GE.1) GO TO 160                     DAT08200
 5160 FORMAT (///7H NINT =,I5,30H IS NOT BETWEEN 1 AND NINTMX =,I2)     DAT08300
      WRITE (IWRITE,5160) NINT,NINTMX                                   DAT08400
      STOP                                                              DAT08500
 5165 FORMAT (///9X,6HTSTART,11X,4HTEND,3X,2HNT)                        DAT08600
  160 WRITE (IWRITE,5165)                                               DAT08700
      N=0                                                               DAT08800
      DO 170 J=1,NINT                                                   DAT08900
      READ (IREAD,5170) TSTART(J),TEND,NT(J)                            DAT09000
 5170 FORMAT (2D15.6,I5)                                                DAT09100
      WRITE (IWRITE,5170) TSTART(J),TEND,NT(J)                          DAT09200
      IF (NT(J) .GE. 2) GO TO 175                                       DAT09300
 5175 FORMAT (/24H ABOVE NT IS LESS THAN 2)                             DAT09400
      WRITE (IWRITE,5175)                                               DAT09500
      STOP                                                              DAT09600
  175 DUB=(TEND-TSTART(J))/DBLE(FLOAT(NT(J)-1))                         DAT09700
      N=N+1                                                             DAT09800
      T(N)=TSTART(J)                                                    DAT09900
      L=NT(J)                                                           DAT10000
      DELTAT(J)=DUB                                                     DAT10100
      DO 180 K=2,L                                                      DAT10200
      N=N+1                                                             DAT10300
  180 T(N)=T(N-1)+DUB                                                   DAT10400
  170 CONTINUE                                                          DAT10500
C-----------------------------------------------------------------------DAT10600
C  READS IN Y VALUES.                                                   DAT10700
C-----------------------------------------------------------------------DAT10800
  200 IF (N.GT.2*NLAMMX+3 .AND. N.LE.NMAX) GO TO 205                    DAT10900
 5205 FORMAT (///4H N =,I5,37H IS NOT BETWEEN 2*NLAMMX+3 AND NMAX =,I5) DAT11000
      WRITE (IWRITE,5205) N,NMAX                                        DAT11100
      STOP                                                              DAT11200
  205 IF (PRY) WRITE (IWRITE,5004)                                      DAT11300
 5004 FORMAT (///1X)                                                    DAT11400
 5207 FORMAT (5E15.6)                                                   DAT11500
      READ (IREAD,5207) (Y(J),J=1,N)                                    DAT11600
      IF (IWT .EQ. 4) GO TO 220                                         DAT11700
      IF (.NOT.PRY) GO TO 500                                           DAT11800
 5210 FORMAT (4(19X,1HT,12X,1HY)/)                                      DAT11900
      WRITE (IWRITE,5210)                                               DAT12000
 5220 FORMAT (4(1PD20.5,E13.5))                                         DAT12100
      WRITE (IWRITE,5220) (T(J),Y(J),J=1,N)                             DAT12200
      GO TO 500                                                         DAT12300
C-----------------------------------------------------------------------DAT12400
C  READS IN SPECIAL WEIGHTS (IWT=4).                                    DAT12500
C-----------------------------------------------------------------------DAT12600
  220 READ (IREAD,5207) (SQRTW(J),J=1,N)                                DAT12700
      DO 227 J=1,N                                                      DAT12800
  227 SQRTW(J)=SQRT(SQRTW(J))                                           DAT12900
      IF (.NOT.PRY) GO TO 500                                           DAT13000
 5230 FORMAT (3(17X,1HT,12X,1HY,5X,8HSQRT(WT))/)                        DAT13100
      WRITE (IWRITE,5230)                                               DAT13200
 5240 FORMAT (3(1PD18.5,2E13.5))                                        DAT13300
      WRITE (IWRITE,5240) (T(J),Y(J),SQRTW(J),J=1,N)                    DAT13400
  500 RETURN                                                            DAT13500
      END                                                               DAT13600
C--------------------- DOUBLE PRECISION VERSION 2B ---------------------WT000100
C  SUBROUTINE WEIGHT.  DOES INITIAL LEAST SQUARES FITS TO RAW DATA TO   WT000200
C      GENERATE WEIGHTS, EVALUATES DATA TRANSFORMS, TRANSFORM AND       WT000300
C      EXPONENTIAL SUMS AT GRID POINTS FOR INTERPOLATION, AND OTHER     WT000400
C      QUANTITIES FOR LATER USE.                                        WT000500
C  IF FINAL=.TRUE. USES LAMF (VALUES FROM FANLYZ).                      WT000600
C  IF FINAL=.FALSE., GENERATES ITS OWN CRUDE STARTING LAMBDAS.          WT000700
C-----------------------------------------------------------------------WT000800
C  CALLS SUBPROGRAMS - LSTSQR, ETHEOR                                   WT000900
C  WHICH IN TURN CALL - EVAR, PIVOT, PIVOT1, ANLERR                     WT001000
C-----------------------------------------------------------------------WT001100
      SUBROUTINE WEIGHT (FINAL,Y,T,SQRTW,TIMU,E,NMAX,NDIME,             WT001200
     1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG)                               WT001300
      EXTERNAL EVAR                                                     WT001400
      LOGICAL FINAL, REGINT, NOBASE, WTED, FAILED, LDUM, LBLOKA, PRPREL,WT001500
     1 LBLOKB, FFAIL                                                    WT001600
      REAL LAMNMX, LAMST, LAMF                                          WT001700
      DOUBLE PRECISION VARBES, RGRIDE, GEST, GFST, MU, DUB, DDUB,       WT001800
     1 DDDUB, RGRIDF, DS, DSS, DSSS, TRBASE, GF, DGFL, DGRIDF,          WT001900
     2 GFSTL, FHAT, VARMIN, SUMWTY, GSE, SUMWT, ZERO, DGRIDE, GESTL,    WT002000
     3 VAR, PALPHA, PLAM, TIMU, E, TSTART, DELTAT, DBLOKA, T            WT002100
      COMMON /BLOK/ TRBASE(19), GF(19,126), DGFL(19,126), DGRIDF,       WT002200
     1 GFSTL, FHAT(19), VARMIN, SUMWTY, SUMWT, ZERO(10), DGRIDE,        WT002300
     2 GESTL, VAR, PALPHA(10), PLAM(10), DBLOKA(1069), DELTA, FHATMX,   WT002400
     3 CONVRG, LAMST(9), LAMF(9,9,2), RBLOKA, LAMNMX(2,2),              WT002500
     4 RBLOKB(292), JGEADD, NBINOM(20,9), NVGRID(9), NPERG2(9),         WT002600
     5 MDIMA, IBASE, IWTSGN, MINTER, MLAMMX, NGRID2(9), NLAMMX          WT002700
      COMMON /BLOK/ IWT, LABEL(80), NINT, NF, N, IWRITE, IBLOKA(16),    WT002800
     1 REGINT, FAILED, PRPREL, NOBASE, WTED, LBLOKA(4), FFAIL(9,2),     WT002900
     2 LBLOKB(19)                                                       WT003000
      DIMENSION E(NDIME,10), WTALP(10), WTLAM(10), MU(19)               WT003100
      DIMENSION Y(NMAX), T(NMAX), SQRTW(NMAX), TIMU(NMAX),              WT003200
     1 TSTART(NINTMX), DELTAT(NINTMX), NT(NINTMX), GSE(NDIMG,4),        WT003300
     2 ALPH(10,2), TNLAMN(2), GAMMA(21)                                 WT003400
      DATA ALPH /1HT, 1HR, 1HA, 1HN, 1HS, 1HF, 1HO, 1HR, 1HM, 1HS,      WT003500
     1 1HR, 1HA, 1HW, 1H , 1HD, 1HA, 1HT, 1HA, 1H , 1H /,               WT003600
     2 TNLAMN/.2,.02/,T1LAMX/2.08/,NGRID1/240/,NG2MAX/20/               WT003700
      IF (FINAL) GO TO 120                                              WT003800
      FHATMX=0.                                                         WT003900
      MDIMA=2*MLAMMX+1                                                  WT004000
      K=MLAMMX+1                                                        WT004100
      DO 105 J=1,K                                                      WT004200
  105 ZERO(J)=0.D0                                                      WT004300
      TN=-1.E+20                                                        WT004400
      DDUM=-1.E+20                                                      WT004500
      DO 106 K=1,5                                                      WT004600
      DUM=1.E+20                                                        WT004700
      DO 107 J=1,N                                                      WT004800
      DDDUM=SNGL(T(J))                                                  WT004900
      IF (K .EQ. 1) TN=AMAX1(TN,DDDUM)                                  WT005000
  107 IF (DDDUM.GT.DDUM .AND. DDDUM.LT.DUM) DUM=DDDUM                   WT005100
      DDUM=DUM+1.E-5*ABS(DUM)                                           WT005200
      IF (K .EQ. 1) T1=DUM                                              WT005300
      IF (K .EQ. 2) T2=DUM                                              WT005400
  106 CONTINUE                                                          WT005500
      T5=DUM                                                            WT005600
      TNL1=-1.E+20                                                      WT005700
      DUM=TN-1.E-5*ABS(TN)                                              WT005800
      DO 108 J=1,N                                                      WT005900
  108 IF (SNGL(T(J)) .LT. DUM) TNL1=AMAX1(TNL1,SNGL(T(J)))              WT006000
      DUM=.25*(T5-T1)                                                   WT006100
      DELTA=DUM-T1                                                      WT006200
C-----------------------------------------------------------------------WT006300
C  LAMNMX(J,K) = MIN. (J=1) AND MAX. (J=2) ALLOWED VALUES OF LAMBDA     WT006400
C      FOR LEAST SQUARES PARAMETERS (K=1) AND FOR STARTING VALUES (K=2).WT006500
C-----------------------------------------------------------------------WT006600
      LAMNMX(2,1)=T1LAMX/(T1+DELTA)                                     WT006700
      LAMNMX(2,2)=AMIN1(LAMNMX(2,1),.693/(T1+DELTA))                    WT006800
      LAMNMX(1,2)=TNLAMN(1)/(TN+DELTA)                                  WT006900
      IF (NOBASE) GO TO 110                                             WT007000
      NPARMX=2*NLAMMX+1                                                 WT007100
      IBASE=1                                                           WT007200
      LAMNMX(1,1)=LAMNMX(1,2)                                           WT007300
      GO TO 115                                                         WT007400
  110 IBASE=0                                                           WT007500
      NPARMX=2*NLAMMX                                                   WT007600
      FHAT(NPARMX+1)=0.D0                                               WT007700
      LAMNMX(1,1)=TNLAMN(2)/(TN+DELTA)                                  WT007800
C-----------------------------------------------------------------------WT007900
C  IF NECESSARY, INCREASE DELTA SO THAT DISTANCE BETWEEN FIRST TWO      WT008000
C      POINTS ON X=ALOG(T+DELTA) AXIS IS LESS THAN A HALF-PERIOD OF     WT008100
C      THE HIGHEST FREQUENCY USED IN THE TRANSFORMS                     WT008200
C-----------------------------------------------------------------------WT008300
  115 XRANGE=ALOG((TN+DELTA)**2/((T1+DELTA)*(TNL1+DELTA)))              WT008400
      IF (ALOG((T2+DELTA)/(T1+DELTA))*FLOAT(2*NLAMMX) .LT. XRANGE)      WT008500
     1 GO TO 150                                                        WT008600
      DELTA=DELTA+DUM                                                   WT008700
      GO TO 115                                                         WT008800
  120 IWT=IABS(IWTSGN)                                                  WT008900
      IF (IWT.EQ.1 .OR. IWT.EQ.4) RETURN                                WT009000
      DO 130 J=1,N                                                      WT009100
  130 Y(J)=Y(J)/SQRTW(J)                                                WT009200
  150 IF (IWT .EQ. 4) GO TO 400                                         WT009300
      SUMWT=DBLE(FLOAT(N))                                              WT009400
      SUMWTY=0.D0                                                       WT009500
      DUM=0.                                                            WT009600
      DO 160 J=1,N                                                      WT009700
      SUMWTY=SUMWTY+DBLE(Y(J))                                          WT009800
      DUM=DUM+Y(J)*Y(J)                                                 WT009900
  160 SQRTW(J)=1.                                                       WT010000
      VARMIN=DBLE(DUM*CONVRG*CONVRG)                                    WT010100
      IF (IWT.EQ.1) GO TO 400                                           WT010200
      WTED=.FALSE.                                                      WT010300
      VARBES=1.D+20                                                     WT010400
      RATIOL=LAMNMX(2,2)/LAMNMX(1,2)                                    WT010500
C-----------------------------------------------------------------------WT010600
C  START OF MAIN LOOP FOR STEPWISE ANALYSIS                             WT010700
C-----------------------------------------------------------------------WT010800
      DO 200 JLAM=1,NLAMMX                                              WT010900
      NF=2*JLAM+IBASE                                                   WT011000
 5200 FORMAT (1H1,25X,80A1//30X,                                        WT011100
     1 41HPRELIMINARY ANALYSIS TO DETERMINE WEIGHTS/30X,                WT011200
     2 17HANALYSIS ASSUMING,I2,11H COMPONENTS)                          WT011300
      IF (PRPREL) WRITE (IWRITE,5200) LABEL,JLAM                        WT011400
      IF (FINAL) GO TO 220                                              WT011500
      FFAIL(JLAM,1)=.TRUE.                                              WT011600
      DUM=RATIOL**(1./FLOAT(JLAM+1))                                    WT011700
      LAMST(1)=LAMNMX(1,2)*DUM                                          WT011800
      IF (JLAM .EQ. 1) GO TO 240                                        WT011900
      DO 210 J=2,JLAM                                                   WT012000
  210 LAMST(J)=LAMST(J-1)*DUM                                           WT012100
      GO TO 240                                                         WT012200
  220 DO 230 J=1,JLAM                                                   WT012300
  230 LAMST(J)=LAMF(J,JLAM,1)                                           WT012400
  240 CALL LSTSQR (LAMST,JLAM,.TRUE.,N,Y,T,JLAM+IBASE,EVAR,GSE,NDIMG,   WT012500
     1 SQRTW,TIMU,E,NMAX,TSTART,DELTAT,NT,NINT,1,FFAIL(JLAM,1),NDIME,   WT012600
     2 .FALSE.,PRPREL)                                                  WT012700
      IF (FAILED .OR. VAR.GE.VARBES) GO TO 300                          WT012800
      VARBES=VAR                                                        WT012900
      NLINWT=JLAM+IBASE                                                 WT013000
      DO 260 J=1,NLINWT                                                 WT013100
      WTALP(J)=SNGL(PALPHA(J))                                          WT013200
  260 WTLAM(J)=SNGL(PLAM(J))                                            WT013300
      IF (IWTSGN .LT. 0) NLINWT=JLAM                                    WT013400
  200 CONTINUE                                                          WT013500
  300 IF (VARBES .LT. 1.D+20) GO TO 310                                 WT013600
 5300 FORMAT (////34H 1-COMPONENT ANALYSIS TO DETERMINE,                WT013700
     1 54H WEIGHTS SOMEHOW FAILED - UNIT WEIGHTS WILL BE USED.  ,       WT013800
     2 22(2H**))                                                        WT013900
      WRITE (IWRITE,5300)                                               WT014000
      IWT=1                                                             WT014100
      GO TO 400                                                         WT014200
C-----------------------------------------------------------------------WT014300
C  CALCULATE WEIGHTS INCLUDING ERRFIT (STANDARD ERROR IN REGION OF      WT014400
C      THE THEORETICAL CURVE WHERE IT IS CLOSEST TO Y=0., WHERE A       WT014500
C      WEIGHT COULD OTHERWISE BE DISASTROUSLY LARGE IF ERRFIT WERE      WT014600
C      NOT ACCOUNTED FOR).                                              WT014700
C-----------------------------------------------------------------------WT014800
  310 DUM=1.E+20                                                        WT014900
      DO 320 K=1,N                                                      WT015000
      S=0.                                                              WT015100
      DDUM=SNGL(T(K))                                                   WT015200
      DO 330 J=1,NLINWT                                                 WT015300
  330 S=S+WTALP(J)*EXP(-DDUM*WTLAM(J))                                  WT015400
      SQRTW(K)=S                                                        WT015500
      IF (ABS(S) .GE. DUM) GO TO 320                                    WT015600
      DUM=ABS(S)                                                        WT015700
      L=K                                                               WT015800
  320 CONTINUE                                                          WT015900
      K=MIN0(N,L+5)                                                     WT016000
      I=MAX0(1,K-9)                                                     WT016100
      DUM=0.                                                            WT016200
      DDUM=0.                                                           WT016300
      IF (IWTSGN.LT.0 .AND. .NOT.NOBASE) DDUM=WTALP(NLINWT+1)           WT016400
      DO 340 J=I,K                                                      WT016500
  340 DUM=DUM+(Y(J)-SQRTW(J)-DDUM)**2                                   WT016600
      ERRFIT=SQRT(DUM/FLOAT(K-I+1))                                     WT016700
      DO 350 K=1,N                                                      WT016800
      SQRTW(K)=1./(ABS(SQRTW(K))+ERRFIT)                                WT016900
      IF (IWT .EQ. 2) SQRTW(K)=SQRT(SQRTW(K))                           WT017000
  350 CONTINUE                                                          WT017100
      L=1                                                               WT017200
      IF (FINAL) L=2                                                    WT017300
 5350 FORMAT (/////41H PARAMETERS USED TO GENERATE WEIGHTS FOR ,        WT017400
     1 10A1//5X,6HLAMBDA,9X,5HALPHA,10X,                                WT017500
     2 49HERRFIT (UNCERTAINTY TERM ADDED TO ABSOLUTE VALUES,            WT017600
     3 24H OF THEORETICAL CURVE) =,1PE9.2/(1X,1PE10.3,E14.3))           WT017700
      WRITE (IWRITE,5350) (ALPH(J,L),J=1,10),ERRFIT,(WTLAM(K),WTALP(K), WT017800
     1 K=1,NLINWT)                                                      WT017900
  400 WTED=IWT.NE.1                                                     WT018000
      IF (FINAL) GO TO 600                                              WT018100
C-----------------------------------------------------------------------WT018200
C  INITIAL COMPUTATION OF QUANTITIES TO BE USED LATER                   WT018300
C-----------------------------------------------------------------------WT018400
      DGRIDE=DBLE(ALOG(LAMNMX(2,2)/LAMNMX(1,2))/FLOAT(NGRID1-1))        WT018500
      RGRIDE=DEXP(DGRIDE)                                               WT018600
      JGEADD=MINTER/2                                                   WT018700
      IF (NOBASE) JGEADD=JGEADD+INT(ALOG(LAMNMX(1,2)/LAMNMX(1,1))/      WT018800
     1 SNGL(DGRIDE)+.5)                                                 WT018900
      GEST=DBLE(LAMNMX(1,2)/SNGL(RGRIDE)**(JGEADD+1))                   WT019000
      GESTL=DLOG(GEST)                                                  WT019100
      NGE1=NGRID1+JGEADD+INT(ALOG(2.*LAMNMX(2,1)/LAMNMX(2,2))/          WT019200
     1 SNGL(DGRIDE)+.5)+MINTER/2                                        WT019300
      IF (NGE1 .LE. 500) GO TO 401                                      WT019400
 5401 FORMAT (////20H FIX LAMNMX.  NGE1 =,I9,1X,51(2H**))               WT019500
      WRITE (IWRITE,5401) NGE1                                          WT019600
      STOP                                                              WT019700
  401 JGEADD=JGEADD+1                                                   WT019800
      IF (NLAMMX .EQ. 1) GO TO 403                                      WT019900
      DO 402 J=2,NLAMMX                                                 WT020000
  402 NPERG2(J)=NGRID1/NGRID2(J)                                        WT020100
  403 GAMMA(1)=1.                                                       WT020200
      DO 404 J=1,NG2MAX                                                 WT020300
  404 GAMMA(J+1)=GAMMA(J)*FLOAT(J)                                      WT020400
      DO 405 J=1,NLAMMX                                                 WT020500
      DO 406 K=J,NG2MAX                                                 WT020600
      L=K-J+1                                                           WT020700
  406 NBINOM(K,J)=INT(GAMMA(K)/(GAMMA(L)*GAMMA(J))+.5)                  WT020800
      IF (J .EQ. 1) GO TO 405                                           WT020900
      L=NGRID2(J)+1                                                     WT021000
      K=L-J                                                             WT021100
      NVGRID(J)=INT(GAMMA(L)/(GAMMA(J+1)*GAMMA(K))+.5)                  WT021200
  405 CONTINUE                                                          WT021300
      NGF=NGRID1/2+2*(MINTER/2)                                         WT021400
      DGRIDF=DBLE(ALOG(LAMNMX(2,1)/LAMNMX(1,1))/FLOAT(NGRID1/2-1))      WT021500
      RGRIDF=DEXP(DGRIDF)                                               WT021600
      GFST=DBLE(LAMNMX(1,1)/SNGL(RGRIDF)**(MINTER/2+1))                 WT021700
      GFSTL=DLOG(GFST)                                                  WT021800
C-----------------------------------------------------------------------WT021900
C  COMPUTE FHAT (TRANSFORMS OF DATA) AND GF AND DGFL (INTERPOLATION     WT022000
C      GRID POINTS FOR TRANSFORMS AND THEIR DERIVATIVES W.R.T.          WT022100
C      ALOG(LAMBDA))                                                    WT022200
C-----------------------------------------------------------------------WT022300
      IF (.NOT.REGINT) GO TO 420                                        WT022400
      DO 410 J=1,NINT                                                   WT022500
  410 TSTART(J)=TSTART(J)+DBLE(DELTA)                                   WT022600
  420 DUB=DBLE(DELTA)                                                   WT022700
      DO 430 J=1,N                                                      WT022800
  430 T(J)=T(J)+DUB                                                     WT022900
      MU(1)=DBLE(6.283185/XRANGE)                                       WT023000
      DO 440 K=1,NPARMX                                                 WT023100
      IF (K .GT. 1) GO TO 450                                           WT023200
      DO 445 J=1,N                                                      WT023300
      FHATMX=FHATMX+SQRTW(J)*ABS(Y(J))                                  WT023400
  445 TIMU(J)=DBLE(SQRTW(J))                                            WT023500
      GO TO 470                                                         WT023600
  450 L=K/2                                                             WT023700
      IF (L+L .LT. K) GO TO 460                                         WT023800
      DUB=MU(1)                                                         WT023900
      IF (L .GT. 1) DUB=MU(L-1)+DUB                                     WT024000
      MU(L)=DUB                                                         WT024100
      DO 455 J=1,N                                                      WT024200
  455 TIMU(J)=DCOS(DLOG(T(J))*DUB)*DBLE(SQRTW(J))                       WT024300
      GO TO 470                                                         WT024400
  460 DUB=MU(L)                                                         WT024500
      DO 465 J=1,N                                                      WT024600
  465 TIMU(J)=DSIN(DLOG(T(J))*DUB)*DBLE(SQRTW(J))                       WT024700
  470 DUB=0.D0                                                          WT024800
      DDUB=0.D0                                                         WT024900
      DO 475 J=1,N                                                      WT025000
      DUB=DUB+Y(J)*TIMU(J)                                              WT025100
  475 DDUB=DDUB+TIMU(J)                                                 WT025200
      FHAT(K)=DUB                                                       WT025300
      TRBASE(K)=DDUB                                                    WT025400
      DUB=GFST                                                          WT025500
      DO 480 L=1,NGF                                                    WT025600
      DUB=DUB*RGRIDF                                                    WT025700
      DS=0.D0                                                           WT025800
      DSS=0.D0                                                          WT025900
      IF (REGINT) GO TO 500                                             WT026000
      DO 490 J=1,N                                                      WT026100
      DDUB=DEXP(-DUB*T(J))*TIMU(J)                                      WT026200
      DS=DS+DDUB                                                        WT026300
  490 DSS=DSS-DDUB*T(J)                                                 WT026400
      GO TO 550                                                         WT026500
C-----------------------------------------------------------------------WT026600
C  USE RECURSION RELATIONS IF REGINT=.TRUE.                             WT026700
C-----------------------------------------------------------------------WT026800
  500 J=0                                                               WT026900
      DO 510 JJ=1,NINT                                                  WT027000
      DDUB=DEXP(-DUB*DELTAT(JJ))                                        WT027100
      DDDUB=DEXP(-DUB*TSTART(JJ))                                       WT027200
      II=NT(JJ)                                                         WT027300
      DO 520 I=1,II                                                     WT027400
      J=J+1                                                             WT027500
      DSSS=DDDUB*TIMU(J)                                                WT027600
      DS=DS+DSSS                                                        WT027700
      DSS=DSS-T(J)*DSSS                                                 WT027800
      DDDUB=DDDUB*DDUB                                                  WT027900
  520 CONTINUE                                                          WT028000
  510 CONTINUE                                                          WT028100
  550 GF(K,L)=DS                                                        WT028200
      DGFL(K,L)=DSS*DUB                                                 WT028300
  480 CONTINUE                                                          WT028400
  440 CONTINUE                                                          WT028500
      IF (PRPREL) WRITE (IWRITE,5560) DELTA,FHAT(1),(MU(K),FHAT(2*K),   WT028600
     1 FHAT(2*K+1),K=1,NLAMMX)                                          WT028700
 5560 FORMAT (////6X,2HMU,13X,                                          WT028800
     1 38HREAL AND IMAGINARY PARTS OF TRANSFORMS,15X,7HDELTA =,1PE9.2   WT028900
     2 /4X,4H.000,1PD17.3/(1X,0PF7.3,1PD17.3,D14.3))                    WT029000
      DUB=DBLE(DELTA)                                                   WT029100
      DO 560 J=1,N                                                      WT029200
  560 T(J)=T(J)-DUB                                                     WT029300
      IF (.NOT.REGINT) GO TO 600                                        WT029400
      DO 570 J=1,NINT                                                   WT029500
  570 TSTART(J)=TSTART(J)-DBLE(DELTA)                                   WT029600
  600 IF (.NOT.WTED) GO TO 620                                          WT029700
      DUM=0.                                                            WT029800
      SUMWT=0.D0                                                        WT029900
      SUMWTY=0.D0                                                       WT030000
      DO 610 J=1,N                                                      WT030100
      Y(J)=Y(J)*SQRTW(J)                                                WT030200
      DUM=DUM+Y(J)*Y(J)                                                 WT030300
      IF (NOBASE) GO TO 610                                             WT030400
      DUB=DBLE(SQRTW(J))                                                WT030500
      SUMWT=SUMWT+DUB*DUB                                               WT030600
      SUMWTY=SUMWTY+DUB*DBLE(Y(J))                                      WT030700
  610 CONTINUE                                                          WT030800
      VARMIN=DBLE(DUM*CONVRG*CONVRG)                                    WT030900
  620 IF (FINAL .OR. .NOT.REGINT) GO TO 900                             WT031000
C-----------------------------------------------------------------------WT031100
C  COMPUTE GSE (INTERPOLATION GRID POINTS FOR EXPONENTIAL SUMS)         WT031200
C-----------------------------------------------------------------------WT031300
      DUB=GEST                                                          WT031400
      L=NGE1-INT(.69315/SNGL(DGRIDE))                                   WT031500
      DO 630 K=1,NGE1                                                   WT031600
      DUB=DUB*RGRIDE                                                    WT031700
      LDUM=K .LE. L                                                     WT031800
      CALL ETHEOR (DUB,.TRUE.,WTED,T,SQRTW,Y,NMAX,TSTART,DELTAT,NT,     WT031900
     1 NINT,DS,DSS,DSSS,DDDUB,.TRUE.,LDUM)                              WT032000
      GSE(K,1)=DS                                                       WT032100
      GSE(K,2)=DSS                                                      WT032200
      GSE(K,3)=DSSS                                                     WT032300
      GSE(K,4)=DDDUB                                                    WT032400
  630 CONTINUE                                                          WT032500
  900 RETURN                                                            WT032600
      END                                                               WT032700
C--------------------- DOUBLE PRECISION VERSION 2B ---------------------FAN00100
C  SUBROUTINE FANLYZ.  DOES CONSTRAINED STEPWISE LEAST SQUARES ANALYSIS FAN00200
C      OF RAW DATA (USING GRID SEARCH IF NECESSARY) TO GET STARTING     FAN00300
C      VALUES FOR ANALYSIS OF TRANSFORMS, WHICH YEILD LAMF (STARTING    FAN00400
C      VALUES FOR FINAL ANALYSIS OF RAW DATA IN YANLYZ).                FAN00500
C-----------------------------------------------------------------------FAN00600
C  CALLS SUBPROGRAMS - LSTSQR, EVAR                                     FAN00700
C  WHICH IN TURN CALL - PIVOT, PIVOT1, ETHEOR, VARF                     FAN00800
C-----------------------------------------------------------------------FAN00900
      SUBROUTINE FANLYZ (T,Y,SQRTW,YLYFIT,E,NMAX,NDIME,                 FAN01000
     1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG)                               FAN01100
      EXTERNAL VARF, EVAR                                               FAN01200
      LOGICAL PRPREL, WTED, FAILED, SEARCH, LDUM, FFAIL, LBLOKA,        FAN01300
     1 LBLOKB, LBLOKC, LBLOKD, SAVEPK                                   FAN01400
      REAL LAMNMX, LAMST, LAMF                                          FAN01500
      DOUBLE PRECISION YLYFIT, E, TSTART, DELTAT, GSE, DUB, DDUB, T,    FAN01600
     1 DGRIDE, GESTL, VAR, PALPHA, PLAM, VARE, VARTRY, DBLOKA, DBLOKB   FAN01700
      COMMON /BLOK/ DBLOKA(4841), DGRIDE, GESTL, VAR, PALPHA(10),       FAN01800
     1 PLAM(10), DBLOKB(1069), RBLOKA(3), LAMST(9), LAMF(9,9,2),        FAN01900
     2 RBLOKB, LAMNMX(2,2), RBLOKC(40), VGRID2(252), JGEADD,            FAN02000
     3 IBLOKA(180), NVGRID(9), NPERG2(9), IBLOKB, IBASE, IBLOKC(2),     FAN02100
     4 MLAMMX, NGRID2(9), NLAMMX, IWT, LABEL(80), NINT, NF, N,          FAN02200
     5 IWRITE, IBLOKD(4), MTRY, IBLOKE(2), NFTRY(9), LBLOKA, FAILED     FAN02300
      COMMON /BLOK/ PRPREL, LBLOKB, WTED, LBLOKC(4), FFAIL(9,2),        FAN02400
     1 LBLOKD(19)                                                       FAN02500
      DIMENSION E(NDIME,10), LENDEQ(9), LEQ(9), JCENTR(20), LSTART(9,45)FAN02600
      DIMENSION T(NMAX), Y(NMAX), SQRTW(NMAX), GSE(NDIMG,4),            FAN02700
     1 YLYFIT(NMAX), TSTART(NINTMX), DELTAT(NINTMX), NT(NINTMX)         FAN02800
      EQUIVALENCE (LENDEQ(1),LEND1), (LENDEQ(2),LEND2), (LENDEQ(3),     FAN02900
     1 LEND3), (LENDEQ(4),LEND4), (LENDEQ(5),LEND5), (LENDEQ(6),LEND6), FAN03000
     2 (LENDEQ(7),LEND7), (LENDEQ(8),LEND8), (LENDEQ(9),LEND9),         FAN03100
     3 (LEQ(1),L1), (LEQ(2),L2), (LEQ(3),L3), (LEQ(4),L4), (LEQ(5),L5), FAN03200
     4 (LEQ(6),L6), (LEQ(7),L7), (LEQ(8),L8), (LEQ(9),L9)               FAN03300
      DATA LLPEAK/1/                                                    FAN03400
      WTED=IWT .NE. 1                                                   FAN03500
      LAMST(1)=SQRT(LAMNMX(2,2)*LAMNMX(1,2))                            FAN03600
C-----------------------------------------------------------------------FAN03700
C  START OF MAIN LOOP FOR STEPWISE ANALYSIS                             FAN03800
C-----------------------------------------------------------------------FAN03900
      DO 200 JLAM=1,NLAMMX                                              FAN04000
C-----------------------------------------------------------------------FAN04100
C  INITIALIZE FOR FITS ASSUMING JLAM COMPONENTS                         FAN04200
C-----------------------------------------------------------------------FAN04300
 5150 FORMAT (1H1,25X,80A1/)                                            FAN04400
      IF (PRPREL) WRITE (IWRITE,5150) LABEL                             FAN04500
      VARE=1.D+20                                                       FAN04600
      VARTRY=1.D+20                                                     FAN04700
      NF=2*JLAM+IBASE                                                   FAN04800
      IF (JLAM .EQ. 1) GO TO 220                                        FAN04900
      LVGRID=NVGRID(JLAM)                                               FAN05000
      DO 210 J=1,LVGRID                                                 FAN05100
  210 VGRID2(J)=0.                                                      FAN05200
  220 JLAML1=JLAM-1                                                     FAN05300
      JLAMP1=JLAM+1                                                     FAN05400
      IF (PRPREL) WRITE (IWRITE,5210) JLAM                              FAN05500
 5210 FORMAT (30X,45HINITIAL ANALYSIS OF THE RAW DATA AND THEN THE,     FAN05600
     2 22H TRANSFORMS - ASSUMING,I2,11H COMPONENTS)                     FAN05700
C-----------------------------------------------------------------------FAN05800
C  START OF LOOP FOR UP TO MTRY JLAM-COMPONENT FITS                     FAN05900
C-----------------------------------------------------------------------FAN06000
      DO 300 JTRY=1,MTRY                                                FAN06100
 5310 FORMAT (////30X,22HFIT USING STARTING SET,I3)                     FAN06200
      IF (PRPREL) WRITE (IWRITE,5310) JTRY                              FAN06300
      IF (JLAM.EQ.1) GO TO 600                                          FAN06400
      IF (JTRY-2) 315,350,500                                           FAN06500
C-----------------------------------------------------------------------FAN06600
C  SELECT FIRST SET OF STARTING VALUES                                  FAN06700
C-----------------------------------------------------------------------FAN06800
  315 DO 320 J=1,JLAML1                                                 FAN06900
  320 LAMST(J)=LAMF(J,JLAML1,2)                                         FAN07000
      IF (LPEAK.EQ.1 .OR. LPEAK.EQ.JLAM) GO TO 345                      FAN07100
      LAMST(JLAM)=SQRT(LAMF(LPEAK-1,JLAML1,2)*LAMF(LPEAK,JLAML1,2))     FAN07200
      GO TO 600                                                         FAN07300
  345 IF (LPEAK .EQ. 1) LAMST(JLAM)=SQRT(LAMNMX(1,2)*LAMF(1,JLAML1,2))  FAN07400
      IF (LPEAK .EQ. JLAM) LAMST(JLAM)=SQRT(LAMF(JLAML1,JLAML1,2)*      FAN07500
     1 LAMNMX(2,2))                                                     FAN07600
      GO TO 600                                                         FAN07700
C-----------------------------------------------------------------------FAN07800
C  DO COARSE GRID SEARCH                                                FAN07900
C-----------------------------------------------------------------------FAN08000
  350 L=NGRID2(JLAM)-JLAM                                               FAN08100
      DO 360 J=1,JLAM                                                   FAN08200
      L=L+1                                                             FAN08300
  360 LENDEQ(J)=L                                                       FAN08400
      IF (JLAM .EQ. MLAMMX) GO TO 368                                   FAN08500
      DO 365 J=JLAMP1,MLAMMX                                            FAN08600
  365 LENDEQ(J)=LENDEQ(JLAM)                                            FAN08700
  368 JCENTR(1)=JGEADD-1+(NPERG2(JLAM)+1)/2                             FAN08800
      L=NGRID2(JLAM)                                                    FAN08900
      DO 370 J=2,L                                                      FAN09000
  370 JCENTR(J)=JCENTR(J-1)+NPERG2(JLAM)                                FAN09100
      SEARCH=.TRUE.                                                     FAN09200
C-----------------------------------------------------------------------FAN09300
C  ENTRY POINT TO FIND STARTING LAMBDA VALUES FROM GRID                 FAN09400
C-----------------------------------------------------------------------FAN09500
  400 LOC=0                                                             FAN09600
      DO 410 L1=1,LEND1                                                 FAN09700
      LST2=L1+1                                                         FAN09800
      DO 420 L2=LST2,LEND2                                              FAN09900
      LST3=L2+1                                                         FAN10000
      IF (JLAM .LT. 3) LST3=LEND3                                       FAN10100
      DO 430 L3=LST3,LEND3                                              FAN10200
      LST4=L3+1                                                         FAN10300
      IF (JLAM .LT. 4) LST4=LEND4                                       FAN10400
      DO 440 L4=LST4,LEND4                                              FAN10500
      LST5=L4+1                                                         FAN10600
      IF (JLAM .LT. 5) LST5=LEND5                                       FAN10700
      DO 450 L5=LST5,LEND5                                              FAN10800
      LST6=L5+1                                                         FAN10900
      IF (JLAM .LT. 6) LST6=LEND6                                       FAN11000
      DO 460 L6=LST6,LEND6                                              FAN11100
      LST7=L6+1                                                         FAN11200
      IF (JLAM .LT. 7) LST7=LEND7                                       FAN11300
      DO 470 L7=LST7,LEND7                                              FAN11400
      LST8=L7+1                                                         FAN11500
      IF (JLAM .LT. 8) LST8=LEND8                                       FAN11600
      DO 480 L8=LST8,LEND8                                              FAN11700
      LST9=L8+1                                                         FAN11800
      IF (JLAM .LT. 9) LST9=LEND9                                       FAN11900
      DO 490 L9=LST9,LEND9                                              FAN12000
      LOC=LOC+1                                                         FAN12100
      IF (SEARCH) GO TO 498                                             FAN12200
      IF (VGRID2(LOC) .GT. VGMAX) GO TO 496                             FAN12300
      MINDIF=1000                                                       FAN12400
      IF (JTRY .LE. 2) GO TO 493                                        FAN12500
      L=JTRY-1                                                          FAN12600
      DO 491 J=2,L                                                      FAN12700
      I=0                                                               FAN12800
      DO 492 K=1,JLAM                                                   FAN12900
  492 I=I+IABS(LEQ(K)-LSTART(K,J))                                      FAN13000
      MINDIF=MIN0(MINDIF,I)                                             FAN13100
  491 CONTINUE                                                          FAN13200
      IF (MINDIF-MAXDIF) 496,493,494                                    FAN13300
  493 IF (VGRID2(LOC) .GE. VGBEST) GO TO 496                            FAN13400
  494 VGBEST=VGRID2(LOC)                                                FAN13500
      MAXDIF=MINDIF                                                     FAN13600
      DO 495 J=1,JLAM                                                   FAN13700
  495 LSTART(J,JTRY)=LEQ(J)                                             FAN13800
  496 IF (LOC .LT. LVGRID) GO TO 490                                    FAN13900
      DO 497 J=1,JLAM                                                   FAN14000
      L=LSTART(J,JTRY)                                                  FAN14100
  497 LAMST(J)=EXP(SNGL(GESTL+FLOAT(JCENTR(L))*DGRIDE))                 FAN14200
      GO TO 600                                                         FAN14300
  498 DO 499 J=1,JLAM                                                   FAN14400
      L=LEQ(J)                                                          FAN14500
  499 PLAM(J)=DEXP(GESTL+FLOAT(JCENTR(L))*DGRIDE)                       FAN14600
      CALL EVAR (0.,JLAM,DUB,T,SQRTW,NMAX,TSTART,DELTAT,NT,             FAN14700
     1 NINT,.FALSE.,LDUM,E,Y,YLYFIT,JLAM+IBASE,1.D+20,0,NDIME,.TRUE.,   FAN14800
     2 GSE,NDIMG,J)                                                     FAN14900
  490 CONTINUE                                                          FAN15000
  480 CONTINUE                                                          FAN15100
  470 CONTINUE                                                          FAN15200
  460 CONTINUE                                                          FAN15300
  450 CONTINUE                                                          FAN15400
  440 CONTINUE                                                          FAN15500
  430 CONTINUE                                                          FAN15600
  420 CONTINUE                                                          FAN15700
  410 CONTINUE                                                          FAN15800
C-----------------------------------------------------------------------FAN15900
C  FIND LOWEST VARIANCE ON COARSE GRID                                  FAN16000
C-----------------------------------------------------------------------FAN16100
  500 VGMAX=1.E+20                                                      FAN16200
      DO 510 J=1,LVGRID                                                 FAN16300
  510 VGMAX=AMIN1(VGMAX,VGRID2(J))                                      FAN16400
      IF (JTRY .GT. 2) VGMAX=10.*VGMAX                                  FAN16500
      MAXDIF=-1                                                         FAN16600
      SEARCH=.FALSE.                                                    FAN16700
      VGBEST=1.E+20                                                     FAN16800
      GO TO 400                                                         FAN16900
C-----------------------------------------------------------------------FAN17000
C  DO LEAST SQUARES FITS                                                FAN17100
C-----------------------------------------------------------------------FAN17200
  600 LDUM=JTRY.EQ.1                                                    FAN17300
      CALL LSTSQR (LAMST,JLAM,.TRUE.,N,Y,T,JLAM+IBASE,                  FAN17400
     1 EVAR,GSE,NDIMG,SQRTW,YLYFIT,E,NMAX,TSTART,DELTAT,NT,NINT,        FAN17500
     2 2, LDUM    ,NDIME,.TRUE.,PRPREL)                                 FAN17600
      FAILED=FAILED .OR. VAR.GT.VARE                                    FAN17700
      IF (VAR .GT. VARE) GO TO 300                                      FAN17800
      VARE=VAR                                                          FAN17900
      DO 605 J=1,JLAM                                                   FAN18000
      LAMST(J)=SNGL(PLAM(J))                                            FAN18100
  605 LAMF(J,JLAM,1)=LAMST(J)                                           FAN18200
      SAVEPK=VARTRY .GE. 1.D+20                                         FAN18300
      IF (FAILED) GO TO 610                                             FAN18400
 5600 FORMAT (/51H END OF FIT TO RAW DATA, START OF FIT TO TRANSFORMS)  FAN18500
      IF (PRPREL) WRITE (IWRITE,5600)                                   FAN18600
      CALL LSTSQR (LAMST,JLAM,.FALSE.,NF,Y,T,JLAM+IBASE,                FAN18700
     1 VARF,GSE,NDIMG,SQRTW,YLYFIT,E,NMAX,TSTART,DELTAT,NT,NINT,        FAN18800
     2 3,.FALSE.,NDIME,.TRUE.,PRPREL)                                   FAN18900
      IF (FAILED .AND. VAR.GE.VARTRY) GO TO 300                         FAN19000
      SAVEPK=.TRUE.                                                     FAN19100
      VARTRY=VAR                                                        FAN19200
      DO 608 J=1,JLAM                                                   FAN19300
  608 LAMF(J,JLAM,2)=SNGL(PLAM(J))                                      FAN19400
  610 IF (.NOT.SAVEPK .AND. JLAM.NE.1) GO TO 300                        FAN19500
      IF (JLAM .EQ. NLAMMX) GO TO 670                                   FAN19600
C-----------------------------------------------------------------------FAN19700
C  FINDS MAX SUM OF ABSOLUTE ADJACENT AMPLITUDES FOR LATER DETERMINATIONFAN19800
C      OF STARTING VALUE FOR LAMST(JLAMP1)                              FAN19900
C-----------------------------------------------------------------------FAN20000
      DUB=DABS(PALPHA(1))+DABS(PALPHA(JLAMP1))                          FAN20100
      LLPEAK=1                                                          FAN20200
      IF (DUB .GT. DABS(PALPHA(JLAM))) GO TO 620                        FAN20300
      LLPEAK=JLAMP1                                                     FAN20400
      DUB=DABS(PALPHA(JLAM))                                            FAN20500
  620 IF (JLAM .EQ. 1) GO TO 690                                        FAN20600
      DO 630 J=2,JLAM                                                   FAN20700
      DDUB=DABS(PALPHA(J-1))+DABS(PALPHA(J))                            FAN20800
      IF (DDUB .LE. DUB) GO TO 630                                      FAN20900
      DUB=DDUB                                                          FAN21000
      LLPEAK=J                                                          FAN21100
  630 CONTINUE                                                          FAN21200
  670 IF (.NOT.FAILED .OR. JLAM.EQ.1) GO TO 690                         FAN21300
  300 CONTINUE                                                          FAN21400
      JTRY=MTRY                                                         FAN21500
  690 LPEAK=LLPEAK                                                      FAN21600
      NFTRY(JLAM)=JTRY                                                  FAN21700
      FFAIL(JLAM,2)=FAILED                                              FAN21800
      FFAIL(JLAM,1)=VARTRY .GE. 1.D+20                                  FAN21900
      IF (.NOT.FFAIL(JLAM,1)) GO TO 200                                 FAN22000
      DO 695 J=1,JLAM                                                   FAN22100
  695 LAMF(J,JLAM,2)=LAMF(J,JLAM,1)                                     FAN22200
  200 CONTINUE                                                          FAN22300
      RETURN                                                            FAN22400
      END                                                               FAN22500
C--------------------- DOUBLE PRECISION VERSION 2B ---------------------YAN00100
C  SUBROUTINE YANLYZ.  FOR CONSTRAINED STEPWISE LEAST SQUARES ANALYSIS  YAN00200
C      OF RAW DATA USING LAMF (STARTING VALUES OBTAINED FROM FANLYZ).   YAN00300
C-----------------------------------------------------------------------YAN00400
C  CALLS SUBPROGRAMS - LSTSQR, FISHER, RESIDU                           YAN00500
C  WHICH IN TURN CALL - EVAR, ETHEOR, PIVOT, PIVOT1, PLPRES             YAN00600
C-----------------------------------------------------------------------YAN00700
      SUBROUTINE YANLYZ (T,Y,SQRTW,YLYFIT,E,NMAX,NDIME,                 YAN00800
     1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG)                               YAN00900
      EXTERNAL EVAR                                                     YAN01000
      LOGICAL NOBASE, PRFINL, PLOTRS, REPEAT, WTED, FAILED, YFAIL, DONE,YAN01100
     1 REGINT, FFAIL, LBLOKA, LBLOKB, LBLOKC                            YAN01200
      REAL LAMNMX, LAMST, LAMF                                          YAN01300
      DOUBLE PRECISION VARBES, VARMIN, VAR, PALPHA, PLAM, YLYFIT, E,    YAN01400
     1 TSTART, DELTAT, GSE, DBLOKA, DBLOKB, T                           YAN01500
      COMMON /BLOK/ DBLOKA(4843), VAR, PALPHA(10), PLAM(10),            YAN01600
     1 DBLOKB(1069), RBLOKA(3), LAMST(9), LAMF(9,9,2), RBLOKB,          YAN01700
     2 LAMNMX(2,2), SIGYY, ENPHI, SIGMAP(19), PCERR(19), RBLOKC(252),   YAN01800
     3 IBLOKA(200), IBASE, IBLOKB(12), NLAMMX, IWT, LABEL(80), NINT,    YAN01900
     4 NF, N, IWRITE, IBLOKC, LINEPG, IBLOKD(3), ITER, JLAMP1,          YAN02000
     5 NFTRY(9), REGINT, FAILED, LBLOKA, NOBASE, WTED, REPEAT           YAN02100
      COMMON /BLOK/ PLOTRS, PRFINL, LBLOKB, FFAIL(9,2), LBLOKC(19)      YAN02200
      DIMENSION E(NDIME,10), YFAIL(9), ASAVE(10,9,6),                   YAN02300
     1 PROBSV(9), NEXBES(6,5), VARSV(9), SIGYSV(9), ITERSV(9),          YAN02400
     2 ENPHSV(9), PROBLN(9), DONE(9)                                    YAN02500
      DIMENSION T(NMAX), Y(NMAX), SQRTW(NMAX), GSE(NDIMG,4),            YAN02600
     1 YLYFIT(NMAX), TSTART(NINTMX), DELTAT(NINTMX),                    YAN02700
     2 NT(NINTMX), PUNCOR(6)                                            YAN02800
      DATA NEXBES/1H , 1H , 1H , 1H , 1H , 1H , 1HS, 1HE, 1HC, 1HO,     YAN02900
     1 1HN, 1HD, 1H , 1HT, 1HH, 1HI, 1HR, 1HD, 1HF, 1HO, 1HU, 1HR, 1HT, YAN03000
     2 1HH, 1H , 1HF, 1HI, 1HF, 1HT, 1HH/                               YAN03100
      WTED=IWT .NE. 1                                                   YAN03200
      JLAMBS=1                                                          YAN03300
      VARBES=1.D+20                                                     YAN03400
      VARMIN=1.D+20                                                     YAN03500
      ABSYMX=0.                                                         YAN03600
      DO 180 J=1,N                                                      YAN03700
  180 ABSYMX=AMAX1(ABSYMX,ABS(Y(J)))                                    YAN03800
C-----------------------------------------------------------------------YAN03900
C  START OF MAIN LOOP FOR STEPWISE ANALYSIS                             YAN04000
C-----------------------------------------------------------------------YAN04100
      DO 200 JLAM=1,NLAMMX                                              YAN04200
      JLAMP1=JLAM+1                                                     YAN04300
      NF=2*JLAM+IBASE                                                   YAN04400
      IF (.NOT.PRFINL) GO TO 210                                        YAN04500
 5200 FORMAT (1H1,25X,80A1/)                                            YAN04600
      WRITE (IWRITE,5200) LABEL                                         YAN04700
 5240 FORMAT (40X,23HFINAL ANALYSIS ASSUMING,I2,11H COMPONENTS)         YAN04800
      WRITE (IWRITE,5240) JLAM                                          YAN04900
  210 DO 220 J=1,JLAM                                                   YAN05000
  220 LAMST(J)=LAMF(J,JLAM,2)                                           YAN05100
      CALL LSTSQR (LAMST,JLAM,.TRUE.,N,Y,T,JLAM+IBASE,EVAR,             YAN05200
     1 GSE,NDIMG,SQRTW,YLYFIT,E,NMAX,TSTART,DELTAT,NT,NINT,             YAN05300
     2 4,FFAIL(JLAM,2),NDIME,.FALSE.,PRFINL)                            YAN05400
      YFAIL(JLAM)=FAILED                                                YAN05500
      IF (.NOT.FAILED) GO TO 225                                        YAN05600
      IF (FFAIL(JLAM,1)) GO TO 250                                      YAN05700
 5222 FORMAT (////43H SECOND ATTEMPT USING STARTING LAMBDAS FROM,       YAN05800
     1 29H INITIAL ANALYSIS OF RAW DATA)                                YAN05900
      IF (PRFINL) WRITE (IWRITE,5222)                                   YAN06000
      DO 222 J=1,JLAM                                                   YAN06100
  222 LAMST(J)=LAMF(J,JLAM,1)                                           YAN06200
      CALL LSTSQR (LAMST,JLAM,.TRUE.,N,Y,T,JLAM+IBASE,EVAR,             YAN06300
     1 GSE,NDIMG,SQRTW,YLYFIT,E,NMAX,TSTART,DELTAT,NT,NINT,             YAN06400
     2 4,.FALSE.,NDIME,.FALSE.,PRFINL)                                  YAN06500
      YFAIL(JLAM)=FAILED                                                YAN06600
      IF (FAILED) GO TO 250                                             YAN06700
C-----------------------------------------------------------------------YAN06800
C  SAVE VALUES FOR LATER USE AND FOR FINAL OUTPUT                       YAN06900
C-----------------------------------------------------------------------YAN07000
  225 SIGYSV(JLAM)=SIGYY                                                YAN07100
      VARSV(JLAM)=SNGL(VAR)                                             YAN07200
      ITERSV(JLAM)=ITER                                                 YAN07300
      ENPHSV(JLAM)=ENPHI                                                YAN07400
      I=JLAM                                                            YAN07500
      DO 230 K=1,JLAM                                                   YAN07600
      ASAVE(K,JLAM,1)=SNGL(PALPHA(K))                                   YAN07700
      ASAVE(K,JLAM,2)=SNGL(PLAM(K))                                     YAN07800
      ASAVE(K,JLAM,3)=SIGMAP(K)                                         YAN07900
      I=I+1                                                             YAN08000
      ASAVE(K,JLAM,4)=SIGMAP(I)                                         YAN08100
      ASAVE(K,JLAM,5)=PCERR(K)                                          YAN08200
      ASAVE(K,JLAM,6)=PCERR(I)                                          YAN08300
  230 CONTINUE                                                          YAN08400
      IF (NOBASE) GO TO 240                                             YAN08500
      ASAVE(JLAMP1,JLAM,2)=0.                                           YAN08600
      ASAVE(JLAMP1,JLAM,3)=0.                                           YAN08700
      ASAVE(JLAMP1,JLAM,5)=0.                                           YAN08800
      ASAVE(JLAMP1,JLAM,1)=SNGL(PALPHA(JLAMP1))                         YAN08900
      ASAVE(JLAMP1,JLAM,4)=SIGMAP(NF)                                   YAN09000
      ASAVE(JLAMP1,JLAM,6)=PCERR(NF)                                    YAN09100
C-----------------------------------------------------------------------YAN09200
C  TEST IF THIS SOLUTION IS THE BEST SO FAR                             YAN09300
C-----------------------------------------------------------------------YAN09400
  240 IF (VARBES .GE. 1.D+20) GO TO 248                                 YAN09500
      IF (VAR .GE. VARMIN) GO TO 200                                    YAN09600
      K=N-NF                                                            YAN09700
      DUM=SNGL(VARBES/VAR-1.D0)                                         YAN09800
      L=2*(JLAM-JLAMBS)                                                 YAN09900
      IF (FISHER(FLOAT(K)*DUM/(FLOAT(L)*(1.+ENPHI*FLOAT((NF+2)*N)/      YAN10000
     1 FLOAT(NF*K))),L,K) .LE. .5) GO TO 250                            YAN10100
  248 JLAMBS=JLAM                                                       YAN10200
      VARBES=VAR                                                        YAN10300
  250 VARMIN=DMIN1(VARMIN,VAR)                                          YAN10400
  200 CONTINUE                                                          YAN10500
C-----------------------------------------------------------------------YAN10600
C  END OF MAIN LOOP FOR STEPWISE ANALYSIS                               YAN10700
C-----------------------------------------------------------------------YAN10800
      IF (VARBES .LT. 1.D+20) GO TO 300                                 YAN10900
 5300 FORMAT (///50H ALL ANALYSES OF THE RAW DATA FAILED - CHECK INPUT, YAN11000
     1 22H DATA FOR GROSS ERRORS   )                                    YAN11100
      WRITE (IWRITE,5300)                                               YAN11200
      RETURN                                                            YAN11300
C-----------------------------------------------------------------------YAN11400
C  CALCULATE PROBSV AND PROBLN (PROBABILITY AND UNCORRECTED PROBABILITY YAN11500
C      (I.E., USING ENPHI=0.) THAT THE BEST SOLUTION IS ACTUALLY        YAN11600
C      BETTER THAN THE OTHERS)                                          YAN11700
C-----------------------------------------------------------------------YAN11800
  300 PNXBES=1.E+20                                                     YAN11900
      JLAMNX=JLAMBS                                                     YAN12000
      DO 310 J=1,NLAMMX                                                 YAN12100
      PROBSV(J)=2.                                                      YAN12200
      IF (J.EQ.JLAMBS .OR. YFAIL(J)) GO TO 310                          YAN12300
      I=2*MAX0(JLAMBS,J)+IBASE                                          YAN12400
      K=N-I                                                             YAN12500
      L=2*(JLAMBS-J)                                                    YAN12600
      DUM=VARSV(J)/SNGL(VARBES)                                         YAN12700
      IF (L .LT. 0) GO TO 315                                           YAN12800
      DDUM=FLOAT(K)*(DUM-1.)/FLOAT(L)                                   YAN12900
      DDDUM=1.+ENPHSV(JLAMBS)*FLOAT((I+2)*N)/FLOAT(I*K)                 YAN13000
      PROBSV(J)=FISHER(DDUM/DDDUM,L,K)                                  YAN13100
      PROBLN(J)=FISHER(DDUM,L,K)                                        YAN13200
      GO TO 325                                                         YAN13300
  315 IF (DUM .LT. 1.) GO TO 320                                        YAN13400
      PROBSV(J)=1.                                                      YAN13500
      PROBLN(J)=1.                                                      YAN13600
      GO TO 325                                                         YAN13700
  320 DDUM=FLOAT(K)*(1.-1./DUM)/FLOAT(L)                                YAN13800
      DDDUM=1.+ENPHSV(J)*FLOAT((I+2)*N)/FLOAT(I*K)                      YAN13900
      PROBLN(J)=1.-FISHER(DDUM,-L,K)                                    YAN14000
      PROBSV(J)=1.-FISHER(DDUM/DDDUM,-L,K)                              YAN14100
  325 IF (PROBSV(J) .GE. PNXBES) GO TO 310                              YAN14200
      PNXBES=PROBSV(J)                                                  YAN14300
      PNXLIN=PROBLN(J)                                                  YAN14400
      JLAMNX=J                                                          YAN14500
  310 CONTINUE                                                          YAN14600
      PROBSV(JLAMBS)=0.                                                 YAN14700
      IF (PNXBES .LT. 1.E+20) GO TO 330                                 YAN14800
      PNXBES=1.                                                         YAN14900
      PNXLIN=1.                                                         YAN15000
C-----------------------------------------------------------------------YAN15100
C  FINAL SUMMARY OF RESULTS - SUMMARIZES UP TO BEST 5 SOLUTIONS.        YAN15200
C-----------------------------------------------------------------------YAN15300
  330 CALL RESIDU (JLAMBS,PLOTRS,T,Y,SQRTW,YLYFIT,NMAX,PUNCOR,ASAVE)    YAN15400
 5330 FORMAT (38H1DISCRETE - Version 2B (December 1990),10X,80A1/)      YAN15500
  340 WRITE (IWRITE,5330) LABEL                                         YAN15600
      DO 350 J=1,NLAMMX                                                 YAN15700
  350 DONE(J)=YFAIL(J)                                                  YAN15800
      L=17+JLAMBS+IBASE                                                 YAN15900
 5350 FORMAT (30X,6A1,14H BEST SOLUTION,I13,11H COMPONENTS)             YAN16000
      WRITE (IWRITE,5350) (NEXBES(J,1),J=1,6),JLAMBS                    YAN16100
      IF (PNXBES .GT. .95) GO TO 360                                    YAN16200
 5360 FORMAT (48H ++++++ NOTE - THE APPROX. PROBABILITY THAT THIS,      YAN16300
     1 34H SOLUTION IS ACTUALLY BEST IS ONLY,F6.3,                      YAN16400
     2 38H.  LOOK AT SECOND BEST SOLUTION ALSO. ,6(1H+))                YAN16500
      WRITE (IWRITE,5360) PNXBES                                        YAN16600
      L=L-1                                                             YAN16700
  360 DO 370 J=1,NLAMMX                                                 YAN16800
      IF (J .GT. 1) GO TO 375                                           YAN16900
      DONE(JLAMBS)=.TRUE.                                               YAN17000
      I=JLAMBS                                                          YAN17100
      GO TO 400                                                         YAN17200
  375 DUM=1.E+20                                                        YAN17300
      DO 380 K=1,NLAMMX                                                 YAN17400
      IF (PROBSV(K).GT.DUM .OR. DONE(K)) GO TO 380                      YAN17500
      DUM=PROBSV(K)                                                     YAN17600
      I=K                                                               YAN17700
  380 CONTINUE                                                          YAN17800
      IF (DUM.GE.1.E+20 .OR. J.GT.5) GO TO 500                          YAN17900
      DONE(I)=.TRUE.                                                    YAN18000
      CALL RESIDU (I,.FALSE.,T,Y,SQRTW,YLYFIT,NMAX,PUNCOR,ASAVE)        YAN18100
      L=L+15+I+IBASE                                                    YAN18200
      IF (L-2 .LE. LINEPG) GO TO 385                                    YAN18300
      L=14+I+IBASE                                                      YAN18400
      WRITE (IWRITE,5999)                                               YAN18500
 5999 FORMAT (1H1)                                                      YAN18600
      GO TO 390                                                         YAN18700
  385 WRITE (IWRITE,5001)                                               YAN18800
      WRITE (IWRITE,5001)                                               YAN18900
 5001 FORMAT (1X)                                                       YAN19000
  390 IF (DUM .GT. .95) WRITE (IWRITE,5350) (NEXBES(K,J),K=1,6),I       YAN19100
      IF (DUM .LE. .95) WRITE (IWRITE,5385) (NEXBES(K,J),K=1,6),I       YAN19200
 5385 FORMAT (30X,6A1,14H BEST SOLUTION,I13,11H COMPONENTS,3X,          YAN19300
     1 30H-   A SIGNIFICANT POSSIBILITY ,25(1H+))                       YAN19400
  400 WRITE (IWRITE,5001)                                               YAN19500
      WRITE (IWRITE,5400)                                               YAN19600
 5400 FORMAT (1X,70(1H*))                                               YAN19700
      WRITE (IWRITE,5402) NFTRY(I)                                      YAN19800
 5402 FORMAT (2H *,5X,25HALPHA +- STD ERR  PERCENT,11X,                 YAN19900
     1 31HLAMBDA +- STD ERR  PERCENT *   ,                              YAN20000
     2 41HSTARTING LAMBDA (FROM FIT TO TRANSFORMS -,                    YAN20100
     3 I4,7H TRIES))                                                    YAN20200
 5404 FORMAT (2H *,1PE10.3,3H +-,E8.1,0PF9.3,1PE17.3,3H +-,             YAN20300
     1 E8.1,0PF9.3,2H *,1PE18.3,                                        YAN20400
     2 39H (NO EXACT FIT TO THE TRANSFORMS FOUND))                      YAN20500
      IF (FFAIL(I,2)) WRITE (IWRITE,5404) (ASAVE(K,I,1),ASAVE(K,I,4),   YAN20600
     1 ASAVE(K,I,6),ASAVE(K,I,2),ASAVE(K,I,3),ASAVE(K,I,5),LAMF(K,I,2), YAN20700
     2 K=1,I)                                                           YAN20800
 5405 FORMAT (2H *,1PE10.3,3H +-,E8.1,0PF9.3,1PE17.3,3H +-,             YAN20900
     1 E8.1,0PF9.3,2H *,1PE18.3)                                        YAN21000
      IF (.NOT.FFAIL(I,2)) WRITE (IWRITE,5405) (ASAVE(K,I,1),           YAN21100
     1 ASAVE(K,I,4),ASAVE(K,I,6),ASAVE(K,I,2),ASAVE(K,I,3),ASAVE(K,I,5),YAN21200
     2 LAMF(K,I,2),K=1,I)                                               YAN21300
      IF (NOBASE) GO TO 410                                             YAN21400
      K=I+1                                                             YAN21500
      WRITE (IWRITE,5404) ASAVE(K,I,1),ASAVE(K,I,4),ASAVE(K,I,6),       YAN21600
     1 ASAVE(K,I,2),ASAVE(K,I,3),ASAVE(K,I,5)                           YAN21700
  410 WRITE (IWRITE,5400)                                               YAN21800
      IF (J.EQ.1 .AND. PNXBES.GT..95) WRITE (IWRITE,5410) JLAMNX,       YAN21900
     1 JLAMBS, PNXBES                                                   YAN22000
 5410 FORMAT (93X,20(1H*)/34H APPROXIMATE PROBABILITY THAT THIS,        YAN22100
     1 56H SOLUTION IS REALLY BETTER THAN THE SECOND BEST SOLUTION,     YAN22200
     2 9H = * PNG(,I1,1H/,I1,3H) =,F6.3,2H */93X,20(1H*))               YAN22300
 5412 FORMAT (/5H PNG(,I1,1H/,I1,1H),12X,1H=,F21.3,13X,4HNPHI,13X,      YAN22400
     1 1H=,1PE11.2,13X,25H(UNCORRECTED PNG WOULD BE,0PF8.3,1H))         YAN22500
      IF (J .GT. 1) WRITE (IWRITE,5412) I,JLAMBS,PROBSV(I),ENPHSV(I),   YAN22600
     1 PROBLN(I)                                                        YAN22700
 5414 FORMAT (22H ITERATIONS IN FIT   =,I21,13X,                        YAN22800
     1 18HSTD. DEV. OF FIT =,1PE11.4,13X,                               YAN22900
     2 27HSIGNAL/NOISE RATIO OF FIT =,0PF7.0)                           YAN23000
      DDUM=ABSYMX/SIGYSV(I)                                             YAN23100
      IF (WTED) WRITE (IWRITE,5414) ITERSV(I),SIGYSV(I)                 YAN23200
      IF (.NOT.WTED) WRITE (IWRITE,5414) ITERSV(I),SIGYSV(I),DDUM       YAN23300
 5416 FORMAT (21H LAMBDA HELD BETEWEEN,1PE9.2,4H AND,E9.2,13X,          YAN23400
     1 4HNPHI,13X,1H=,E11.2)                                            YAN23500
      IF (J .EQ. 1) WRITE (IWRITE,5416) LAMNMX(1,1),LAMNMX(2,1),        YAN23600
     1 ENPHSV(JLAMBS)                                                   YAN23700
 5418 FORMAT (/4H LAG,25X,5(4X,3HK =,I2),23X,16(1H*)/                   YAN23800
     1 29H PROB. RESIDUALS UNCORRELATED,5F9.3,4X,                       YAN23900
     2 27HWEIGHTED AVERAGE = * PUNC =,F6.3,2H *)                        YAN24000
      WRITE (IWRITE,5418) (K,K=1,5),PUNCOR                              YAN24100
 5422 FORMAT (1X,131(1H=))                                              YAN24200
      IF (L-1 .LE. LINEPG) WRITE (IWRITE,5422)                          YAN24300
      IF (L .LE. LINEPG) WRITE (IWRITE,5422)                            YAN24400
  370 CONTINUE                                                          YAN24500
  500 IF (.NOT.REPEAT) RETURN                                           YAN24600
      REPEAT=.FALSE.                                                    YAN24700
      CALL RESIDU (JLAMBS,.FALSE.,T,Y,SQRTW,YLYFIT,NMAX,PUNCOR,ASAVE)   YAN24800
      GO TO 340                                                         YAN24900
      END                                                               YAN25000
C--------------------- DOUBLE PRECISION VERSION 2B ---------------------LSQ00100
C  SUBROUTINE LSTSQR.  CONSTRAINED LEAST SQUARES FIT USING STEPWISE     LSQ00200
C      REGRESSION APPROACH.                                             LSQ00300
C  N = NO. OF DATA POINTS TO BE FIT, I.E., =N (RAWDAT=.TRUE.) OR        LSQ00400
C      =2*JLAM+IBASE (RAWDAT=.FALSE.).                                  LSQ00500
C  Y,T = RAW DATA, REGARDLESS OF RAWDAT.                                LSQ00600
C  ITYPE=1,2,3,4 FOR PRELIMINARY DETERMINATION OF WEIGHTS, FIT TO RAW   LSQ00700
C      DATA WITH INTERPOLATION, FIT TO TRANSFORMS, AND FINAL FIT TO     LSQ00800
C      RAW DATA, RESPECTIVELY.                                          LSQ00900
C-----------------------------------------------------------------------LSQ01000
C  CALLS SUBPROGRAMS - VARDUM (A DUMMY NAME), PIVOT, ANLERR             LSQ01100
C  WHICH IN TURN CALL - PIVOT1, ETHEOR                                  LSQ01200
C-----------------------------------------------------------------------LSQ01300
      SUBROUTINE LSTSQR (LAMST,JLAM,RAWDAT,N,Y,T,NLIN,VARDUM,           LSQ01400
     1 GSE,NDIMG,SQRTW,YLYFIT,E,NMAX,TSTART,DELTAT,NT,NINT,             LSQ01500
     2 ITYPE,SPLIT,NDIME,INTER,PRLSQ)                                   LSQ01600
      EXTERNAL VARDUM                                                   LSQ01700
      REAL LAMNMX, LAMST                                                LSQ01800
      LOGICAL RAWDAT, NOBASE, PRLSQ, FAILED, ALLPIV(2), PIVLAM,         LSQ01900
     1 LDUM, INTER, PIVALP, LBLOKA, LBLOKB, LBLOKC, SPLIT               LSQ02000
      DOUBLE PRECISION GSE, YLYFIT, E, TSTART, DELTAT, VARMIN, VAR,     LSQ02100
     1 PALPHA, PLAM, PLMTRY, DELTAP, DFCAPL, SFDE, ADUB, DDSE,          LSQ02200
     2 DTOLER, F, DF, VAROLD, VARG, R, DFLAT, DALPL, DUB, DDUB, VT,     LSQ02300
     3 DBLOKA, DBLOKB, DBLOKC, T                                        LSQ02400
      COMMON /BLOK/ DBLOKA(4828), VARMIN, DBLOKB(14), VAR, PALPHA(10),  LSQ02500
     1 PLAM(10), DBLOKC(100), PLMTRY(10), DELTAP(19), DFCAPL(10,9),     LSQ02600
     2 SFDE(9), ADUB(19,19), DDSE(9,9), DTOLER(19), F(19,10),           LSQ02700
     3 DF(19,10), DELTA, FHATMX, CONVRG, RBLOKA(171), PRECIS,           LSQ02800
     4 LAMNMX(2,2), RBLOKB(2), SIGMAP(19), PCERR(19), RBLOKC(252),      LSQ02900
     5 IBLOKA(199), MDIMA, IBLOKB(96), NF, IBLOKC, IWRITE, IBLOKD(5)    LSQ03000
      COMMON /BLOK/ ITER, JLAMP1, IBLOKE(9), LBLOKA, FAILED, LBLOKB,    LSQ03100
     1 NOBASE, LBLOKC(23), PIVALP(19)                                   LSQ03200
      DIMENSION PIVLAM(9), R(10), DALPL(10,9)                           LSQ03300
      DIMENSION LAMST(JLAM), Y(NMAX), T(NMAX),                          LSQ03400
     1 HOLRTH(2), SQRTW(NMAX), YLYFIT(NMAX),                            LSQ03500
     2 E(NDIME,NLIN), TSTART(NINT), QT(3), VT(3), GSE(NDIMG,4),         LSQ03600
     3 DELTAT(NINT), NT(NINT), MXITER(4), RLAMMN(2)                     LSQ03700
      DATA QMIN/.001/, HOLRTH/1H , 1H*/, MCONV/3/, NABORT/4/,           LSQ03800
     1 MXITER/30, 30, 40, 40/,RLAMMN/1.1, 2./,DFLAT/-1.D-2/             LSQ03900
      ISPLIT=1                                                          LSQ04000
      IF (SPLIT) ISPLIT=2                                               LSQ04100
      JLAMP1=JLAM+1                                                     LSQ04200
      NCONV=0                                                           LSQ04300
      NVARUP=0                                                          LSQ04400
      NFLAT=0                                                           LSQ04500
      ALLPIV(2)=.TRUE.                                                  LSQ04600
      IF (NOBASE) PALPHA(JLAMP1)=0.D0                                   LSQ04700
      PLAM(JLAMP1)=0.D0                                                 LSQ04800
      PLMTRY(JLAMP1)=0.D0                                               LSQ04900
      PIVALP(JLAMP1)=.TRUE.                                             LSQ05000
C-----------------------------------------------------------------------LSQ05100
C  PUT PLAM IN ASCENDING ORDER AND SEPARATE THEM IF NECESSARY           LSQ05200
C-----------------------------------------------------------------------LSQ05300
      DO 150 J=1,JLAM                                                   LSQ05400
      PIVLAM(J)=.FALSE.                                                 LSQ05500
  150 DELTAP(J)=0.D0                                                    LSQ05600
      DO 160 J=1,JLAM                                                   LSQ05700
      DUM=1.E+20                                                        LSQ05800
      DO 165 K=1,JLAM                                                   LSQ05900
      IF (PIVLAM(K) .OR. LAMST(K).GT.DUM) GO TO 165                     LSQ06000
      L=K                                                               LSQ06100
      DUM=LAMST(K)                                                      LSQ06200
  165 CONTINUE                                                          LSQ06300
      PIVLAM(L)=.TRUE.                                                  LSQ06400
      IF (J .EQ. 1) PLAM(J)=DBLE(AMAX1(DUM,LAMNMX(1,ISPLIT)))           LSQ06500
      IF (J .GT. 1) PLAM(J)=DBLE(AMIN1(AMAX1(DUM,RLAMMN(ISPLIT)*        LSQ06600
     1 SNGL(PLAM(J-1))),LAMNMX(2,1)))                                   LSQ06700
  160 CONTINUE                                                          LSQ06800
      VAROLD=1.D+20                                                     LSQ06900
      VARG=1.D+20                                                       LSQ07000
      ITER=-1                                                           LSQ07100
      Q=0.                                                              LSQ07200
      IF (PRLSQ) WRITE (IWRITE,5205)                                    LSQ07300
 5205 FORMAT (/4H ITR,2X,8HVARIANCE,2X,9HDAMPING Q,3X,8HBASELINE,       LSQ07400
     1 5(5X,5HALPHA,3X,6HLAMBDA))                                       LSQ07500
C-----------------------------------------------------------------------LSQ07600
C  START OF MAIN LOOP FOR THE LEAST SQUARES FIT                         LSQ07700
C-----------------------------------------------------------------------LSQ07800
  210 NTRY=1                                                            LSQ07900
  215 ITER=ITER+1                                                       LSQ08000
  220 CALL VARDUM (Q,JLAM,VAR,T,SQRTW,NMAX,TSTART,DELTAT,NT,            LSQ08100
     1 NINT,.TRUE.,LDUM,E,Y,YLYFIT,NLIN,VARG,NTRY,NDIME,INTER,          LSQ08200
     2 GSE,NDIMG,IERROR)                                                LSQ08300
      IF (IERROR .EQ. 5) GO TO 800                                      LSQ08400
      ALLPIV(1)=LDUM                                                    LSQ08500
      IF (VAR.LE.VARG .OR. NTRY.EQ.2) GO TO 240                         LSQ08600
      NTRY=2                                                            LSQ08700
      Q=QG                                                              LSQ08800
      GO TO 220                                                         LSQ08900
  240 DO 242 J=1,JLAM                                                   LSQ09000
      R(J)=PALPHA(J)                                                    LSQ09100
      IF (.NOT.RAWDAT) R(J)=PALPHA(J)*EXP(-SNGL(PLMTRY(J))*DELTA)       LSQ09200
  242 PLAM(J)=PLMTRY(J)                                                 LSQ09300
      DUB=VAR-VAROLD                                                    LSQ09400
      IF (.NOT.PRLSQ) GO TO 250                                         LSQ09500
      NU=MIN0(5,JLAM)                                                   LSQ09600
      DDUM=HOLRTH(1)                                                    LSQ09700
      IF (DUB .GT. 0.D0) DDUM=HOLRTH(2)                                 LSQ09800
      DO 244 J=1,JLAM                                                   LSQ09900
      PCERR(J)=HOLRTH(1)                                                LSQ10000
      IF (.NOT.PIVLAM(J)) PCERR(J)=HOLRTH(2)                            LSQ10100
  244 CONTINUE                                                          LSQ10200
 5240 FORMAT (1X,I2,1PD11.4,A1,E10.2,A1,D10.2,1X,5(1P2D9.2,A1))         LSQ10300
      WRITE (IWRITE,5240) ITER,VAR,DDUM,Q,HOLRTH(NTRY),PALPHA(JLAMP1),  LSQ10400
     1 (R(J),PLAM(J),PCERR(J),J=1,NU)                                   LSQ10500
      IF (JLAM .LE. 5) GO TO 250                                        LSQ10600
 5250 FORMAT (37X,4(1P2D9.2,A1))                                        LSQ10700
      WRITE (IWRITE,5250) (R(J),PLAM(J),PCERR(J),J=6,JLAM)              LSQ10800
 5001 FORMAT (1X)                                                       LSQ10900
      WRITE (IWRITE,5001)                                               LSQ11000
C-----------------------------------------------------------------------LSQ11100
C  TESTS FOR CONVERGENCE AND DIVERGENCE                                 LSQ11200
C-----------------------------------------------------------------------LSQ11300
  250 IF ((DABS(DUB).GT.DMAX1(VARMIN,VAROLD*CONVRG) .AND. RAWDAT) .OR.  LSQ11400
     1 (VAR.GT.DBLE(FLOAT(NF)*(CONVRG*FHATMX)**2) .AND. .NOT.RAWDAT))   LSQ11500
     2 GO TO 270                                                        LSQ11600
      NCONV=NCONV+1                                                     LSQ11700
      IF (NCONV.GE.MCONV .OR. .NOT.RAWDAT) GO TO 800                    LSQ11800
      NVARUP=0                                                          LSQ11900
      GO TO 300                                                         LSQ12000
  270 NCONV=0                                                           LSQ12100
      IF (DUB .LE. 0.D0) GO TO 275                                      LSQ12200
      NVARUP=NVARUP+1                                                   LSQ12300
      IF (NVARUP.LT.NABORT .OR. ITER.LT.JLAM) GO TO 280                 LSQ12400
      IERROR=1                                                          LSQ12500
      GO TO 800                                                         LSQ12600
  275 NVARUP=0                                                          LSQ12700
  280 IF (ITER.LT.MXITER(ITYPE)) GO TO 290                              LSQ12800
      IERROR=2                                                          LSQ12900
      GO TO 800                                                         LSQ13000
  290 IF (DUB/VAROLD.LT.DFLAT .OR. RAWDAT) GO TO 295                    LSQ13100
      NFLAT=NFLAT+1                                                     LSQ13200
      IF (NFLAT .LT. NABORT) GO TO 300                                  LSQ13300
      IERROR=3                                                          LSQ13400
      GO TO 800                                                         LSQ13500
  295 NFLAT=0                                                           LSQ13600
C-----------------------------------------------------------------------LSQ13700
C  FINISH COMPUTATION OF DFCAPL (PARTIAL OF FCAP W.R.T. LAMBDA)         LSQ13800
C-----------------------------------------------------------------------LSQ13900
  300 DO 310 K=1,JLAM                                                   LSQ14000
      DO 320 J=1,NLIN                                                   LSQ14100
      DFCAPL(J,K)=DFCAPL(J,K)*PALPHA(K)                                 LSQ14200
  320 CONTINUE                                                          LSQ14300
      DFCAPL(K,K)=DFCAPL(K,K)+SFDE(K)                                   LSQ14400
  310 CONTINUE                                                          LSQ14500
C-----------------------------------------------------------------------LSQ14600
C  CALCULATE DALPL (PARTIAL OF PALPHA W.R.T. LAMBDA) FROM MATRIX        LSQ14700
C  PRODUCT ADUB*DFCAPL                                                  LSQ14800
C-----------------------------------------------------------------------LSQ14900
      DO 400 J=1,NLIN                                                   LSQ15000
      DO 400 K=1,JLAM                                                   LSQ15100
  400 DALPL(J,K)=0.D0                                                   LSQ15200
      DO 410 K=1,JLAM                                                   LSQ15300
      IF (.NOT.PIVALP(K)) GO TO 410                                     LSQ15400
      DO 420 J=1,NLIN                                                   LSQ15500
      IF (.NOT.PIVALP(J)) GO TO 420                                     LSQ15600
      DUB=0.D0                                                          LSQ15700
      DO 430 I=1,NLIN                                                   LSQ15800
      IF (PIVALP(I)) DUB=DUB+ADUB(J,I)*DFCAPL(I,K)                      LSQ15900
  430 CONTINUE                                                          LSQ16000
      DALPL(J,K)=DUB                                                    LSQ16100
  420 CONTINUE                                                          LSQ16200
  410 CONTINUE                                                          LSQ16300
C-----------------------------------------------------------------------LSQ16400
C  CALCULATE ADUB (NORMAL EQUATION MATRIX FOR NONLINEAR LEAST SQUARES)  LSQ16500
C  AND INVERT IT                                                        LSQ16600
C-----------------------------------------------------------------------LSQ16700
      IF (.NOT.RAWDAT) GO TO 550                                        LSQ16800
      DO 510 K=1,JLAM                                                   LSQ16900
      DO 520 J=1,K                                                      LSQ17000
      DUB=0.D0                                                          LSQ17100
      DO 530 L=1,NLIN                                                   LSQ17200
  530 DUB=DUB-DALPL(L,J)*DFCAPL(L,K)                                    LSQ17300
      IF (DABS(PALPHA(K)) .GT. DABS(PALPHA(J))) GO TO 540               LSQ17400
      DUB=DUB+DALPL(J,K)*SFDE(J)                                        LSQ17500
      DUB=DUB+DALPL(K,J)*SFDE(K)                                        LSQ17600
      GO TO 545                                                         LSQ17700
  540 DUB=DUB+DALPL(K,J)*SFDE(K)                                        LSQ17800
      DUB=DUB+DALPL(J,K)*SFDE(J)                                        LSQ17900
  545 DDUB=PALPHA(K)*PALPHA(J)*DDSE(J,K)                                LSQ18000
      ADUB(J,K)=DUB+DDUB                                                LSQ18100
  520 CONTINUE                                                          LSQ18200
      DTOLER(K)=DDUB                                                    LSQ18300
      ADUB(K,JLAMP1)=PALPHA(K)*SFDE(K)                                  LSQ18400
  510 CONTINUE                                                          LSQ18500
      GO TO 590                                                         LSQ18600
  550 DO 560 J=1,JLAM                                                   LSQ18700
      DTOLER(J)=0.D0                                                    LSQ18800
      R(J)=0.D0                                                         LSQ18900
      DO 565 K=J,JLAMP1                                                 LSQ19000
  565 ADUB(J,K)=0.D0                                                    LSQ19100
  560 CONTINUE                                                          LSQ19200
      DO 570 NU=1,N                                                     LSQ19300
      DO 575 K=1,JLAM                                                   LSQ19400
      IF (.NOT.PIVALP(K)) GO TO 575                                     LSQ19500
      DUB=0.D0                                                          LSQ19600
      DO 578 J=1,NLIN                                                   LSQ19700
  578 DUB=DUB+F(NU,J)*DALPL(J,K)                                        LSQ19800
      DDUB=DF(NU,K)*PALPHA(K)                                           LSQ19900
      DTOLER(K)=DTOLER(K)+DDUB*DDUB                                     LSQ20000
      R(K)=DUB+DDUB                                                     LSQ20100
  575 CONTINUE                                                          LSQ20200
      DDUB=YLYFIT(NU)                                                   LSQ20300
      DO 580 K=1,JLAM                                                   LSQ20400
      IF (.NOT.PIVALP(K)) GO TO 580                                     LSQ20500
      DUB=R(K)                                                          LSQ20600
      DO 585 J=1,K                                                      LSQ20700
  585 ADUB(J,K)=ADUB(J,K)+R(J)*DUB                                      LSQ20800
      ADUB(K,JLAMP1)=ADUB(K,JLAMP1)+R(K)*DDUB                           LSQ20900
  580 CONTINUE                                                          LSQ21000
  570 CONTINUE                                                          LSQ21100
  590 DO 595 J=1,JLAM                                                   LSQ21200
      DTOLER(J)=PRECIS*DMAX1(DTOLER(J),ADUB(J,J))                       LSQ21300
  595 R(J)=ADUB(J,JLAMP1)                                               LSQ21400
      ADUB(JLAMP1,JLAMP1)=0.D0                                          LSQ21500
      CALL PIVOT (ADUB,MDIMA,JLAM,.FALSE.,.FALSE.,.TRUE.,               LSQ21600
     1 DBLE(LAMNMX(1,1)),DBLE(LAMNMX(2,1)),PLAM,DTOLER,DELTAP,LDUM,     LSQ21700
     2 QG,PIVLAM,JLAM,QMAX,IERROR)                                      LSQ21800
      IF (IERROR .EQ. 5) GO TO 800                                      LSQ21900
      ALLPIV(2)=LDUM                                                    LSQ22000
C-----------------------------------------------------------------------LSQ22100
C  USE MODIFICATION-A OF G. E. P. BOX TO ESTIMATE BEST Q (FRACTIONAL    LSQ22200
C  STEP SIZE)                                                           LSQ22300
C-----------------------------------------------------------------------LSQ22400
      VAROLD=VAR                                                        LSQ22500
      DDUM=0.                                                           LSQ22600
      DO 600 J=1,JLAM                                                   LSQ22700
  600 DDUM=DDUM+SNGL(DELTAP(J)*R(J))                                    LSQ22800
      IF (DDUM .GT. 0.) GO TO 610                                       LSQ22900
 5610 FORMAT (42H GAUSS VECTOR POINTING IN WRONG DIRECTION.,2X,         LSQ23000
     1 44(2H**))                                                        LSQ23100
      IF (PRLSQ) WRITE (IWRITE,5610)                                    LSQ23200
      NTRY=2                                                            LSQ23300
      Q=AMIN1(QMIN,QMAX)                                                LSQ23400
      GO TO 215                                                         LSQ23500
  610 CALL VARDUM (QG,JLAM,VARG,T,SQRTW,NMAX,TSTART,DELTAT,NT,          LSQ23600
     1 NINT,.FALSE.,LDUM,E,Y,YLYFIT,NLIN,VAR,NTRY,NDIME,INTER,          LSQ23700
     2 GSE,NDIMG,IERROR)                                                LSQ23800
      IF (IERROR .EQ. 5) GO TO 800                                      LSQ23900
      ALLPIV(1)=LDUM                                                    LSQ24000
      DUM=SNGL(VARG-VAROLD)+2.*DDUM*QG                                  LSQ24100
      IF (DUM .LE. 0.) DUM=.25*DDUM*QG*QG                               LSQ24200
      Q=AMIN1(DDUM*QG*QG/DUM,4.,QMAX)                                   LSQ24300
      IF (VARG.LE.VAROLD .OR. (Q.GE..125*QG .AND. RAWDAT)) GO TO 210    LSQ24400
C-----------------------------------------------------------------------LSQ24500
C  PUT QUADRATIC PARABOLA THRU ACTUAL VARIANCE SURFACE WHEN A           LSQ24600
C      STRONG NONLINEARITY IS INDICATED BY A SMALL Q                    LSQ24700
C-----------------------------------------------------------------------LSQ24800
      QT(1)=AMIN1(.125*QG,AMAX1(Q,.02))                                 LSQ24900
      IF (Q .GT. .125*QG) QT(1)=Q                                       LSQ25000
      CALL VARDUM (QT(1),JLAM,DUB,T,SQRTW,NMAX,TSTART,DELTAT,NT,        LSQ25100
     1 NINT,.FALSE.,LDUM,E,Y,YLYFIT,NLIN,VAR,NTRY,NDIME,INTER,          LSQ25200
     2 GSE,NDIMG,IERROR)                                                LSQ25300
      IF (IERROR .EQ. 5) GO TO 800                                      LSQ25400
      VT(1)=DUB                                                         LSQ25500
      IF (VT(1) .GE. VAROLD) GO TO 615                                  LSQ25600
      QT(2)=4.*QT(1)                                                    LSQ25700
      IF (QT(1) .LE. .125*QG) GO TO 617                                 LSQ25800
      QT(2)=0.                                                          LSQ25900
      VT(2)=VAROLD                                                      LSQ26000
      QT(3)=QG                                                          LSQ26100
      VT(3)=VARG                                                        LSQ26200
      GO TO 630                                                         LSQ26300
  615 QT(2)=.5*QT(1)                                                    LSQ26400
      QT(2)=AMIN1(QT(2),.05)                                            LSQ26500
  617 CALL VARDUM (QT(2),JLAM,DUB,T,SQRTW,NMAX,TSTART,DELTAT,NT,        LSQ26600
     1 NINT,.FALSE.,LDUM,E,Y,YLYFIT,NLIN,VAR,NTRY,NDIME,INTER,          LSQ26700
     2 GSE,NDIMG,IERROR)                                                LSQ26800
      IF (IERROR .EQ. 5) GO TO 800                                      LSQ26900
      VT(2)=DUB                                                         LSQ27000
      IF (VT(1).GE.VAROLD .OR. VT(2).GE.VT(1)) GO TO 620                LSQ27100
      QT(3)=QG                                                          LSQ27200
      VT(3)=VARG                                                        LSQ27300
      GO TO 630                                                         LSQ27400
  620 QT(3)=0.                                                          LSQ27500
      VT(3)=VAROLD                                                      LSQ27600
  630 VARG=DMIN1(VT(1),VT(2),VARG)                                      LSQ27700
      IF (VT(1) .LE. VARG) QG=QT(1)                                     LSQ27800
      IF (VT(2) .LE. VARG) QG=QT(2)                                     LSQ27900
      DUM=SNGL((VT(1)-VT(3)))/(QT(1)-QT(3))                             LSQ28000
      DDUM=(DUM-SNGL((VT(2)-VT(3)))/(QT(2)-QT(3)))/(QT(1)-QT(2))        LSQ28100
      IF (DDUM .GT. 0.) GO TO 650                                       LSQ28200
  640 Q=AMIN1(Q,.125*QT(2),.01)                                         LSQ28300
      GO TO 210                                                         LSQ28400
  650 DUM=.5*((QT(1)+QT(3))-DUM/DDUM)                                   LSQ28500
      IF (DUM .LE. 0.) GO TO 640                                        LSQ28600
      Q=AMIN1(DUM,8.,QMAX)                                              LSQ28700
      GO TO 210                                                         LSQ28800
C-----------------------------------------------------------------------LSQ28900
C  CHECK FOR ABNORMAL EXIT AND CALL ANLERR                              LSQ29000
C-----------------------------------------------------------------------LSQ29100
  800 FAILED=IERROR.NE.0 .OR. .NOT.ALLPIV(1)                            LSQ29200
      IF (.NOT.RAWDAT) GO TO 890                                        LSQ29300
      FAILED=FAILED  .OR. .NOT.ALLPIV(2)                                LSQ29400
      IF (.NOT.FAILED .OR. ITYPE.EQ.4) CALL ANLERR (JLAM,T,SQRTW,NMAX,E,LSQ29500
     1 NLIN,NDIME,TSTART,DELTAT,NT,NINT,N,ITYPE,INTER,PRLSQ)            LSQ29600
  890 IF (.NOT.(PRLSQ .AND. FAILED)) GO TO 950                          LSQ29700
 5889 FORMAT (/47H VARIANCE DID NOT SIGNIFICANTLY DECREASE IN THE,      LSQ29800
     1 5H LAST,I2,7H TIMES.,5X,6(11H***TOO SLOW))                       LSQ29900
      IF (IERROR .EQ. 3) WRITE (IWRITE,5889) NABORT                     LSQ30000
 5890 FORMAT (/41H MAXIMUM ITERATIONS WITHOUT CONVERGENCE. ,            LSQ30100
     1 7(13H***ITERATIONS))                                             LSQ30200
      IF (IERROR .EQ. 2) WRITE (IWRITE,5890)                            LSQ30300
 5892 FORMAT (/28H VARIANCE INCREASED THE LAST,I2,8H TIMES. ,           LSQ30400
     1 3X,7(13H***DIVERGENCE))                                          LSQ30500
      IF (IERROR .EQ. 1) WRITE (IWRITE,5892) NVARUP                     LSQ30600
 5894 FORMAT (/42H NEARLY SINGULAR NORMAL EQUATIONS MATRIX. ,           LSQ30700
     1 6(15H****SINGULARITY))                                           LSQ30800
      IF (IERROR .EQ. 5) WRITE (IWRITE,5894)                            LSQ30900
C-----------------------------------------------------------------------LSQ31000
C  PUT PLAM IN ASCENDING ORDER                                          LSQ31100
C-----------------------------------------------------------------------LSQ31200
  950 DO 960 J=1,JLAM                                                   LSQ31300
      DUB=1.E+20                                                        LSQ31400
      DO 985 K=J,JLAM                                                   LSQ31500
      IF (PLAM(K) .GT. DUB) GO TO 985                                   LSQ31600
      L=K                                                               LSQ31700
      DUB=PLAM(K)                                                       LSQ31800
  985 CONTINUE                                                          LSQ31900
      PLAM(L)=PLAM(J)                                                   LSQ32000
      PLAM(J)=DUB                                                       LSQ32100
      DUB=PALPHA(L)                                                     LSQ32200
      PALPHA(L)=PALPHA(J)                                               LSQ32300
      PALPHA(J)=DUB                                                     LSQ32400
      IF (.NOT.RAWDAT) GO TO 960                                        LSQ32500
      DUM=SIGMAP(L)                                                     LSQ32600
      SIGMAP(L)=SIGMAP(J)                                               LSQ32700
      SIGMAP(J)=DUM                                                     LSQ32800
      K=L+JLAM                                                          LSQ32900
      JK=J+JLAM                                                         LSQ33000
      DUM=SIGMAP(K)                                                     LSQ33100
      SIGMAP(K)=SIGMAP(JK)                                              LSQ33200
      SIGMAP(JK)=DUM                                                    LSQ33300
      PCERR(J)=100.*SIGMAP(J)/SNGL(PLAM(J))                             LSQ33400
      IF (DABS(PALPHA(J)) .GT. 0.D0) PCERR(JK)=100.*SIGMAP(JK)/ABS(     LSQ33500
     1 SNGL(PALPHA(J)))                                                 LSQ33600
  960 CONTINUE                                                          LSQ33700
      RETURN                                                            LSQ33800
      END                                                               LSQ33900
C--------------------- DOUBLE PRECISION VERSION 2B ---------------------EVR00100
C  SUBROUTINE EVAR.  FOR FRACTIONAL STEP SIZE Q, EVALUATES PALPHA       EVR00200
C      (LINEAR LEAST SQUARES PARAMETERS), VAR (WEIGHTED VARIANCE) AND   EVR00300
C      ACCUMULATES IT IN VGRID2 (COARSE GRID SUM), AND, IF              EVR00400
C      INVERT=.TRUE., EVALUATES ARRAYS SFDE, DFCAPL, AND FULL UPPER     EVR00500
C      TRIANGLES OF DDSE AND SE FOR LATER USE IN NONLINEAR LEAST        EVR00600
C      SQUARES.                                                         EVR00700
C-----------------------------------------------------------------------EVR00800
C  CALLS SUBPROGRAMS - ETHEOR, PIVOT                                    EVR00900
C  WHICH IN TURN CALL - PIVOT1                                          EVR01000
C-----------------------------------------------------------------------EVR01100
      SUBROUTINE EVAR (Q,JLAM,VAR,T,SQRTW,NMAX,TSTART,DELTAT,NT,NINT,   EVR01200
     1 INVERT,ALLPIV,E,Y,YLYDUM,NLIN,VARG,NTRY,NDIME,INTER,GSE,NDIMG,   EVR01300
     2 IERROR)                                                          EVR01400
      LOGICAL INVERT, ALLPIV, LDDUM, INTER, LBLOKA, LBLOKB, LBLOKC,     EVR01500
     1 NOTINV, NOBASE, WTED, REGINT, NONNEG, PIVALP                     EVR01600
      REAL LAMNMX                                                       EVR01700
      DOUBLE PRECISION TSTART, DELTAT, E, YLYDUM, GSE, SUMWTY, SUMWT,   EVR01800
     1 ZERO, DGRIDE, GESTL, PALPHA, PLAM, SE, PLMTRY, DELTAP, DFCAPL,   EVR01900
     2 SFDE, ADUB, DDSE, DTOLER, DBLOKA, DBLOKB, DBLOKC, ADUM, DUM,     EVR02000
     3 DDUM, DDDUM, XLAM, HZ, H, HH, DX, DDX, S, FACTOR, XEXP,          EVR02100
     4 YLYFIT, VAR, VARG, T                                             EVR02200
      COMMON /BLOK/ DBLOKA(4829), SUMWTY, SUMWT, ZERO(10), DGRIDE,      EVR02300
     1 GESTL, DBLOKB, PALPHA(10), PLAM(10), SE(10,10), PLMTRY(10),      EVR02400
     2 DELTAP(19), DFCAPL(10,9), SFDE(9), ADUB(19,19), DDSE(9,9),       EVR02500
     3 DTOLER(19), DBLOKC(380), RBLOKA(174), PRECIS, LAMNMX(2,2),       EVR02600
     4 RBLOKB(40), VGRID2(252), JGEADD, NBINOM(20,9), IBLOKA(9),        EVR02700
     5 NPERG2(9), MDIMA, IBLOKB(4), NGRID2(9), IBLOKC(84), N            EVR02800
      COMMON /BLOK/ IBLOKD(17), REGINT, LBLOKA(2), NOBASE, WTED,        EVR02900
     1 LBLOKB(3), NONNEG, LBLOKC(18), PIVALP(19)                        EVR03000
      DIMENSION ADUM(10), FACTOR(9), XEXP(9), JG2(9)                    EVR03100
      DIMENSION T(NMAX), SQRTW(NMAX), TSTART(NINT), DELTAT(NINT),       EVR03200
     1 NT(NINT), E(NDIME,NLIN), Y(NMAX), YLYDUM(1), GSE(NDIMG,4)        EVR03300
      DO 110 J=1,JLAM                                                   EVR03400
  110 PLMTRY(J)=DMAX1(DBLE(LAMNMX(1,1)),                                EVR03500
     1                DMIN1(DBLE(LAMNMX(2,1)),PLAM(J)+Q*DELTAP(J)))     EVR03501
      NLINP=NLIN+1                                                      EVR03600
      DO 160 K=1,NLINP                                                  EVR03700
  160 ADUB(K,NLINP)=0.D0                                                EVR03800
      IF (NOBASE) GO TO 170                                             EVR03900
      PLMTRY(NLIN)=0.D0                                                 EVR04000
      ADUB(NLIN,NLIN)=SUMWT                                             EVR04100
      ADUB(NLIN,NLINP)=SUMWTY                                           EVR04200
  170 NOTINV=.NOT.INVERT                                                EVR04300
      IF (REGINT) GO TO 300                                             EVR04400
C-----------------------------------------------------------------------EVR04500
C  ORDINARY COMPUTATION OF E (EXPONENTIALS), ADUB (NORMAL EQUATIONS     EVR04600
C      MATRIX FOR LINEAR LEAST SQUARES) AND, IF INVERT=.TRUE., PART     EVR04700
C      OF DFCAPL (PARTIAL OF FCAP W.R.T. LAMBDA) WHEN REGINT=.FALSE.    EVR04800
C-----------------------------------------------------------------------EVR04900
      DO 190 K=1,NLIN                                                   EVR05000
      DO 195 J=1,K                                                      EVR05100
      IF (J .GT. JLAM) GO TO 190                                        EVR05200
      ADUB(J,K)=0.D0                                                    EVR05300
      DFCAPL(K,J)=0.D0                                                  EVR05400
  195 IF (K .LE. JLAM) DDSE(J,K)=0.D0                                   EVR05500
  190 CONTINUE                                                          EVR05600
      DO 200 J=1,N                                                      EVR05700
      DUM=-T(J)                                                         EVR05800
      DO 210 K=1,JLAM                                                   EVR05900
      ADUM(K)=DBLE(SQRTW(J))*DEXP(PLMTRY(K)*DUM)                        EVR06000
  210 E(J,K)=ADUM(K)                                                    EVR06100
      IF (NOBASE) GO TO 215                                             EVR06200
      ADUM(NLIN)=DBLE(SQRTW(J))                                         EVR06300
      E(J,NLIN)=ADUM(NLIN)                                              EVR06400
  215 DO 220 K=1,NLIN                                                   EVR06500
      DDUM=ADUM(K)                                                      EVR06600
      LDDUM=K .LE. JLAM                                                 EVR06700
      DO 230 L=1,K                                                      EVR06800
      IF (L .GT. JLAM) GO TO 220                                        EVR06900
      DDDUM=ADUM(L)*DDUM                                                EVR07000
      ADUB(L,K)=ADUB(L,K)+DDDUM                                         EVR07100
      IF (NOTINV) GO TO 230                                             EVR07200
      DDDUM=DDDUM*DUM                                                   EVR07300
      DFCAPL(K,L)=DFCAPL(K,L)-DDDUM                                     EVR07400
      IF (LDDUM) DDSE(L,K)=DDSE(L,K)+DDDUM*DUM                          EVR07500
  230 CONTINUE                                                          EVR07600
      ADUB(K,NLINP)=ADUB(K,NLINP)+DDUM*DBLE(Y(J))                       EVR07700
  220 CONTINUE                                                          EVR07800
  200 CONTINUE                                                          EVR07900
      GO TO 500                                                         EVR08000
C-----------------------------------------------------------------------EVR08100
C  EVALUATE ADUB AND, IF INVERT=.TRUE., PART OF DFCAPL AND DDSE         EVR08200
C      (SECOND DERIVATIVES OF EXPONENTIAL SUMS W.R.T. LAMBDA)           EVR08300
C      BY HERMITE AND LAGRANGE INTERPOLATION, IF REGINT=INTER=.TRUE.    EVR08400
C-----------------------------------------------------------------------EVR08500
  300 IF (.NOT.INTER) GO TO 400                                         EVR08600
      DO 310 K=1,NLIN                                                   EVR08700
      LDDUM=K .LE. JLAM                                                 EVR08800
      DO 320 J=1,K                                                      EVR08900
      IF (J .GT. JLAM) GO TO 310                                        EVR09000
      XLAM=PLMTRY(K)+PLMTRY(J)                                          EVR09100
      DUM=(DLOG(XLAM)-GESTL)/DGRIDE                                     EVR09200
      JZ=IDINT(DUM+.5D0)                                                EVR09300
      HZ=DUM-DBLE(FLOAT(JZ))                                            EVR09400
      IF (DABS(HZ) .GT. 1.D-5) GO TO 330                                EVR09500
      ADUB(J,K)=GSE(JZ,1)                                               EVR09600
      IF (NOTINV) GO TO 320                                             EVR09700
      DFCAPL(K,J)=-GSE(JZ,2)/XLAM                                       EVR09800
      IF (LDDUM) DDSE(J,K)=(GSE(JZ,3)-GSE(JZ,2))/(XLAM*XLAM)            EVR09900
      GO TO 320                                                         EVR10000
  330 JM=IDINT(DUM)                                                     EVR10100
      JP=JM+1                                                           EVR10200
      H=DUM-DBLE(FLOAT(JM))                                             EVR10300
      HH=1.D0-H                                                         EVR10400
      DDUM=1.D0+3.D0*H                                                  EVR10500
      DDDUM=1.D0+3.D0*HH                                                EVR10600
      DX=H*DGRIDE                                                       EVR10700
      DDX=-HH*DGRIDE                                                    EVR10800
      ADUB(J,K)=HH*HH*HH*((DDUM+6.D0*H*H)*GSE(JM,1)+DX*(DDUM*GSE(JM,2)+ EVR10900
     1 .5D0*DX*GSE(JM,3)))+H*H*H*((DDDUM+6.D0*HH*HH)*GSE(JP,1)+         EVR11000
     2 DDX*(DDDUM*GSE(JP,2)+.5D0*DDX*GSE(JP,3)))                        EVR11100
      IF (NOTINV) GO TO 320                                             EVR11200
      HH=HZ*HZ                                                          EVR11300
      DX=(HH-1.D0)*HZ*.5D0                                              EVR11400
      DUM=1.D0/(HZ+1.D0)                                                EVR11500
      DDUM=1.D0/HZ                                                      EVR11600
      DDDUM=1.D0/(HZ-1.D0)                                              EVR11700
      H=DX*DX*(((3.D0*HZ+4.D0)*GSE(JZ-1,2)*DUM+DGRIDE*GSE(JZ-1,3))*DUM+ EVR11800
     1 4.D0*DDUM*(DDUM*GSE(JZ,2)+DGRIDE*GSE(JZ,3))+DDDUM*((4.D0-3.D0*   EVR11900
     2 HZ)*GSE(JZ+1,2)*DDDUM+DGRIDE*GSE(JZ+1,3)))                       EVR12000
      DFCAPL(K,J)=-H/XLAM                                               EVR12100
      IF (.NOT.LDDUM) GO TO 320                                         EVR12200
      DUM=DX*(HH-4.D0)*(HH-9.D0)*(GSE(JZ+3,3)/(HZ-3.D0)+GSE(JZ-3,3)/    EVR12300
     1 (HZ+3.D0)-6.D0*(GSE(JZ+2,3)/(HZ-2.D0)+GSE(JZ-2,3)/(HZ+2.D0))+    EVR12400
     2 15.D0*(GSE(JZ+1,3)*DDDUM+GSE(JZ-1,3)*DUM)-20.D0*GSE(JZ,3)*DDUM)/ EVR12500
     3 360.D0                                                           EVR12600
      DDSE(J,K)=(DUM-H)/(XLAM*XLAM)                                     EVR12700
  320 CONTINUE                                                          EVR12800
      DUM=(DLOG(PLMTRY(K))-GESTL)/DGRIDE                                EVR12900
      JZ=IDINT(DUM+.5D0)                                                EVR13000
      HZ=DUM-DBLE(FLOAT(JZ))                                            EVR13100
      IF (DABS(HZ) .GT. 1.D-5) GO TO 340                                EVR13200
      ADUB(K,NLINP)=GSE(JZ,4)                                           EVR13300
      GO TO 310                                                         EVR13400
  340 HH=HZ*HZ                                                          EVR13500
      ADUB(K,NLINP)=HZ*(HH-1.D0)*(HH-4.D0)*(HH-9.D0)*(GSE(JZ+3,4)/(HZ-  EVR13600
     1 3.D0)+GSE(JZ-3,4)/(HZ+3.D0)-6.D0*(GSE(JZ+2,4)/(HZ-2.D0)+         EVR13700
     2 GSE(JZ-2,4)/(HZ+2.D0))+15.D0*(GSE(JZ+1,4)/(HZ-1.D0)+GSE(JZ-1,4)/ EVR13800
     3 (HZ+1.D0))-20.D0*GSE(JZ,4)/HZ)/720.D0                            EVR13900
  310 CONTINUE                                                          EVR14000
      GO TO 500                                                         EVR14100
C-----------------------------------------------------------------------EVR14200
C  EVALUATE ADUB, PART OF DFCAPL, AND DDSE BY CALLING ETHEOR            EVR14300
C      IF REGINT=.TRUE. AND INTER=.FALSE.                               EVR14400
C-----------------------------------------------------------------------EVR14500
  400 DO 410 K=1,NLIN                                                   EVR14600
      LDDUM=K .LE. JLAM                                                 EVR14700
      DO 420 J=1,K                                                      EVR14800
      IF (J .GT. JLAM) GO TO 410                                        EVR14900
      XLAM=PLMTRY(K)+PLMTRY(J)                                          EVR15000
      CALL ETHEOR (XLAM,INVERT,WTED,T,SQRTW,Y,NMAX,TSTART,DELTAT,       EVR15100
     1 NT,NINT,DUM,DDUM,DDDUM,S,.TRUE.,.FALSE.)                         EVR15200
      ADUB(J,K)=DUM                                                     EVR15300
      IF (NOTINV) GO TO 420                                             EVR15400
      DFCAPL(K,J)=-DDUM/XLAM                                            EVR15500
      IF (LDDUM) DDSE(J,K)=(DDDUM-DDUM)/(XLAM*XLAM)                     EVR15600
  420 CONTINUE                                                          EVR15700
      CALL ETHEOR (PLMTRY(K),.FALSE.,WTED,T,SQRTW,Y,NMAX,TSTART,DELTAT, EVR15800
     1 NT,NINT,DUM,DDUM,DDDUM,S,.FALSE.,.TRUE.)                         EVR15900
      ADUB(K,NLINP)=S                                                   EVR16000
  410 CONTINUE                                                          EVR16100
  500 DO 510 K=1,NLIN                                                   EVR16200
      DTOLER(K)=10.*PRECIS*ADUB(K,K)                                    EVR16300
      DO 515 J=1,K                                                      EVR16400
  515 SE(J,K)=ADUB(J,K)                                                 EVR16500
  510 CONTINUE                                                          EVR16600
C-----------------------------------------------------------------------EVR16700
C  SOLVE NORMAL EQUATIONS FOR PALPHA USING STEPWISE REGRESSION          EVR16800
C-----------------------------------------------------------------------EVR16900
      CALL PIVOT (ADUB,MDIMA,NLIN,.FALSE.,INVERT,NONNEG,0.D0,           EVR17000
     1 1.D+20,ZERO,DTOLER,PALPHA,ALLPIV,SP,PIVALP,JLAM,SP,IERROR)       EVR17100
      IF (IERROR .NE. 5) GO TO 518                                      EVR17200
      VAR=1.D+20                                                        EVR17300
      RETURN                                                            EVR17400
C-----------------------------------------------------------------------EVR17500
C  CALCULATE VAR AND, IF INVERT=.TRUE., SFDE (SUMS OF RESIDUALS         EVR17600
C      TIMES DERIVATIVES OF NONLINEAR TERMS)                            EVR17700
C-----------------------------------------------------------------------EVR17800
  518 VAR=0.D0                                                          EVR17900
      IF (NOTINV) GO TO 528                                             EVR18000
      DO 520 K=1,JLAM                                                   EVR18100
      SFDE(K)=0.D0                                                      EVR18200
      DO 525 J=1,K                                                      EVR18300
  525 DFCAPL(J,K)=DFCAPL(K,J)                                           EVR18400
  520 CONTINUE                                                          EVR18500
  528 IF (.NOT.REGINT) GO TO 600                                        EVR18600
      K=0                                                               EVR18700
      DO 530 J=1,NINT                                                   EVR18800
      DO 540 L=1,JLAM                                                   EVR18900
      FACTOR(L)=DEXP(-PLMTRY(L)*DELTAT(J))                              EVR19000
  540 XEXP(L)=DEXP(-PLMTRY(L)*TSTART(J))                                EVR19100
      KK=NT(J)                                                          EVR19200
      DO 550 LL=1,KK                                                    EVR19300
      K=K+1                                                             EVR19400
      DUM=DBLE(SQRTW(K))                                                EVR19500
      YLYFIT=DBLE(Y(K))                                                 EVR19600
      IF (.NOT.NOBASE) YLYFIT=YLYFIT-PALPHA(NLIN)*DUM                   EVR19700
      DO 560 L=1,JLAM                                                   EVR19800
      ADUM(L)=XEXP(L)*DUM                                               EVR19900
      YLYFIT=YLYFIT-ADUM(L)*PALPHA(L)                                   EVR20000
      XEXP(L)=XEXP(L)*FACTOR(L)                                         EVR20100
  560 CONTINUE                                                          EVR20200
      VAR=VAR+YLYFIT*YLYFIT                                             EVR20300
      IF (NOTINV) GO TO 550                                             EVR20400
      DDUM=T(K)                                                         EVR20500
      DO 570 L=1,JLAM                                                   EVR20600
  570 SFDE(L)=SFDE(L)-DDUM*YLYFIT*ADUM(L)                               EVR20700
  550 CONTINUE                                                          EVR20800
  530 CONTINUE                                                          EVR20900
      GO TO 700                                                         EVR21000
  600 DO 610 J=1,N                                                      EVR21100
      YLYFIT=DBLE(Y(J))                                                 EVR21200
      DO 620 K=1,NLIN                                                   EVR21300
  620 YLYFIT=YLYFIT-PALPHA(K)*E(J,K)                                    EVR21400
      VAR=VAR+YLYFIT*YLYFIT                                             EVR21500
      IF (NOTINV) GO TO 610                                             EVR21600
      DUM=T(J)                                                          EVR21700
      DO 630 K=1,JLAM                                                   EVR21800
  630 SFDE(K)=SFDE(K)-DUM*YLYFIT*E(J,K)                                 EVR21900
  610 CONTINUE                                                          EVR22000
  700 IF (JLAM.EQ.1 .OR. .NOT.INTER) GO TO 800                          EVR22100
C-----------------------------------------------------------------------EVR22200
C  ADD VAR INTO VGRID2 (COARSE GRID SUM)                                EVR22300
C-----------------------------------------------------------------------EVR22400
      J2MAX=0                                                           EVR22500
      DO 710 L=1,JLAM                                                   EVR22600
      IF (PLMTRY(L).LT.DBLE(LAMNMX(1,2)) .OR. PLMTRY(L).GT.             EVR22700
     1 DBLE(LAMNMX(2,2))) GO TO 800                                     EVR22800
      JG2(L)=(IDINT((DLOG(PLMTRY(L))-GESTL)/DGRIDE+.5D0)-JGEADD)/       EVR22900
     1 NPERG2(JLAM)+1                                                   EVR23000
      IF (JG2(L) .LT. J2MAX) GO TO 710                                  EVR23100
      J2MAX=JG2(L)                                                      EVR23200
      LMAX=L                                                            EVR23300
  710 CONTINUE                                                          EVR23400
      JG2(LMAX)=0                                                       EVR23500
      LOC=1                                                             EVR23600
      DO 720 J=1,JLAM                                                   EVR23700
      J2NEXT=0                                                          EVR23800
      DO 730 L=1,JLAM                                                   EVR23900
      IF (JG2(L) .LT. J2NEXT) GO TO 730                                 EVR24000
      J2NEXT=JG2(L)                                                     EVR24100
      LNEXT=L                                                           EVR24200
  730 CONTINUE                                                          EVR24300
      L=J2NEXT+1                                                        EVR24400
      IF (J2MAX-L) 800,750,740                                          EVR24500
  740 LL=J2MAX-1                                                        EVR24600
      DO 745 K=L,LL                                                     EVR24700
      KK=NGRID2(JLAM)-K+1                                               EVR24800
  745 LOC=LOC+NBINOM(KK,J)                                              EVR24900
  750 J2MAX=J2NEXT                                                      EVR25000
      JG2(LNEXT)=0                                                      EVR25100
  720 CONTINUE                                                          EVR25200
      SP=SNGL(VAR)                                                      EVR25300
      LOC=MIN0(LOC,252)                                                 EVR25400
      IF (.NOT.INVERT .OR. (VAR.GT.VARG .AND. NTRY.EQ.1)) SP=SP*1.E-6   EVR25500
      VGRID2(LOC)=VGRID2(LOC)+SP                                        EVR25600
  800 RETURN                                                            EVR25700
      END                                                               EVR25800
C--------------------- DOUBLE PRECISION VERSION 2B ---------------------VRF00100
C  SUBROUTINE VARF.  SAME AS EVAR, EXCEPT FOR TRANSFORMS INSTEAD OF     VRF00200
C      EXPONENTIALS AND DOES NO ACCUMULATION IN A COARSE GRID.          VRF00300
C-----------------------------------------------------------------------VRF00400
C  CALLS SUBPROGRAMS - PIVOT                                            VRF00500
C  WHICH IN TURN CALL - PIVOT1                                          VRF00600
C-----------------------------------------------------------------------VRF00700
      SUBROUTINE VARF (Q,JLAM,VAR,T,SQRTW,NMAX,TSTART,DELTAT,NT,NINT,   VRF00800
     1 INVERT,ALLPIV,E,Y,YLYFIT,NLIN,VARG,NTRY,NDIME,INTDUM,GSE,NDIMG,  VRF00900
     2 IERROR)                                                          VRF01000
      LOGICAL INVERT, ALLPIV, INTDUM, LBLOKA, LBLOKB, LBLOKC,           VRF01100
     1 NOBASE, NONNEG, PIVALP, NOTINV                                   VRF01200
      DOUBLE PRECISION TRBASE, GF, DGFL, DGRIDF, GFSTL, FHAT, ZERO,     VRF01300
     1 PALPHA, PLAM, PLMTRY, DELTAP, DFCAPL, ADUB, DTOLER, F, DF,       VRF01400
     2 DBLOKA, DBLOKB, DBLOKC, DBLOKD, DBLOKE, VAR, TSTART, DELTAT,     VRF01500
     3 E, YLYFIT, VARG, GSE, RLAM, DUM, ADUM, HH, HZ, DDUM, DDDUM, T    VRF01600
      COMMON /BLOK/ TRBASE(19), GF(19,126), DGFL(19,126), DGRIDF,       VRF01700
     1 GFSTL, FHAT(19), DBLOKA(3), ZERO(10), DBLOKB(3), PALPHA(10),     VRF01800
     2 PLAM(10), DBLOKC(100), PLMTRY(10), DELTAP(19), DFCAPL(10,9),     VRF01900
     3 DBLOKD(9), ADUB(19,19), DBLOKE(9,9), DTOLER(19), F(19,10),       VRF02000
     4 DF(19,10), RBLOKA(174), PRECIS, RBLOKB(296), IBLOKA(199),        VRF02100
     5 MDIMA, IBLOKB(96), NF, IBLOKC(18), LBLOKA(3), NOBASE             VRF02200
      COMMON /BLOK/ LBLOKB(4), NONNEG, LBLOKC(18), PIVALP(19)           VRF02300
      DIMENSION T(1), SQRTW(1), TSTART(1), DELTAT(1), NT(1), E(1,1),    VRF02400
     1 Y(1), YLYFIT(NMAX), ADUM(9), GSE(1,1)                            VRF02500
      NLINP=NLIN+1                                                      VRF02600
      NOTINV=.NOT.INVERT                                                VRF02700
      IF (NOBASE) GO TO 190                                             VRF02800
      PLMTRY(NLIN)=0.D0                                                 VRF02900
      DO 130 J=1,NF                                                     VRF03000
  130 F(J,NLIN)=TRBASE(J)                                               VRF03100
C-----------------------------------------------------------------------VRF03200
C  COMPUTE F (NONLINEAR TERMS) AND, IF INVERT=.TRUE., DF (THEIR         VRF03300
C      DERIVATIVES W.R.T. LAMBDA) BY HERMITE AND LAGRANGE               VRF03400
C      INTERPOLATION.                                                   VRF03500
C-----------------------------------------------------------------------VRF03600
  190 DO 200 J=1,JLAM                                                   VRF03700
      PLMTRY(J)=DMAX1(DBLE(RBLOKB(1)),                                  VRF03800
     1                DMIN1(DBLE(RBLOKB(2)),PLAM(J)+Q*DELTAP(J)))       VRF03801
      RLAM=1.D0/PLMTRY(J)                                               VRF03900
      DUM=(DLOG(PLMTRY(J))-GFSTL)/DGRIDF                                VRF04000
      JZ=IDINT(DUM+.5D0)                                                VRF04100
      HZ=DUM-DBLE(FLOAT(JZ))                                            VRF04200
      IF (DABS(HZ) .GT. 1.D-5) GO TO 220                                VRF04300
      DO 210 K=1,NF                                                     VRF04400
      F(K,J)=GF(K,JZ)                                                   VRF04500
  210 IF (INVERT) DF(K,J)=DGFL(K,JZ)*RLAM                               VRF04600
      GO TO 200                                                         VRF04700
  220 ADUM(7)=1.D0/HZ                                                   VRF04800
      ADUM(1)=4.D0*ADUM(7)*ADUM(7)                                      VRF04900
      ADUM(2)=4.D0*DGRIDF*ADUM(7)                                       VRF05000
      ADUM(8)=1.D0/(HZ+1.D0)                                            VRF05100
      ADUM(3)=(3.D0*HZ+4.D0)*ADUM(8)*ADUM(8)                            VRF05200
      ADUM(4)=DGRIDF*ADUM(8)                                            VRF05300
      ADUM(9)=1.D0/(HZ-1.D0)                                            VRF05400
      ADUM(5)=(4.D0-3.D0*HZ)*ADUM(9)*ADUM(9)                            VRF05500
      ADUM(6)=DGRIDF*ADUM(9)                                            VRF05600
      HH=HZ*HZ                                                          VRF05700
      DDUM=(HH-1.D0)*HZ*.5D0                                            VRF05800
      DDDUM=DDUM*DDUM                                                   VRF05900
      DO 230 K=1,NF                                                     VRF06000
  230 F(K,J)=DDDUM*(ADUM(1)*GF(K,JZ)+ADUM(2)*DGFL(K,JZ)+ADUM(3)*        VRF06100
     1 GF(K,JZ-1)+ADUM(4)*DGFL(K,JZ-1)+ADUM(5)*GF(K,JZ+1)+ADUM(6)*      VRF06200
     2 DGFL(K,JZ+1))                                                    VRF06300
      IF (NOTINV) GO TO 200                                             VRF06400
      ADUM(1)=-20.D0*ADUM(7)                                            VRF06500
      ADUM(2)=15.D0*ADUM(8)                                             VRF06600
      ADUM(3)=15.D0*ADUM(9)                                             VRF06700
      ADUM(4)=-6.D0/(HZ+2.D0)                                           VRF06800
      ADUM(5)=-6.D0/(HZ-2.D0)                                           VRF06900
      ADUM(6)=1.D0/(HZ+3.D0)                                            VRF07000
      ADUM(7)=1.D0/(HZ-3.D0)                                            VRF07100
      DDDUM=DDUM*(HH-4.D0)*(HH-9.D0)*RLAM/360.D0                        VRF07200
      DO 240 K=1,NF                                                     VRF07300
  240 DF(K,J)=DDDUM*(ADUM(1)*DGFL(K,JZ)+ADUM(2)*DGFL(K,JZ-1)+ADUM(3)*   VRF07400
     1 DGFL(K,JZ+1)+ADUM(4)*DGFL(K,JZ-2)+ADUM(5)*DGFL(K,JZ+2)+          VRF07500
     2 ADUM(6)*DGFL(K,JZ-3)+ADUM(7)*DGFL(K,JZ+3))                       VRF07600
  200 CONTINUE                                                          VRF07700
C-----------------------------------------------------------------------VRF07800
C  COMPUTE AND INVERT ADUB (NORMAL EQUATIONS MATRIX FOR LINEAR          VRF07900
C      LEAST SQUARES                                                    VRF08000
C-----------------------------------------------------------------------VRF08100
      ADUB(NLINP,NLINP)=0.D0                                            VRF08200
      DO 300 K=1,NLIN                                                   VRF08300
      DO 310 J=1,K                                                      VRF08400
      DUM=0.D0                                                          VRF08500
      DO 320 L=1,NF                                                     VRF08600
  320 DUM=DUM+F(L,J)*F(L,K)                                             VRF08700
      ADUB(J,K)=DUM                                                     VRF08800
  310 CONTINUE                                                          VRF08900
      DTOLER(K)=10.*PRECIS*ADUB(K,K)                                    VRF09000
      DUM=0.D0                                                          VRF09100
      DO 330 L=1,NF                                                     VRF09200
  330 DUM=DUM+F(L,K)*FHAT(L)                                            VRF09300
      ADUB(K,NLINP)=DUM                                                 VRF09400
  300 CONTINUE                                                          VRF09500
      CALL PIVOT (ADUB,MDIMA,NLIN,.FALSE.,INVERT,NONNEG,0.D0,           VRF09600
     1 1.D+20,ZERO,DTOLER,PALPHA,ALLPIV,SP,PIVALP,JLAM,SP,IERROR)       VRF09700
          IF (IERROR .EQ. 5) RETURN                                     VRF09800
C-----------------------------------------------------------------------VRF09900
C  CALCULATE VAR (WEIGHTED VARIANCE OF FIT) AND YLYFIT (WEIGHTED        VRF10000
C      RESIDUALS), AND,IF INVERT=.TRUE., PART OF DFCAPL FOR             VRF10100
C      LATER USE IN NONLINEAR LEAST SQUARES                             VRF10200
C-----------------------------------------------------------------------VRF10300
      VAR=0.D0                                                          VRF10400
      DO 410 L=1,NF                                                     VRF10500
      DUM=FHAT(L)                                                       VRF10600
      DO 420 J=1,NLIN                                                   VRF10700
  420 DUM=DUM-PALPHA(J)*F(L,J)                                          VRF10800
      YLYFIT(L)=DUM                                                     VRF10900
      VAR=VAR+DUM*DUM                                                   VRF11000
  410 CONTINUE                                                          VRF11100
      IF (NOTINV .OR. (VAR.GT.VARG .AND. NTRY.EQ.1)) RETURN             VRF11200
      DO 450 K=1,JLAM                                                   VRF11300
      DO 470 J=1,NLIN                                                   VRF11400
      DUM=0.D0                                                          VRF11500
      DO 480 L=1,NF                                                     VRF11600
  480 DUM=DUM-F(L,J)*DF(L,K)                                            VRF11700
      DFCAPL(J,K)=DUM                                                   VRF11800
  470 CONTINUE                                                          VRF11900
  450 CONTINUE                                                          VRF12000
      RETURN                                                            VRF12100
      END                                                               VRF12200
C--------------------- DOUBLE PRECISION VERSION 2B ---------------------ETH00100
C  SUBROUTINE ETHEOR.  IF DOE=.TRUE., PUTS EXPONENTIAL SUM AND, IF      ETH00200
C      INVERT=.TRUE., ITS FIRST AND SECOND DERIVATIVES W.R.T.           ETH00300
C      ALOG(LAMBDA) IN SE, DSE, DDSE.  IF DOYE=.TRUE., PUTS SUM OF Y    ETH00400
C      TIMES EXPONENTIALS IN SYE.  DONE IN DOUBLE PRECISION.            ETH00500
C  CAN ONLY BE USED WHEN REGINT=.TRUE..                                 ETH00600
C  IF WTED=.TRUE., USES RECURSION RELATIONS, AND IF WTED=.FALSE. USES   ETH00700
C      GEOMETRIC SUM FORMULAS FOR SE, DSE, AND DDSE.                    ETH00800
C  NO COMMON BLOCKS NEEDED.                                             ETH00900
C-----------------------------------------------------------------------ETH01000
C  CALLS NO SUBPROGRAMS.                                                ETH01100
C-----------------------------------------------------------------------ETH01200
      SUBROUTINE ETHEOR (PLAMBD,INVERT,WTED,T,SQRTW,Y,NMAX,TSTART,      ETH01300
     1 DELTAT,NT,NINT,SE,DSE,DDSE,SYE,DOE,DOYE)                         ETH01400
      DOUBLE PRECISION EXPSMA, AN, B, C, DDDUM, DDSEI, DDUM,            ETH01500
     1 DEDNL, DL, DNL, DSEI, DUM, EDL, ELTS, ENTS, RDEDL, SEI, S,       ETH01600
     2 SS, SYEI, PLAMBD, SE, DSE, DDSE, SYE, TSTART, DELTAT, T          ETH01700
      LOGICAL INVERT, NOTINV, WTED, DOE, DOYE, LDDUM                    ETH01800
      DIMENSION SQRTW(NMAX), TSTART(NINT), DELTAT(NINT), NT(NINT),      ETH01900
     1 T(NMAX), Y(NMAX)                                                 ETH02000
      EXPSMA(S)=S*(720.D0-S*(360.D0-S*(120.D0-S*(30.D0-S*(6.D0-S)))))/  ETH02100
     1 720.D0                                                           ETH02200
      IF (DOE) SEI=0.D0                                                 ETH02300
      IF (INVERT) DSEI=0.D0                                             ETH02400
      IF (INVERT) DDSEI=0.D0                                            ETH02500
      NOTINV=.NOT.INVERT                                                ETH02600
      IF (.NOT.(DOYE .OR. WTED)) GO TO 200                              ETH02700
      IF (DOYE) SYEI=0.D0                                               ETH02800
      LDDUM=.NOT.(DOE .AND. WTED)                                       ETH02900
      J=0                                                               ETH03000
      DO 110 L=1,NINT                                                   ETH03100
      DDUM=DEXP(-PLAMBD*DELTAT(L))                                      ETH03200
      DDDUM=DEXP(-PLAMBD*TSTART(L))                                     ETH03300
      KK=NT(L)                                                          ETH03400
      DO 120 K=1,KK                                                     ETH03500
      J=J+1                                                             ETH03600
      DUM=DDDUM*DBLE(SQRTW(J))                                          ETH03700
      DDDUM=DDDUM*DDUM                                                  ETH03800
      IF (DOYE) SYEI=SYEI+DUM*DBLE(Y(J))                                ETH03900
      IF (LDDUM) GO TO 120                                              ETH04000
      DUM=DUM*DBLE(SQRTW(J))                                            ETH04100
      SEI=SEI+DUM                                                       ETH04200
      IF (NOTINV) GO TO 120                                             ETH04300
      DUM=DUM*T(J)                                                      ETH04400
      DSEI=DSEI-DUM                                                     ETH04500
      DDSEI=DDSEI+DUM*T(J)                                              ETH04600
  120 CONTINUE                                                          ETH04700
  110 CONTINUE                                                          ETH04800
      IF (DOYE) SYE=SYEI                                                ETH04900
      IF (.NOT.DOE) RETURN                                              ETH05000
      IF (WTED) GO TO 300                                               ETH05100
C-----------------------------------------------------------------------ETH05200
C  USE GEOMETRIC SUM FORMULAS WHEN WTED=.FALSE.                         ETH05300
C-----------------------------------------------------------------------ETH05400
  200 DO 210 L=1,NINT                                                   ETH05500
      DL=PLAMBD*DELTAT(L)                                               ETH05600
      AN=DBLE(FLOAT(NT(L)))                                             ETH05700
      DNL=AN*DL                                                         ETH05800
      IF (DABS(DL) .LE. 1.D-3) GO TO 214                                ETH05900
      RDEDL=1.D0/(1.D0-DEXP(-DL))                                       ETH06000
      DEDNL=1.D0-DEXP(-DNL)                                             ETH06100
      GO TO 218                                                         ETH06200
  214 RDEDL=1.D0/EXPSMA(DL)                                             ETH06300
      IF (DABS(DNL) .LE. 1.D-3) DEDNL=EXPSMA(DNL)                       ETH06400
      IF (DABS(DNL) .GT. 1.D-3) DEDNL=1.D0-DEXP(-DNL)                   ETH06500
  218 ELTS=DEXP(-PLAMBD*TSTART(L))                                      ETH06600
      S=ELTS*DEDNL*RDEDL                                                ETH06700
      SEI=SEI+S                                                         ETH06800
      IF (NOTINV) GO TO 210                                             ETH06900
      B=DELTAT(L)*RDEDL                                                 ETH07000
      EDL=DEXP(-DL)                                                     ETH07100
      ENTS=ELTS*DEXP(-DNL)                                              ETH07200
      C=AN*ENTS-S*EDL                                                   ETH07300
      SS=B*C-TSTART(L)*S                                                ETH07400
      DSEI=DSEI+SS                                                      ETH07500
      DDSEI=DDSEI-TSTART(L)*SS-B*(EDL*(B*C+SS-S*DELTAT(L))+AN*(AN*      ETH07600
     1 DELTAT(L)+TSTART(L))*ENTS)                                       ETH07700
  210 CONTINUE                                                          ETH07800
  300 SE=SEI                                                            ETH07900
      IF (NOTINV) RETURN                                                ETH08000
      DSE=PLAMBD*DSEI                                                   ETH08100
      DDSE=DSE+DDSEI*PLAMBD*PLAMBD                                      ETH08200
      RETURN                                                            ETH08300
      END                                                               ETH08400
C--------------------- DOUBLE PRECISION VERSION 2B ---------------------PIV00100
C  SUBROUTINE PIVOT.  DOES GAUSS-JORDAN PIVOTS ON ALL NP DIAGONAL       PIV00200
C      ELEMENTS OF AUGMENTED MATRIX ADUB EXCEPT THOSE THAT HAVE A       PIV00300
C      TOLERANCE LESS THAN DTOLER OR THAT CAUSE A VIOLATION OF THE      PIV00400
C      CONSTRAINT PMIN .LE. P .LE. PMAX.                                PIV00500
C  ONLY NEEDS FULL UPPER TRIANGLE OF MATRIX ON INPUT, AND OVERWRITES IT PIV00600
C      WITH INVERSE AND SOLUTION.                                       PIV00700
C  NEEDS NO COMMON BLOCKS.                                              PIV00800
C-----------------------------------------------------------------------PIV00900
C  CALLS SUBPROGRAMS - PIVOT1                                           PIV01000
C-----------------------------------------------------------------------PIV01100
      SUBROUTINE PIVOT (ADUB,MDIMA,NP,ONLYIN,INVERT,CONSTR,PMIN,        PIV01200
     1 PMAX,P,DTOLER,DELTAP,ALLPIV,QG,PIV,JLAM,QMAX,IERROR)             PIV01300
      LOGICAL ONLYIN, ALLPIV, CONSTR, INVERT, PIV(NP)                   PIV01400
      DOUBLE PRECISION ADUB, PMIN, PMAX, P, DTOLER, DELTAP, DDUM, DUM,  PIV01500
     1 RDELTA                                                           PIV01600
      DIMENSION DELTAP(NP), P(NP), ADUB(MDIMA,MDIMA), DTOLER(NP)        PIV01700
      IERROR=0                                                          PIV01800
      QG=1.                                                             PIV01900
      QMAX=999.                                                         PIV02000
      NPP=NP+1                                                          PIV02100
      IF (ONLYIN) NPP=NP                                                PIV02200
      DO 110 J=1,NP                                                     PIV02300
  110 PIV(J)=.FALSE.                                                    PIV02400
      IF (NPP .EQ. 1) GO TO 190                                         PIV02500
      DO 120 J=2,NPP                                                    PIV02600
      L=J-1                                                             PIV02700
      DO 120 K=1,L                                                      PIV02800
      ADUB(J,K)=ADUB(K,J)                                               PIV02900
  120 CONTINUE                                                          PIV03000
C-----------------------------------------------------------------------PIV03100
C  MAIN LOOP FOR COMPLETE PIVOTING                                      PIV03200
C-----------------------------------------------------------------------PIV03300
  190 DO 200 NPIV=1,NP                                                  PIV03400
C-----------------------------------------------------------------------PIV03500
C  FIND NEXT PIVOT ELEMENT                                              PIV03600
C-----------------------------------------------------------------------PIV03700
      DUM=0.D0                                                          PIV03800
      DO 210 J=1,NP                                                     PIV03900
      IF (PIV(J) .OR. ADUB(J,J).LE.DTOLER(J)) GO TO 210                 PIV04000
      DDUM=ADUB(J,NPP)*ADUB(J,NPP)/ADUB(J,J)                            PIV04100
      IF (DDUM .LE. DUM) GO TO 210                                      PIV04200
      DUM=DDUM                                                          PIV04300
      L=J                                                               PIV04400
  210 CONTINUE                                                          PIV04500
      IF (DUM .LE. 0.D0) GO TO 300                                      PIV04600
      PIV(L)=.TRUE.                                                     PIV04700
      CALL PIVOT1 (ADUB,MDIMA,NPP,L)                                    PIV04800
  200 CONTINUE                                                          PIV04900
      NPIV=NP+1                                                         PIV05000
  300 NPIV=NPIV-1                                                       PIV05100
      IF (ONLYIN .OR. .NOT.CONSTR) GO TO 400                            PIV05200
C-----------------------------------------------------------------------PIV05300
C  ENFORCE CONSTRAINTS BY MAKING QG.LT.1. OR BY UNPIVOTING              PIV05400
C-----------------------------------------------------------------------PIV05500
      DO 330 J=1,JLAM                                                   PIV05600
      DUM=1.D+20                                                        PIV05700
      DO 340 K=1,JLAM                                                   PIV05800
      IF (.NOT.PIV(K)) GO TO 340                                        PIV05900
      RDELTA=1.D0/ADUB(K,NPP)                                           PIV06000
      DDUM=DMAX1((PMIN-P(K))*RDELTA,(PMAX-P(K))*RDELTA)                 PIV06100
      IF (DDUM .GE. DUM) GO TO 340                                      PIV06200
      DUM=DDUM                                                          PIV06300
      L=K                                                               PIV06400
  340 CONTINUE                                                          PIV06500
      IF (DUM .GT. 1.D-5) GO TO 350                                     PIV06600
      NPIV=NPIV-1                                                       PIV06700
      PIV(L)=.FALSE.                                                    PIV06800
      CALL PIVOT1 (ADUB,MDIMA,NPP,L)                                    PIV06900
  330 CONTINUE                                                          PIV07000
C-----------------------------------------------------------------------PIV07100
C  IF ALL ELEMENTS ARE SOMEHOW UNPIVOTED, TAKE ERROR EXIT               PIV07200
C-----------------------------------------------------------------------PIV07300
  350 IF (NPIV .GT. 0) GO TO 390                                        PIV07400
      IERROR=5                                                          PIV07500
      RETURN                                                            PIV07600
  390 QG=AMIN1(SNGL(DUM),1.)                                            PIV07700
      QMAX=SNGL(DUM)                                                    PIV07800
  400 ALLPIV=NPIV .EQ. NP                                               PIV07900
      IF (ONLYIN) RETURN                                                PIV08000
      DO 430 J=1,NP                                                     PIV08100
      DELTAP(J)=0.D0                                                    PIV08200
  430 IF (PIV(J)) DELTAP(J)=ADUB(J,NPP)                                 PIV08300
      RETURN                                                            PIV08400
      END                                                               PIV08500
C--------------------- DOUBLE PRECISION VERSION 2B ---------------------P1000100
C  SUBROUTINE PIVOT1. DOES A GAUSS-JORDAN PIVOT ON DIAGONAL ELEMENT LPIVP1000200
C  NEEDS NO COMMON BLOCKS.                                              P1000300
C-----------------------------------------------------------------------P1000400
C  CALLS NO SUBPROGRAMS                                                 P1000500
C-----------------------------------------------------------------------P1000600
      SUBROUTINE PIVOT1 (ADUB,MDIMA,NPP,LPIV)                           P1000700
      DOUBLE PRECISION ADUB, DUB, DDUB                                  P1000800
      DIMENSION ADUB(MDIMA,MDIMA)                                       P1000900
      DUB=1.D0/ADUB(LPIV,LPIV)                                          P1001000
      DO 110 J=1,NPP                                                    P1001100
      IF (J .EQ. LPIV) GO TO 110                                        P1001200
      DDUB=ADUB(J,LPIV)*DUB                                             P1001300
      DO 120 K=1,NPP                                                    P1001400
  120 IF (K .NE. LPIV) ADUB(J,K)=ADUB(J,K)-ADUB(LPIV,K)*DDUB            P1001500
  110 CONTINUE                                                          P1001600
      DO 130 J=1,NPP                                                    P1001700
      ADUB(J,LPIV)=-ADUB(J,LPIV)*DUB                                    P1001800
  130 ADUB(LPIV,J)=ADUB(LPIV,J)*DUB                                     P1001900
      ADUB(LPIV,LPIV)=DUB                                               P1002000
      RETURN                                                            P1002100
      END                                                               P1002200
C--------------------- DOUBLE PRECISION VERSION 2B ---------------------ANL00100
C  SUBROUTINE ANLERR.  DOES FINAL ERROR ANALYSIS OF FIT AND CALCULATES  ANL00200
C      NONLINEARITY TERM ENPHI OF BEALE FOR LATER USE.                  ANL00300
C-----------------------------------------------------------------------ANL00400
C  CALLS SUBPROGRAMS - PIVOT                                            ANL00500
C  WHICH IN TURN CALL - PIVOT1                                          ANL00600
C-----------------------------------------------------------------------ANL00700
      SUBROUTINE ANLERR (JLAM,T,SQRTW,NMAX,E,NLIN,NDIME,TSTART,DELTAT,  ANL00800
     1 NT,NINT,N,ITYPE,INTER,PRLSQ)                                     ANL00900
      LOGICAL LDUM, NOBASE, PIVALP, FAILED, PRLSQ, REGINT,              ANL01000
     1 LBLOKA, LBLOKB, INTER                                            ANL01100
      DOUBLE PRECISION SUMWT, VAR, PALPHA, PLAM, SE, DFCAPL, ADUB,      ANL01200
     1 DDSE, DTOLER, DBLOKA, DBLOKB, DBLOKC, DBLOKD, DBLOKE,            ANL01300
     2 E, TSTART, DELTAT, DUB, T                                        ANL01400
      COMMON /BLOK/ DBLOKA(4830), SUMWT, DBLOKB(12), VAR, PALPHA(10),   ANL01500
     1 PLAM(10), SE(10,10), DBLOKC(29), DFCAPL(10,9), DBLOKD(9),        ANL01600
     2 ADUB(19,19), DDSE(9,9), DTOLER(19), DBLOKE(380), RBLOKA(174),    ANL01700
     3 PRECIS, RBLOKB(4), SIGYY, ENPHI, SIGMAP(19), PCERR(19),          ANL01800
     4 RBLOKC(252), IBLOKA(199), MDIMA, IBLOKB, IWTSGN, IBLOKC(94),     ANL01900
     5 NF, IBLOKD, IWRITE, IBLOKE(6), JLAMP1, IBLOKF(9), REGINT         ANL02000
      COMMON /BLOK/ FAILED, LBLOKA, NOBASE, LBLOKB(23), PIVALP(19)      ANL02100
      DIMENSION CUI(19), CSTAR(19), DELT(18,18),                        ANL02200
     1 CCSUM(19,18), ADUM(10), ADSP(19,19)                              ANL02300
      DIMENSION SUM(4), T(NMAX), E(NDIME,NLIN), SQRTW(NMAX),            ANL02400
     1 TSTART(NINT), DELTAT(NINT), NT(NINT), HOLRTH(4,3)                ANL02500
      EQUIVALENCE (ADUB(1,1),ADSP(1,1))                                 ANL02600
      DATA HOLRTH/1HL, 1HA, 1HM, 1H , 1HA, 1HL, 1HP, 1H , 1HB, 1HA,     ANL02700
     1 1HS, 1HE/                                                        ANL02800
C-----------------------------------------------------------------------ANL02900
C  COMPUTE CORRELATION COEFFICIENTS AND STANDARD ERROR ESTIMATES        ANL03000
C      USING FULL NORMAL EQUATIONS MATRIX (ADUB) WITH ALL PARAMETERS    ANL03100
C      ASSUMED INDEPENDENT                                              ANL03200
C-----------------------------------------------------------------------ANL03300
      DO 110 K=1,JLAM                                                   ANL03400
      KK=K+JLAM                                                         ANL03500
      DO 120 J=1,K                                                      ANL03600
      JJ=J+JLAM                                                         ANL03700
      ADUB(J,K)=PALPHA(J)*PALPHA(K)*DDSE(J,K)                           ANL03800
      ADUB(J,KK)=-DFCAPL(K,J)*PALPHA(J)                                 ANL03900
      ADUB(K,JJ)=-DFCAPL(J,K)*PALPHA(K)                                 ANL04000
      ADUB(JJ,KK)=SE(J,K)                                               ANL04100
  120 CONTINUE                                                          ANL04200
      DTOLER(K)=ADUB(K,K)*PRECIS                                        ANL04300
      DTOLER(KK)=ADUB(KK,KK)*PRECIS                                     ANL04400
      IF (NOBASE) GO TO 110                                             ANL04500
      ADUB(K,NF)=-DFCAPL(NLIN,K)*PALPHA(K)                              ANL04600
      ADUB(KK,NF)=SE(K,NLIN)                                            ANL04700
  110 CONTINUE                                                          ANL04800
      IF (NOBASE) GO TO 150                                             ANL04900
      ADUB(NF,NF)=SUMWT                                                 ANL05000
      DTOLER(NF)=SUMWT*PRECIS                                           ANL05100
  150 CALL PIVOT (ADUB,MDIMA,NF,.TRUE.,.TRUE.,.FALSE.,0.D0,             ANL05200
     1 0.D0,DTOLER,DTOLER,DTOLER,LDUM,DUM,PIVALP,JLAM,DDUM,J)           ANL05300
      SIGYY=SQRT(SNGL(VAR)/FLOAT(N-NF))                                 ANL05400
      IF (.NOT.LDUM .OR. J.EQ.5) GO TO 172                              ANL05500
      DO 170 J=1,NF                                                     ANL05600
      IF (ADUB(J,J) .LE. 0.D0) GO TO 172                                ANL05700
      CSTAR(J)=1./SQRT(SNGL(ADUB(J,J)))                                 ANL05800
      SIGMAP(J)=SIGYY/CSTAR(J)                                          ANL05900
  170 CONTINUE                                                          ANL06000
      GO TO 174                                                         ANL06100
  172 FAILED=.TRUE.                                                     ANL06200
 5172 FORMAT (//51H0SINGULARITY IN INVERTING FULL LEAST SQUARES MATRIX, ANL06300
     1 30H - NO CORRELATIONS CALCULATED.,3X,24(2H**))                   ANL06400
      IF (PRLSQ) WRITE (IWRITE,5172)                                    ANL06500
      RETURN                                                            ANL06600
  174 IF (INTER) RETURN                                                 ANL06700
      K=JLAM                                                            ANL06800
      DO 176 J=1,JLAM                                                   ANL06900
      PCERR(J)=100.*SIGMAP(J)/SNGL(PLAM(J))                             ANL07000
      K=K+1                                                             ANL07100
      PCERR(K)=100.*SIGMAP(K)/ABS(SNGL(PALPHA(J)))                      ANL07200
  176 CONTINUE                                                          ANL07300
      IF (.NOT.NOBASE) PCERR(NF)=100.*SIGMAP(NF)/ABS(SNGL(PALPHA(       ANL07400
     1 JLAMP1)))                                                        ANL07500
 5176 FORMAT (//10X,24HCORRELATION COEFFICIENTS/5X,18(3X,3A1,I1))       ANL07600
      IF (.NOT.PRLSQ) GO TO 190                                         ANL07700
      K=JLAM                                                            ANL07800
      IF (NOBASE) K=K-1                                                 ANL07900
      IF (K .GT. 0) WRITE (IWRITE,5176) ((HOLRTH(J,1),J=1,3),L,L=1,JLAM)ANL08000
     1 ,((HOLRTH(JJ,2),JJ=1,3),JK,JK=1,K)                               ANL08100
      IF (K .EQ. 0) WRITE (IWRITE,5176) ((HOLRTH(J,1),J=1,3),L,L=1,JLAM)ANL08200
      DO 179 J=2,NF                                                     ANL08300
      DUM=CSTAR(J)                                                      ANL08400
      JK=J-1                                                            ANL08500
      DO 180 K=1,JK                                                     ANL08600
  180 DTOLER(K)=ADUB(J,K)*DUM*CSTAR(K)                                  ANL08700
 5180 FORMAT (1X,3A1,I1,18F7.3)                                         ANL08800
      L=(J-1)/JLAM+1                                                    ANL08900
      JJ=J                                                              ANL09000
      IF (L .EQ. 2) JJ=JJ-JLAM                                          ANL09100
      IF (L .LT. 3) WRITE(IWRITE,5180) (HOLRTH(KK,L),KK=1,3),JJ,        ANL09200
     1 (DTOLER(K),K=1,JK)                                               ANL09300
  179 CONTINUE                                                          ANL09400
 5181 FORMAT (1X,4A1,18F7.3)                                            ANL09500
      IF (L .EQ. 3) WRITE (IWRITE,5181) (HOLRTH(KK,3),KK=1,4),          ANL09600
     1 (DTOLER(K),K=1,JK)                                               ANL09700
C-----------------------------------------------------------------------ANL09800
C  TO PREVENT A VERY SMALL COMPONENT (EFFECTIVELY A SECOND BASELINE)    ANL09900
C      THAT IS HIGHLY CORRELATED WITH THE BASELINE FROM INTERFERING     ANL10000
C      WITH THE CALCULATION OF THE WEIGHTS                              ANL10100
C-----------------------------------------------------------------------ANL10200
  190 IF (IWTSGN.GT.0 .OR. ITYPE.NE.1) GO TO 200                        ANL10300
      DUB=1.D+20                                                        ANL10400
      DO 195 J=1,JLAM                                                   ANL10500
      IF (PLAM(J) .GT. DUB) GO TO 195                                   ANL10600
      DUB=PLAM(J)                                                       ANL10700
      L=J                                                               ANL10800
  195 CONTINUE                                                          ANL10900
      IF (PCERR(L) .LT. 1000.) GO TO 200                                ANL11000
      FAILED=.TRUE.                                                     ANL11100
 5195 FORMAT (//47H SMALLEST LAMBDA HAS TOO HIGH A STANDARD ERROR.,     ANL11200
     1 60H  THIS COULD LEAD TO INACCURATE WEIGHTS - SOLUTION REJECTED., ANL11300
     2 3X,11(2H**))                                                     ANL11400
      IF (PRLSQ) WRITE (IWRITE,5195)                                    ANL11500
  200 IF (FAILED .OR. ITYPE.NE.4) GO TO 282                             ANL11600
C-----------------------------------------------------------------------ANL11700
C  CALCULATION OF ENPHI (EQ. 3.25 OF E. M. L. BEALE FOR F RATIO         ANL11800
C  CORRECTION FACTOR DUE TO NONLINEARITY)                               ANL11900
C-----------------------------------------------------------------------ANL12000
      DO 205 K=1,NF                                                     ANL12100
      DO 205 J=1,NF                                                     ANL12200
  205 ADSP(J,K)=SNGL(ADUB(J,K))                                         ANL12300
      NUJ=JLAM                                                          ANL12400
      DO 210 J=1,JLAM                                                   ANL12500
      DUM=ADSP(J,J)                                                     ANL12600
      NUJ=NUJ+1                                                         ANL12700
      NUK=JLAM                                                          ANL12800
      DO 220 K=1,JLAM                                                   ANL12900
      NUK=NUK+1                                                         ANL13000
      DELT(J,K)=DUM*ADSP(K,K)+2.*ADSP(J,K)*ADSP(J,K)                    ANL13100
      DELT(J,NUK)=DUM*ADSP(K,NUK)+2.*ADSP(J,K)*ADSP(J,NUK)              ANL13200
      DELT(NUJ,K)=ADSP(J,NUJ)*ADSP(K,K)+2.*ADSP(J,K)*ADSP(NUJ,K)        ANL13300
      DELT(NUJ,NUK)=ADSP(J,NUJ)*ADSP(K,NUK)+ADSP(J,K)*ADSP(NUJ,NUK)     ANL13400
     1 +ADSP(J,NUK)*ADSP(NUJ,K)                                         ANL13500
      DO 230 L=1,4                                                      ANL13600
  230 SUM(L)=0.                                                         ANL13700
      NU=0                                                              ANL13800
      DO 240 L=1,NINT                                                   ANL13900
      CUI(1)=EXP(-SNGL(PLAM(J)*DELTAT(L)))                              ANL14000
      CUI(2)=EXP(-SNGL(PLAM(K)*DELTAT(L)))                              ANL14100
      CUI(3)=EXP(-SNGL(PLAM(J)*TSTART(L)))                              ANL14200
      CUI(4)=EXP(-SNGL(PLAM(K)*TSTART(L)))                              ANL14300
      II=NT(L)                                                          ANL14400
      DO 242 I=1,II                                                     ANL14500
      NU=NU+1                                                           ANL14600
      TNU=SNGL(T(NU))                                                   ANL14700
      IF (REGINT) GO TO 244                                             ANL14800
      SS=SNGL(E(NU,J))                                                  ANL14900
      DDDUM=TNU*SNGL(E(NU,K))                                           ANL15000
      GO TO 246                                                         ANL15100
  244 SS=SQRTW(NU)*CUI(3)                                               ANL15200
      DDDUM=TNU*SQRTW(NU)*CUI(4)                                        ANL15300
      CUI(3)=CUI(3)*CUI(1)                                              ANL15400
      CUI(4)=CUI(4)*CUI(2)                                              ANL15500
  246 S=SNGL(PALPHA(J))*TNU*SS                                          ANL15600
      DDUM=SNGL(PALPHA(K))*TNU*DDDUM                                    ANL15700
      SUM(1)=SUM(1)-S*DDUM                                              ANL15800
      SUM(2)=SUM(2)+S*DDDUM                                             ANL15900
      SUM(3)=SUM(3)+SS*DDUM                                             ANL16000
      SUM(4)=SUM(4)-SS*DDDUM                                            ANL16100
  242 CONTINUE                                                          ANL16200
  240 CONTINUE                                                          ANL16300
      CCSUM(J,K)=SUM(1)                                                 ANL16400
      CCSUM(J,NUK)=SUM(2)                                               ANL16500
      CCSUM(NUJ,K)=SUM(3)                                               ANL16600
      CCSUM(NUJ,NUK)=SUM(4)                                             ANL16700
  220 CONTINUE                                                          ANL16800
      IF (NOBASE) GO TO 210                                             ANL16900
      S=0.                                                              ANL17000
      SS=0.                                                             ANL17100
      NU=0                                                              ANL17200
      DO 250 L=1,NINT                                                   ANL17300
      DUM=EXP(-SNGL(PLAM(J)*DELTAT(L)))                                 ANL17400
      DDUM=EXP(-SNGL(PLAM(J)*TSTART(L)))                                ANL17500
      II=NT(L)                                                          ANL17600
      DO 252 I=1,II                                                     ANL17700
      NU=NU+1                                                           ANL17800
      TNU=SNGL(T(NU))                                                   ANL17900
      IF (REGINT) GO TO 254                                             ANL18000
      SSS=SNGL(E(NU,J))                                                 ANL18100
      GO TO 256                                                         ANL18200
  254 SSS=DDUM*SQRTW(NU)                                                ANL18300
      DDUM=DDUM*DUM                                                     ANL18400
  256 DDDUM=TNU*SSS*SQRTW(NU)                                           ANL18500
      SS=SS-DDDUM                                                       ANL18600
      S=S+SNGL(PALPHA(J))*TNU*DDDUM                                     ANL18700
  252 CONTINUE                                                          ANL18800
  250 CONTINUE                                                          ANL18900
      CCSUM(NF,J)=S                                                     ANL19000
      CCSUM(NF,NUJ)=SS                                                  ANL19100
  210 CONTINUE                                                          ANL19200
      LDUM=.NOT.NOBASE                                                  ANL19300
      ENPHI=0.                                                          ANL19400
      DO 260 NU=1,N                                                     ANL19500
      TNU=SNGL(T(NU))                                                   ANL19600
      DO 261 J=1,NLIN                                                   ANL19700
  261 ADUM(J)=SQRTW(NU)*EXP(-SNGL(PLAM(J))*TNU)                         ANL19800
      NUJ=JLAM                                                          ANL19900
      DO 262 J=1,JLAM                                                   ANL20000
      NUJ=NUJ+1                                                         ANL20100
      SUM(1)=0.                                                         ANL20200
      SUM(2)=0.                                                         ANL20300
      NUK=JLAM                                                          ANL20400
      DO 264 K=1,JLAM                                                   ANL20500
      S=0.                                                              ANL20600
      IF (LDUM) S=SQRTW(NU)*ADSP(K,NF)                                  ANL20700
      NUI=JLAM                                                          ANL20800
      DO 266 I=1,JLAM                                                   ANL20900
      NUI=NUI+1                                                         ANL21000
      CUI(NUI)=ADUM(I)                                                  ANL21100
      CUI(I)=-SNGL(PALPHA(I))*TNU*CUI(NUI)                              ANL21200
      S=S+ADSP(K,I)*CUI(I)+ADSP(K,NUI)*CUI(NUI)                         ANL21300
  266 CONTINUE                                                          ANL21400
      SUM(1)=SUM(1)+CCSUM(K,J)*S                                        ANL21500
      SUM(2)=SUM(2)+CCSUM(K,NUJ)*S                                      ANL21600
      NUK=NUK+1                                                         ANL21700
      S=0.                                                              ANL21800
      IF (LDUM) S=SQRTW(NU)*ADSP(NUK,NF)                                ANL21900
      NUI=JLAM                                                          ANL22000
      DO 268 I=1,JLAM                                                   ANL22100
      NUI=NUI+1                                                         ANL22200
      S=S+ADSP(NUK,I)*CUI(I)+ADSP(NUK,NUI)*CUI(NUI)                     ANL22300
  268 CONTINUE                                                          ANL22400
      SUM(1)=SUM(1)+CCSUM(NUK,J)*S                                      ANL22500
      SUM(2)=SUM(2)+CCSUM(NUK,NUJ)*S                                    ANL22600
  264 CONTINUE                                                          ANL22700
      IF (NOBASE) GO TO 272                                             ANL22800
      S=0.                                                              ANL22900
      IF (LDUM) S=SQRTW(NU)*ADSP(NF,NF)                                 ANL23000
      NUI=JLAM                                                          ANL23100
      DO 270 I=1,JLAM                                                   ANL23200
      NUI=NUI+1                                                         ANL23300
      S=S+ADSP(NF,I)*CUI(I)+ADSP(NF,NUI)*CUI(NUI)                       ANL23400
  270 CONTINUE                                                          ANL23500
      SUM(1)=SUM(1)+CCSUM(NF,J)*S                                       ANL23600
      SUM(2)=SUM(2)+CCSUM(NF,NUJ)*S                                     ANL23700
  272 DUM=TNU*ADUM(J)                                                   ANL23800
      CSTAR(J)=SNGL(PALPHA(J))*TNU*DUM-SUM(1)                           ANL23900
      CSTAR(NUJ)=-2.*(DUM+SUM(2))                                       ANL24000
  262 CONTINUE                                                          ANL24100
      SUM(1)=0.                                                         ANL24200
      NUJ=JLAM                                                          ANL24300
      DO 280 J=1,JLAM                                                   ANL24400
      NUJ=NUJ+1                                                         ANL24500
      S=0.                                                              ANL24600
      SS=0.                                                             ANL24700
      NUK=JLAM                                                          ANL24800
      DO 281 K=1,JLAM                                                   ANL24900
      NUK=NUK+1                                                         ANL25000
      S=S+CSTAR(K)*DELT(J,K)+CSTAR(NUK)*DELT(J,NUK)                     ANL25100
      SS=SS+CSTAR(K)*DELT(NUJ,K)+CSTAR(NUK)*DELT(NUJ,NUK)               ANL25200
  281 CONTINUE                                                          ANL25300
      SUM(1)=SUM(1)+CSTAR(J)*S+CSTAR(NUJ)*SS                            ANL25400
  280 CONTINUE                                                          ANL25500
      ENPHI=ENPHI+SUM(1)                                                ANL25600
  260 CONTINUE                                                          ANL25700
      ENPHI=ENPHI*SIGYY**2/FLOAT(NF+2)                                  ANL25800
 5282 FORMAT (36H STANDARD DEVIATION OF FIT = SIGYY =,1PE12.5//         ANL25900
     1 14X,32HLAMBDA  +-  STD. ERROR   PERCENT,15X,                     ANL26000
     2 31HALPHA  +-  STD. ERROR   PERCENT)                              ANL26100
 5002 FORMAT (/1X)                                                      ANL26200
  282 IF (.NOT.PRLSQ) RETURN                                            ANL26300
      WRITE (IWRITE,5002)                                               ANL26400
      IF (.NOT. (FAILED .OR. ITYPE.NE.4)) WRITE (IWRITE,5283) ENPHI     ANL26500
 5283 FORMAT ( 7H NPHI =,1PE10.3)                                       ANL26600
      WRITE (IWRITE,5282) SIGYY                                         ANL26700
      K=JLAM                                                            ANL26800
      DO 284 J=1,JLAM                                                   ANL26900
      K=K+1                                                             ANL27000
 5284 FORMAT (2(1PD20.4,4H  +-,E12.3,0PF10.3))                          ANL27100
      WRITE (IWRITE,5284) PLAM(J),SIGMAP(J),PCERR(J),PALPHA(J),         ANL27200
     1 SIGMAP(K),PCERR(K)                                               ANL27300
  284 CONTINUE                                                          ANL27400
 5286 FORMAT (12X,8HBASELINE,26X,1PD20.4,4H  +-,E12.3,0PF10.3)          ANL27500
      IF (.NOT.NOBASE) WRITE(IWRITE,5286) PALPHA(JLAMP1),SIGMAP(NF),    ANL27600
     1 PCERR(NF)                                                        ANL27700
      RETURN                                                            ANL27800
      END                                                               ANL27900
C-------------------------- VERSION 2A OR 2B ---------------------------FSH00100
C  FUNCTION FISHER.  CALCULATES PROBABILITY THAT FS(NU1,NU2) IS LESS    FSH00200
C  THAN F FOR FISHER F-TEST.                                            FSH00300
C  NU2 MUST BE .GE. 3                                                   FSH00400
C  SUMS UNTIL A TERM IS FOUND LESS THAN TOL*(CURRENT VALUE OF SUM).     FSH00500
C  EQUATION NUMBERS REFER TO ABRAMOWITZ AND STEGUN, P. 946.             FSH00600
C  NO COMMON BLOCKS NEEDED.                                             FSH00700
C-----------------------------------------------------------------------FSH00800
C  CALLS NO SUBPROGRAMS.                                                FSH00900
C-----------------------------------------------------------------------FSH01000
      FUNCTION FISHER (FS,NU1,NU2)                                      FSH01100
      DATA TOL/1.E-9/                                                   FSH01200
      RNU1=FLOAT(NU1)                                                   FSH01300
      RNU2=FLOAT(NU2)                                                   FSH01400
      X=RNU2/(RNU2+RNU1*FS)                                             FSH01500
      IF (X.LT.1. .OR. FS.LT.0.) GO TO 105                              FSH01600
      FISHER=0.                                                         FSH01700
      RETURN                                                            FSH01800
  105 IF (MOD(NU1,2) .NE. 0) GO TO 200                                  FSH01900
C-----------------------------------------------------------------------FSH02000
C  EVALUATES EQ. 26.6.4                                                 FSH02100
C-----------------------------------------------------------------------FSH02200
      Q=1.                                                              FSH02300
      IF (NU1 .EQ. 2) GO TO 150                                         FSH02400
      TERM=1.                                                           FSH02500
      TX=1.-X                                                           FSH02600
      TN=RNU2-2.                                                        FSH02700
      K=NU1-2                                                           FSH02800
      DO 110 J=2,K,2                                                    FSH02900
      TN=TN+2.                                                          FSH03000
      TERM=TERM*TN*TX/FLOAT(J)                                          FSH03100
      Q=Q+TERM                                                          FSH03200
      IF (TERM .LE. TOL*Q) GO TO 150                                    FSH03300
  110 CONTINUE                                                          FSH03400
  150 FISHER=1.-X**(.5*RNU2)*Q                                          FSH03500
      RETURN                                                            FSH03600
  200 IF (MOD(NU2,2) .EQ. 0) GO TO 300                                  FSH03700
C-----------------------------------------------------------------------FSH03800
C  EVALUATES EQ. 26.6.8                                                 FSH03900
C-----------------------------------------------------------------------FSH04000
      TH=ATAN(SQRT(RNU1*FS/RNU2))                                       FSH04100
      CTH=COS(TH)                                                       FSH04200
      CC=CTH*CTH                                                        FSH04300
      STH=SIN(TH)                                                       FSH04400
      SS=STH*STH                                                        FSH04500
      TERM=1.                                                           FSH04600
      Q=1.                                                              FSH04700
      K=NU2-2                                                           FSH04800
      PR=1.                                                             FSH04900
      TN=0.                                                             FSH05000
      IF (K .LT. 3) GO TO 220                                           FSH05100
      DO 210 J=3,K,2                                                    FSH05200
      TN=TN+2.                                                          FSH05300
      DUM=TN/FLOAT(J)                                                   FSH05400
      PR=PR*DUM                                                         FSH05500
      TERM=TERM*CC*DUM                                                  FSH05600
      Q=Q+TERM                                                          FSH05700
      IF (TERM .LE. TOL*Q) GO TO 220                                    FSH05800
  210 CONTINUE                                                          FSH05900
      J=K                                                               FSH06000
  220 FISHER=TH+STH*CTH*Q                                               FSH06100
      Q=0.                                                              FSH06200
      IF (NU1 .EQ. 1) GO TO 290                                         FSH06300
      TN=TN+2.                                                          FSH06400
      PR=PR*TN                                                          FSH06500
      IF (J .EQ. K) GO TO 240                                           FSH06600
      L=J+2                                                             FSH06700
      DO 230 J=L,K,2                                                    FSH06800
      TN=TN+2.                                                          FSH06900
  230 PR=PR*TN/FLOAT(J)                                                 FSH07000
  240 TN=RNU2-1.                                                        FSH07100
      Q=1.                                                              FSH07200
      TERM=1.                                                           FSH07300
      K=NU1-2                                                           FSH07400
      IF (K .LT. 3) GO TO 260                                           FSH07500
      DO 250 J=3,K,2                                                    FSH07600
      TN=TN+2.                                                          FSH07700
      TERM=TERM*TN*SS/FLOAT(J)                                          FSH07800
      Q=Q+TERM                                                          FSH07900
      IF (TERM .LE. TOL*Q) GO TO 260                                    FSH08000
  250 CONTINUE                                                          FSH08100
  260 FISHER=FISHER-Q*PR*STH*CTH**NU2                                   FSH08200
  290 FISHER=FISHER*.6366198                                            FSH08300
      RETURN                                                            FSH08400
C-----------------------------------------------------------------------FSH08500
C  EVALUATES EQ. 26.6.5                                                 FSH08600
C-----------------------------------------------------------------------FSH08700
  300 TERM=1.                                                           FSH08800
      Q=1.                                                              FSH08900
      K=NU2-2                                                           FSH09000
      TN=RNU1-2.                                                        FSH09100
      DO 310 J=2,K,2                                                    FSH09200
      TN=TN+2.                                                          FSH09300
      TERM=TERM*TN*X/FLOAT(J)                                           FSH09400
      Q=Q+TERM                                                          FSH09500
      IF (TERM .LE. TOL*Q) GO TO 320                                    FSH09600
  310 CONTINUE                                                          FSH09700
  320 FISHER=SQRT(1.-X)*Q                                               FSH09800
      RETURN                                                            FSH09900
      END                                                               FSH10000
C--------------------- DOUBLE PRECISION VERSION 2B ---------------------RES00100
C  SUBROUTINE RESIDU. CALCULATES WEIGHTED RESIDUALS FROM ASAVE(J,JLAM,1)RES00200
C      AND ASAVE(J,JLAM,2), WHICH CORRESPOND TO ALPHA AND LAMBDA,       RES00300
C      RETURNS PUNCOR (APPROX. PROBABILITY THAT THE RESIDUALS ARE       RES00400
C      UNCORRELATED FOR LAGS 1 THRU 5) AND WEIGHTED AVERAGE IN          RES00500
C      PUNCOR(6), AND, IF PLOT=.TRUE., CALLS PLPRES TO PLOT WEIGHTED    RES00600
C      RESIDUALS ON PRINTER.                                            RES00700
C-----------------------------------------------------------------------RES00800
C  CALLS SUBPROGRAMS - FISHER, PLPRES                                   RES00900
C-----------------------------------------------------------------------RES01000
      SUBROUTINE RESIDU (JLAM,PLOT,T,Y,SQRTW,YLYFIT,NMAX,PUNCOR,ASAVE)  RES01100
      LOGICAL PLOT, NOBASE, WTED, LDUM, LBLOKA, LBLOKB                  RES01200
      DOUBLE PRECISION DBLOKA, YLYFIT, RSUMST, RSUMEN, STBAR, ENBAR,    RES01300
     1 SS, SE, EE, S, E, T                                              RES01400
      COMMON /BLOK/ DBLOKA(5933), RBLOKA(471), IBLOKA(297), N, IWRITE,  RES01500
     1 IBLOKB, LINEPG, IBLOKC(14), LBLOKA(3), NOBASE, WTED, LBLOKB(41)  RES01600
      DIMENSION ASAVE(10,9,6)                                           RES01700
      DIMENSION T(NMAX), Y(NMAX), SQRTW(NMAX), YLYFIT(NMAX), PUNCOR(6)  RES01800
      JLAMP1=JLAM+1                                                     RES01900
      LDUM=.NOT.NOBASE                                                  RES02000
      IF (LDUM) DUM=ASAVE(JLAMP1,JLAM,1)                                RES02100
      DO 110 K=1,N                                                      RES02200
      YLYFIT(K)=DBLE(Y(K))                                              RES02300
  110 IF (LDUM) YLYFIT(K)=YLYFIT(K)-DBLE(DUM*SQRTW(K))                  RES02400
      DO 120 J=1,JLAM                                                   RES02500
      DUM=ASAVE(J,JLAM,1)                                               RES02600
      DDUM=ASAVE(J,JLAM,2)                                              RES02700
      DO 125 K=1,N                                                      RES02800
  125 YLYFIT(K)=YLYFIT(K)-DBLE(DUM*EXP(-DDUM*SNGL(T(K)))*SQRTW(K))      RES02900
  120 CONTINUE                                                          RES03000
      RSUMST=0.D0                                                       RES03100
      IF (LDUM) GO TO 140                                               RES03200
      DO 130 K=1,N                                                      RES03300
  130 RSUMST=RSUMST+YLYFIT(K)                                           RES03400
  140 RSUMEN=RSUMST                                                     RES03500
C-----------------------------------------------------------------------RES03600
C  CALCULATE AUTOCOVARIANCES AND PROBABILITIES THAT RESIDUALS ARE       RES03700
C  UNCORRELATED (HAMILTON, P. 183)                                      RES03800
C-----------------------------------------------------------------------RES03900
      PUNCOR(6)=0.                                                      RES04000
      DO 150 L=1,5                                                      RES04100
      NEND=N-L                                                          RES04200
      RSUMST=RSUMST-YLYFIT(NEND+1)                                      RES04300
      RSUMEN=RSUMEN-YLYFIT(L)                                           RES04400
      STBAR=RSUMST/FLOAT(NEND)                                          RES04500
      ENBAR=RSUMEN/FLOAT(NEND)                                          RES04600
      SS=0.D0                                                           RES04700
      SE=0.D0                                                           RES04800
      EE=0.D0                                                           RES04900
      J=L                                                               RES05000
      DO 152 K=1,NEND                                                   RES05100
      S=YLYFIT(K)-STBAR                                                 RES05200
      J=J+1                                                             RES05300
      E=YLYFIT(J)-ENBAR                                                 RES05400
      SS=SS+S*S                                                         RES05500
      SE=SE+S*E                                                         RES05600
      EE=EE+E*E                                                         RES05700
  152 CONTINUE                                                          RES05800
      RSQ=SNGL(SE*SE/(SS*EE))                                           RES05900
      PUNCOR(L)=1.-FISHER(FLOAT(NEND-2)*RSQ/(1.-RSQ),1,NEND-2)          RES06000
      PUNCOR(6)=PUNCOR(6)+PUNCOR(L)/FLOAT(L)                            RES06100
  150 CONTINUE                                                          RES06200
      PUNCOR(6)=PUNCOR(6)/2.283333                                      RES06300
      IF (PLOT) CALL PLPRES (YLYFIT,NMAX,N,WTED,IWRITE,LINEPG)          RES06400
      RETURN                                                            RES06500
      END                                                               RES06600
C--------------------- DOUBLE PRECISION VERSION 2B ---------------------PLP00100
C  SUBROUTINE PLPRES. PLOTS RESIDUALS ON PRINTER WITH X-AXIS HORIZONTAL.PLP00200
C  NEEDS NO COMMON BLOCKS                                               PLP00300
C-----------------------------------------------------------------------PLP00400
C  CALLS NO EXTERNAL SUBPROGRAMS.                                       PLP00500
C-----------------------------------------------------------------------PLP00600
      SUBROUTINE PLPRES (YLYFIT,NMAX,N,WTED,IWRITE,LINEPG)              PLP00700
      LOGICAL WTED                                                      PLP00800
      DOUBLE PRECISION RMIN, RMAX, DELTA, BOUND, BMAX, BMIN, YLYFIT     PLP00900
      DIMENSION YLYFIT(NMAX), JCHAR(8), LINE(131), LABEL(6), BOUND(21), PLP01000
     1 LCHARJ(20), LINE1(20)                                            PLP01100
      DATA JCHAR/1H*, 1H-, 1HU, 1HL, 1H , 1H0, 1H=, 1H+/                PLP01200
      IF (LINEPG .GE. 17) GO TO 100                                     PLP01300
 5100 FORMAT (/////47H NO PLOT OF RESIDUALS POSSIBLE BECAUSE LINEPG =,  PLP01400
     1 I2,2X,40(2H**)//)                                                PLP01500
      WRITE (IWRITE,5100) LINEPG                                        PLP01600
      RETURN                                                            PLP01700
  100 MPLOT=(LINEPG-1)/16                                               PLP01800
      MLINE=MIN0(20,(LINEPG-1)/MPLOT-3)                                 PLP01900
      RMIN=1.D+20                                                       PLP02000
      RMAX=-1.D+20                                                      PLP02100
      DO 110 J=1,N                                                      PLP02200
      RMIN=DMIN1(RMIN,YLYFIT(J))                                        PLP02300
  110 RMAX=DMAX1(RMAX,YLYFIT(J))                                        PLP02400
 5200 FORMAT(55H1PLOT OF RESIDUALS MULTIPLIED BY SQUARE-ROOT OF WEIGHTS PLP02500
     1 ,21H - FOR BEST SOLUTION.,4X,15HMAX RESIDUAL=U=,1PD8.1,4X,       PLP02600
     2 15HMIN RESIDUAL=L=,1PD8.1)                                       PLP02700
      IF (WTED) WRITE (IWRITE,5200) RMAX,RMIN                           PLP02800
 5202 FORMAT (18H1PLOT OF RESIDUALS                                     PLP02900
     1 ,21H - FOR BEST SOLUTION.,4X,15HMAX RESIDUAL=U=,1PD8.1,4X,       PLP03000
     2 15HMIN RESIDUAL=L=,1PD8.1)                                       PLP03100
      IF (.NOT.WTED) WRITE (IWRITE,5202) RMAX,RMIN                      PLP03200
      DELTA=(RMAX-RMIN)/FLOAT(MLINE-1)                                  PLP03300
      BOUND(1)=RMAX+.5*DELTA                                            PLP03400
      K=MLINE+1                                                         PLP03500
      DO 120 J=2,K                                                      PLP03600
      BOUND(J)=BOUND(J-1)-DELTA                                         PLP03700
  120 IF (BOUND(J)*BOUND(J-1) .LT. 0.) JAXIS=J-1                        PLP03800
      LABEL(1)=-110                                                     PLP03900
      DO 130 J=2,6                                                      PLP04000
  130 LABEL(J)=LABEL(J-1)+20                                            PLP04100
      K=MLINE-1                                                         PLP04200
      DO 140 J=2,K                                                      PLP04300
      LINE1(J)=JCHAR(5)                                                 PLP04400
  140 LCHARJ(J)=5                                                       PLP04500
      LINE1(1)=JCHAR(3)                                                 PLP04600
      LINE1(JAXIS)=JCHAR(6)                                             PLP04700
      LINE1(MLINE)=JCHAR(4)                                             PLP04800
      LCHARJ(1)=7                                                       PLP04900
      LCHARJ(JAXIS)=2                                                   PLP05000
      LCHARJ(MLINE)=7                                                   PLP05100
      NPOINT=0                                                          PLP05200
      DO 200 NPAGE=1,30                                                 PLP05300
      IF (NPAGE .GT. 1) WRITE(IWRITE,5999)                              PLP05400
 5999 FORMAT (1H1)                                                      PLP05500
      DO 210 NPLOT=1,MPLOT                                              PLP05600
      NST=NPOINT+1                                                      PLP05700
      NEND=NPOINT+130                                                   PLP05800
      NPOINT=NEND                                                       PLP05900
      NLIM=MIN0(NEND,N)                                                 PLP06000
      DO 220 NLINE=1,MLINE                                              PLP06100
      LCHAR=LCHARJ(NLINE)                                               PLP06200
      LINE(1)=LINE1(NLINE)                                              PLP06300
      BMAX=BOUND(NLINE)                                                 PLP06400
      BMIN=BOUND(NLINE+1)                                               PLP06500
      K=1                                                               PLP06600
      DO 230 J=NST,NLIM                                                 PLP06700
      K=K+1                                                             PLP06800
      LINE(K)=JCHAR(LCHAR)                                              PLP06900
      IF (YLYFIT(J).LT.BMAX .AND. YLYFIT(J).GE.BMIN) LINE(K)=JCHAR(1)   PLP07000
  230 CONTINUE                                                          PLP07100
      K=NLIM-NST+2                                                      PLP07200
      IF (NLINE .NE. MLINE) GO TO 235                                   PLP07300
      DO 232 J=21,K,20                                                  PLP07400
  232 IF (LINE(J) .NE. JCHAR(1)) LINE(J)=JCHAR(8)                       PLP07500
 5001 FORMAT (1X)                                                       PLP07600
  235 IF (NLINE .EQ. 1) WRITE (IWRITE,5001)                             PLP07700
 5230 FORMAT (1X,131A1)                                                 PLP07800
      WRITE (IWRITE,5230) (LINE(J),J=1,K)                               PLP07900
  220 CONTINUE                                                          PLP08000
      DO 240 J=1,6                                                      PLP08100
  240 LABEL(J)=LABEL(J)+130                                             PLP08200
 5240 FORMAT (3X,6(16X,I4)/)                                            PLP08300
      WRITE (IWRITE,5240) LABEL                                         PLP08400
      IF (NLIM .EQ. N) RETURN                                           PLP08500
  210 CONTINUE                                                          PLP08600
  200 CONTINUE                                                          PLP08700
 5250 FORMAT (49H NO MORE THAN 30 PAGES OF RESIDUAL PLOTS ALLOWED.,     PLP08800
     1 2X,38(2H**))                                                     PLP08900
      WRITE (IWRITE,5250)                                               PLP09000
      RETURN                                                            PLP09100
      END                                                               PLP09200
