C
C----------------MAIN--EXAFIT--PROGRAM--------------------------------
C
C * PROGRAM TO FIT MODEL EXAFS FUNCTION CHI TO ACTUAL EXAFS DATA
C * WITH ALMOST ALL COMBINATIONS OF PARAMETERS TO VARY IN A NON-LINEAR
C * LEAST SQUARES ITERATING FIT
      PROGRAM EXAFIT
      DIMENSION CHI(1024),CAY(1024),Y(1024),CAYO(1024),RHS(111),
     C    IP(111),ISH(111),SS(2,12,10),S(2)
C	COMMON /TKTRNX/ TMINVX,TMINVY,TMAXVX,TMAXVY,TREALX,TREALY,
C     1	TIMAGX,TIMAGY,TRCOSF,TRSINF,TRSCAL,TRFACX,TRFACY,
C     2	TRPAR1,TRPAR2,TRPAR3,TRPAR4,TRPAR5,TRPAR6,KMOFLG(2),
C     -	KGNMOD,KPADV,KACHAR,KOBLEN,KTRAIL,KLEVEL,KPAD2,
C     3	KBAUDR,KGNFLG,KGRAFL,KHOMEY,KKMODE,KHORSZ,KVERSZ,KTBLSZ,
C     4	KSIZEF,KLMRGN,KRMRGN,KFACTR,KTERM,KLINE,KZAXIS,KBEAMX,KBEAMY,
C     5	KMOVEF,KPCHAR(5),KDASHT,KMINSX,KMINSY,KMAXSX,KMAXSY,KEYCON,
C     6	KINLFT,KOTLFT,KUNIT
      COMMON /CARD1/ MEMBER(10),M,NDS,NITER,NSTEP,DAMPER,DESIRE,
     C    Q,POWER,KPOWER,TFLAG,CHKPRT,CHKOUT
      LOGICAL*1 CHANGE,MEMBER,TFLAG,CHKPRT,CHKOUT,SS,S,YN
      REAL*8 CHI,RHS
      NITER=0
      IREAD=0
10     DO 30 I=1,111
   30 RHS(I)=0.0D0
      RTURN1=0.
      MP=0
      IERR=0
      CHANGE=.TRUE.
      IF(IREAD.NE.2)NPTS=0
      EZERO=0.0
      VPREV=1.E10
      CALL INPUT(CAY,Y,CAYO,NPTS,IP,ISH,MP,EZERO,SS,S,IREAD,IERR)
      IF(IERR.EQ.1)GO TO 10
      CALL OUTPUT(VPREV,SS,S)
      DO 40 N=1,NITER
      NAN=N
      CALL ITER8(CHI,CAY,Y,CAYO,NPTS,NAN,RHS,MP,IP,ISH,SS,S,
     C    VPREV,CHANGE,EZERO,RTURN1)
      IF(NITER.EQ.0) GO TO 50
      IF(RTURN1.EQ.1.)GO TO 50
      IF(.NOT.CHANGE) GO TO 50
      IF(VPREV.EQ.0.)GO TO 50
   40 CONTINUE
   50 TYPE 70
      ACCEPT 60, YN
      IF(YN.NE.'Y')GO TO 80
      IREAD=2
      GO TO 10
   60 FORMAT(A1)
   70 FORMAT($' ANOTHER FIT? (Y/N) ')
   80 CALL EXIT
      END
C--------------------BLOCK--DATA-----------------------------------
C
      BLOCK DATA
      LOGICAL*1 MEMBER,TFLAG,CHKPRT,CHKOUT,FLAG,MMBR,TSTNAM
      COMMON /CARD1/ MEMBER(10),M,NDS,NITER,NSTEP,DAMPER,DESIRE,
     C    Q,POWER,KPOWER,TFLAG,CHKPRT,CHKOUT
      COMMON /CARD2/ CAYLO,CAYHI
      COMMON /CARD3/ FLAG(122),IFLGVR(13)
      COMMON /CARD4/ VAL(131)
      COMMON /CARD5/ NSHELL(10),NENV(10),MMBR(10,10),QV(10),POWERV(10)
      COMMON /CARD6/ TSTNAM(10)
      COMMON /CARD7/ ENVLOP(100,2,10)
      COMMON /CARD8/ NPS(10),XKAY(20,10),PSI(20,10)
      COMMON /LIMITS/ XMIN,XMAX,YMIN,YMAX
      COMMON /YSET/ YSMIN,YSMAX
      COMMON /TERMS/ TERM(140)
      COMMON /PLOT/ XDEL,NEX,NXTX,YDEL,NEY,NYTX,IMIN,IMAX,JMIN,JMAX
      DATA MEMBER /10*' '/
      DATA CHKPRT,CHKOUT,TFLAG/.FALSE.,.FALSE.,.FALSE./
      DATA TERM /140*0.0/
      DATA VAL /131*0.0/
      DATA IFLGVR/'N ','A ','R ','SI','B ','C ','P0','P1','P2',
     C    'P3','G ','H ','E0'/
      DATA DESIRE,M,DAMPER,KPOWER/1.0,0,0.0,0/
      DATA TSTNAM/10*' '/
      DATA MMBR/100*' '/
      DATA NPS/10*0/
      END
C--------------------SUBROUTINE--ADJUST------------------------------
C
      SUBROUTINE ADJUST(SOLN,NBEST,CHI,CAY,Y,NPTS,EZERO,SIGMIN)
C * SUBROUTINE TO ADJUST THE PARAMETER VALUES AFTER EACH ITERATION
C * WHEN THE EDGE ENERGY IS A VARIABLE THE RAW DATA IS ADJUSTED
      COMMON /CARD3/ FLAG(122),IFLGVR(13)
      COMMON /CARD4/ VALS(131)
      DIMENSION CHI(1024),CAY(1024),Y(1024),SOLN(121,2)
      REAL*8 CHI
      LOGICAL*1 FLAG
      EZ=VALS(131)
      DO 10 J=11,131
   10 VALS(J)=SOLN(J-10,NBEST)
      DO 20 J=1,10
      IF(VALS(J+40).EQ.SIGMIN) WRITE(6,50) J,SIGMIN
   20 VALS(J)=VALS(J+10)*VALS(J+20)
      IF(.NOT.FLAG(121))GO TO 40
C * THE EDGE ENERGY IS VARYING; THE K ARRAY IS MODIFIED AT 
C * EACH ITERATION
      DELTA=SOLN(121,NBEST)-EZ
      EZERO=SOLN(121,NBEST)
      VSQDEL=.262467*DELTA
      DO 30 L=1,NPTS
      CAYSQ=CAY(L)*CAY(L)
   30 CAY(L)=SQRT(CAYSQ-VSQDEL)
   40 CONTINUE
      RETURN 
   50 FORMAT('0***  SIGMA(',I1,') IS SET TO ',1PE10.3,'  ***'/)
      END
C----------------SUBROUTINE--BSTPLT----------------------------------
C
      SUBROUTINE BSTPLT(ITERNO,VPREV,EZERO,Y,NPTS,CHI,CAY,DAMBST)
C * SUBROUTINE TO PLOT THE FIT FOR EACH ITERATION
      DIMENSION Y(1024),CAY(1024),CHI(1024)
      COMMON /CARD1/ MEMBER(10),M,NDS,NITER,NSTEP,DAMPER,DESIRE,Q,
     C    POWER,KPOWER,TFLAG,CHKPRT,CHKOUT
      COMMON /LIMITS/ XMIN,XMAX,YMIN,YMAX,SIGMIN
      COMMON /PLOT/ XDEL,NEX,NXTX,YDEL,NEY,NYTX,IMIN,IMAX,JMIN,JMAX
      COMMON /YSET/ YSMIN,YSMAX
      LOGICAL*1 MEMBER,TFLAG,CHKPRT,CHKOUT,YN
      REAL*8 CHI
      J=0
      STDV=0.0
      IF(VPREV.NE.0.0) STDV=SQRT(VPREV)
      EZN=EZERO
      XMIN=CAY(1)
      XMAX=CAY(NPTS)+1.
      YMIN=YSMIN
      YMAX=YSMAX
	CALL VECMOD
      CALL NEWPAG
      IMIN=200
      IMAX=1020
      JMIN=150
      JMAX=770
      CALL SCPLT(1)
      CALL AXES(XMIN,XMAX,XDEL,NEX,NXTX,XT)
      CALL AXES(YMIN,YMAX,YDEL,NEY,NYTX,YT)
      CALL SCPLT(2)
      CALL SCPLT(4)
      CALL TICKS(XT,XV,YT,YV)
	CALL TITLE
      CALL XLABEL(XT,XV)
      CALL YLABEL(YT,YV,KPOWER)
      CALL GRAPH(CHI,CAY,Y,NPTS,J)
      CALL MOVABS(LINWDT(1),LINHGT(2))
      CALL ALFMOD
      TYPE 30, ITERNO,STDV,EZERO,DAMBST
   30 FORMAT(3X,'ITERNO:',I2,2X,'STDEV: ',1PE10.3,2X,'EZERO: ',
     C    1PE11.3,2X,'DMPFAC: ',1PE10.3)
	TYPE *,' '
	TYPE *,' '
	TYPE *,' '
      ACCEPT 10, YN
   10 FORMAT (1A1)
      IF(YN.NE.'Y')GO TO 20
      CALL SCPLT(3)	!FOR HARDCOPY
      CALL IOWAIT(100)
20      RETURN
      END
C------------------SUBROUTINE--CARD1-----------------------------------
C
      SUBROUTINE CARD1
C * SUBROUTINE TO SETUP FIT PARAMETERS
      LOGICAL*1 MEMBER,MEMB(10),TF,TFL,TFLAG,CP,CKP,CHKPRT,CO,CKO,
     C    CHKOUT
      COMMON /CARD1/ MEMBER(10),M,NDS,NITER,NSTEP,DAMPER,DESIRE,Q,
     C    POWER,KPOWER,TFLAG,CHKPRT,CHKOUT
      DO 10 I=1,10
   10 IF(MEMBER(I).EQ.' ' .OR. MEMBER(I).EQ."0)MEMBER(I)=' '
      TYPE 60, MEMBER
      ACCEPT 70, MEMB
      IF(MEMB(1).EQ.' '.OR.MEMB(1).EQ.'/') GO TO 20
      DECODE(10,70,MEMB) MEMBER
   20 TYPE 80, M
      ACCEPT *, M
      IF(M.EQ.0)M=1
      TYPE 90, NDS
      ACCEPT *, NDS
      TYPE 100, NITER
      ACCEPT *, NITER
      TYPE 110, NSTEP
      ACCEPT *, NSTEP
      TYPE 120, DAMPER
      ACCEPT *, DAMPER
      TYPE 130, DESIRE
      ACCEPT *, DESIRE
      TYPE 140, Q
      ACCEPT *, Q
      IF(Q.EQ.0.) Q=1.0
      TYPE 150, POWER
      ACCEPT *, POWER
      TYPE 160, KPOWER
      ACCEPT *, KPOWER
      TF='N'
      IF(TFLAG) TF='Y'
      TYPE 170, TF
      ACCEPT 180, TFL
      IF(TFL.EQ.' '.OR.TFL.EQ.'/') GO TO 30
      TFLAG=.FALSE.
      IF(TFL.EQ.'Y')TFLAG=.TRUE.
   30 CP='N'
      IF(CHKPRT) CP='Y'
      TYPE 190, CP
      ACCEPT 180, CKP
      IF(CKP.EQ.' '.OR.CKP.EQ.'/') GO TO 40
      CHKPRT=.FALSE.
      IF(CKP.EQ.'Y') CHKPRT=.TRUE.
   40 CO='N'
      IF(CHKOUT) CO='Y'
      TYPE 200, CO
      ACCEPT 180, CKO
      IF(CKO.EQ.' '.OR.CKO.EQ.'/') GO TO 50
      CHKOUT=.FALSE.
      IF(CKO.EQ.'Y') CHKOUT=.TRUE.
   50 GO TO 210
   60 FORMAT('$CASENAME: ',10A1,' TYPE CHANGE: ')
   70 FORMAT(10A1)
   80 FORMAT('$NO. OF SHELLS: ',I2,' TYPE CHANGE: ')
   90 FORMAT('$NO. OF DAMPS PER ITERATION: ',I2,' TYPE CHANGE: ')
  100 FORMAT('$NO. OF ITERATIONS: ',I2,' TYPE CHANGE: ')
  110 FORMAT('$FIT INCREMENT: ',I2,' TYPE CHANGE: ')
  120 FORMAT('$DAMPING FACTOR: ',F10.3,' TYPE CHANGE: ')
  130 FORMAT('$PERCENT CHANGE STOP: ',F10.3,' TYPE CHANGE: ')
  140 FORMAT(' MULTIPLY DATA BY  Q*K**N ',/,
     C    '$Q :',F10.3,' TYPE CHANGE: ')
  150 FORMAT('$N: ',F10.3,' TYPE CHANGE: ')
  160 FORMAT(' WEIGHTING  CHI=K**M*SUM[ ] ',/,
     C    '$M: ',I2,' TYPE CHANGE: ')
  170 FORMAT('$TESTCASE? ',A1,' TYPE CHANGE: ')
  180 FORMAT(A1)
  190 FORMAT('$SETUP MAP? ',A1,' TYPE CHANGE: ')
  200 FORMAT('$STORE SETUP? ',A1,' TYPE CHANGE: ')
  210 RETURN 
      END
C------------------SUBROUTINE--CARD2-----------------------------------
C 
      SUBROUTINE CARD2(XKMIN,XKMAX)
C * SUBROUTINE TO SETUP THE K-RANGE TO FIT
      COMMON /CARD2/ CAYLO,CAYHI
      TYPE 10, XKMIN,XKMAX
      TYPE 20, CAYLO
      ACCEPT *, CAYLO
      TYPE 30, CAYHI
      ACCEPT *, CAYHI
      GO TO 40
   10 FORMAT(/,' MAXIMUM K-RANGE AVAILABLE ',/,
     C    ' XKMIN: ',F10.3,/,' XKMAX: ',F10.3)
   20 FORMAT(' K-RANGE TO FIT ',/,'$KMIN: ',F10.3,' TYPE CHANGE: ')
   30 FORMAT('$KMAX: ',F10.3,' TYPE CHANGE: ')
   40 RETURN 
      END
C--------------------SUBROUTINE--CARD3-----------------------------------
C
      SUBROUTINE CARD3
C * SUBROUTINE TO DETERMINE THE PARAMETERS TO VARY 
      LOGICAL*1 FLAG,FL,FLG,MEMBER
      COMMON /CARD1/ MEMBER(10),M
      COMMON /CARD3/ FLAG(122),IFLGVR(13)
      J=1
10      TYPE 40, J
      DO 20 I=1,12
      K=(I-1)*10+J
      FL='N'
      IF(FLAG(K)) FL='Y'
      TYPE 50, IFLGVR(I),FL
      ACCEPT 60, FLG
      IF(FLG.EQ.' '.OR.FLG.EQ.'/')GO TO 20
      FLAG(K)=.FALSE.
      IF(FLG.EQ.'Y') FLAG(K)=.TRUE.
   20 CONTINUE
      IF(J.GT.1) GO TO 30
      FL='N'
      IF(FLAG(121)) FL='Y'
      TYPE 50, IFLGVR(13),FL
      ACCEPT 60,FLG
      IF(FLG.EQ.' '.OR.FLG.EQ.'/')GO TO 30
      FLAG(121)=.FALSE.
      IF(FLG.EQ.'Y') FLAG(121)=.TRUE.
   30 IF(M.EQ.J) GO TO 70
      J=J+1
      GO TO 10
   40 FORMAT(/,' PARAMETERS TO VARY IN THE FIT ',/,' SHELL: ',I2)
   50 FORMAT('$',A2,'? ',A1,'   TYPE CHANGE: ')
   60 FORMAT(A1)
   70 RETURN 
      END
C-----------------SUBROUTINE--CARD4----------------------------------
C
      SUBROUTINE CARD4
C * SUBROUTINE TO LOAD TRAIL PARAMETERS FOR FITS
      LOGICAL*1 FLAG,MEMBER
      COMMON /CARD1/ MEMBER(10),M
      COMMON /CARD3/ FLAG(122),IFLGVR(13)
      COMMON /CARD4/ ZA(10),Z(10),A(10),R(10),SIGMA(10),B(10),
     C    GAMMA(10),P(10,4),G(10),H(10),EZ
      J=1
10      TYPE 30, J
      TYPE 40, IFLGVR(1),Z(J)
      ACCEPT *, Z(J)
      TYPE 40, IFLGVR(2),A(J)
      ACCEPT *, A(J)
      TYPE 40, IFLGVR(3),R(J)
      ACCEPT *, R(J)
      TYPE 40, IFLGVR(4),SIGMA(J)
      ACCEPT *, SIGMA(J)
      TYPE 40, IFLGVR(5),B(J)
      ACCEPT *, B(J)
      TYPE 40, IFLGVR(6),GAMMA(J)
      ACCEPT *,GAMMA(J)
      TYPE 40, IFLGVR(7),P(J,1)
      ACCEPT *, P(J,1)
      TYPE 40, IFLGVR(8),P(J,2)
      ACCEPT *, P(J,2)
      TYPE 40, IFLGVR(9),P(J,3)
      ACCEPT *, P(J,3)
      TYPE 40, IFLGVR(10),P(J,4)
      ACCEPT *, P(J,4)
      TYPE 40, IFLGVR(11),G(J)
      ACCEPT *, G(J)
      TYPE 40, IFLGVR(12),H(J)
      ACCEPT *, H(J)
      IF(J.GT.1) GO TO 20
      TYPE 40, IFLGVR(13),EZ
      ACCEPT *, EZ
   20 IF(J.EQ.M) GO TO 50
      J=J+1
      GO TO 10
   30 FORMAT(/,' TRIAL PARAMETERS ',/,' SHELL: ',I2)
   40 FORMAT('$',A2,': ',1PE12.5,'   TYPE CHANGE: ')
   50 RETURN
      END
C----------------SUBROUTINE--CARD5--------------------------------------
C
      SUBROUTINE CARD5
C * SUBROUTINE TO SET UP A RECALL OF REFERENCE ENVELOPES FROM DISK FILE
      LOGICAL*1 MEMBER,MMBR,MEMBR(10)
      COMMON /CARD1/ MEMBER(10),M
      COMMON /CARD5/ NSHELL(10),NENV(10),MMBR(10,10),QV(10),POWERV(10)
      J=1
10      TYPE 70, J
      NSHELL(J)=J
   20 TYPE 80, NENV(J)
      ACCEPT *, NENV(J)
      TYPE 90, (MMBR(I,J),I=1,10)
      ACCEPT 100, MEMBR
      IF(MEMBR(1).EQ.'/') GO TO 60
      ITEST=0
      DO 30 I=1,10
      IF(MEMBR(I).EQ.' '.OR. MEMBR(I).EQ."0)MEMBR(I)=' '
   30 IF(MEMBR(I).NE.' ')ITEST=1
      IF(ITEST.EQ.0)GO TO 40
      DECODE(10,100,MEMBR) (MMBR(I,J),I=1,10)
      GO TO 60
   40 DO 50 I=1,10
   50 MMBR(I,J)=' '
   60 TYPE 110, J,J,J,QV(J)
      ACCEPT *, QV(J)
      TYPE 120, J,POWERV(J)
      ACCEPT *, POWERV(J)
      IF(J.EQ.M) GO TO 130
      J=J+1
      GO TO 10
   70 FORMAT(/,' REFERENCE ENVELOPE DATA ',/,' SHELL: ',I2)
   80 FORMAT('$NO. OF ENVELOPE POINTS [0 TO 100]: ',I3,' TYPE CHANGE: ')
   90 FORMAT('$CASE MEMBER: ',10A1,' TYPE CHANGE: ')
  100 FORMAT(10A1)
  110 FORMAT(' WEIGHTED BY Q',I1,'*K**N',I1,/,
     C    '$Q',I1,': ',F10.3,' TYPE CHANGE: ')
  120 FORMAT('$N',I1,': ',F10.3,' TYPE CHANGE: ')
  130 RETURN
      END
C-----------------SUBROUTINE--CARD6-------------------------------------
C
      SUBROUTINE CARD6
C * SUBROUTINE TO ENTER A TESTCASE NAME FOR THE GENERATION OF A TEST
C * FUNCTION AND ITS STORAGE ON DISK
      LOGICAL*1 TSTNAM,TESTNM(10)
      COMMON /CARD6/ TSTNAM(10)
      TYPE 20, (TSTNAM(I),I=1,10)
      ACCEPT 30, (TESTNM(I),I=1,10)
      IF(TESTNM(1).EQ.' '.OR.TESTNM(1).EQ.'/')GO TO 10
      DECODE(10,30,TESTNM) TSTNAM
   10 GO TO 40
   20 FORMAT(/,'$TESTCASE NAME: ',10A1,' TYPE CHANGE: ')
   30 FORMAT(10A1)
   40 RETURN
      END
C------------------SUBROUTINE--CARD7----------------------------------
C
      SUBROUTINE CARD7
C * SUBROUTINE TO PROMPT KEYBOARD INPUT OF ENVELOPE (UP TO 100 PAIRS OF
C * POINTS)
      LOGICAL*1 MEMBER,MMBR
      COMMON /CARD1/ MEMBER(10),M
      COMMON /CARD5/ NSHELL(10),NENV(10),MMBR(10,10),QV(10),POWERV(10)
      COMMON /CARD7/ ENVLOP(100,2,10)
      ITST=0
      J=1
10      TYPE 60, NSHELL(J)
      DO 20 I=1,10
   20 IF(MMBR(I,J).NE.' ')ITST=1
      IF(ITST.EQ.1) TYPE 80,NSHELL(J)
      IF(ITST.EQ.1) GO TO 50
      DO 40 I=1,NENV(J)
      K=1
      IF(I.NE.30*K)GO TO 30
      K=K+1
   30 TYPE 70, I,ENVLOP(I,1,J),ENVLOP(I,2,J)
      ACCEPT *, ENVLOP(I,1,J),ENVLOP(I,2,J)
   40 CONTINUE
   50 IF(J.EQ.M) GO TO 90
      J=J+1
      GO TO 10
   60 FORMAT(/,' INPUT OF ENVELOPE PAIRS ',/,' SHELL: ',I2)
   70 FORMAT(' ENVELOPE PAIR',I3,3X,E10.4,',',E10.4,/,
     C    '$TYPE CHANGES: ')
   80 FORMAT(' SHELL',I2,' ENVELOPE IS SPECIFIED AS A FILE MEMBER IN',
     C    ' CARD5',/,' CHANGE CARD5 IF A TABULAR ENVELOPE IS DESIRED.')
   90 RETURN
      END
C------------------SUBROUTINE--CARD8-------------------------------------
C
      SUBROUTINE CARD8
C * PROGRAM TO PROMPT FOR KEYBOARD INPUT OF PHASE SHIFT (NOT > 20 POINTS)
      LOGICAL*1 MEMBER
      COMMON /CARD1/ MEMBER(10),M
      COMMON /CARD8/ NPS(10),XKAY(20,10),PSI(20,10)
      I=1
10      TYPE 50, I,NPS(I)
      ACCEPT *, NPS(I)
      IF(NPS(I).NE.0) GO TO 30
      DO 20 K=1,20
      PSI(K,I)=0.0
      XKAY(K,I)=0.0
   20 CONTINUE
      GO TO 70
   30 DO 40 KK=1,NPS(I)
      TYPE 60, KK,XKAY(KK,I),PSI(KK,I)
      ACCEPT *, XKAY(KK,I),PSI(KK,I)
   40 CONTINUE
      IF(I.EQ.M) GO TO 70
      I=I+1
      GO TO 10
   50 FORMAT(/,' INPUT PHASE SHIFT ',/,'$SHELL',I2,' HAS ',I2,
     C    ' PAIRS OF POINTS   TYPE CHANGE: ')
   60 FORMAT(' PHASE SHIFT PAIR:',I2,3X,E13.5,',',E13.5,/,
     C    '$TYPE CHANGES: ')
   70 RETURN
      END
C------------------SUBROUTINE--CASEFL-----------------------------------
C
      SUBROUTINE CASEFL(V,W,NPTS,MEMBER,XKMIN,XKMAX,IERR)
C * THIS PROGRAM READS IN A CASEFILE FROM DISK AND PASSES THE DATA
C * ON TO THE FITTING SETUP ROUTINE
      DIMENSION V(1024),W(1024)
      LOGICAL*1 KOPT,LOPT,NAMES(10),BLANK,ZERO,MEMBER(10)
      DATA IOLD/'O '/,BLANK/' '/,ZERO/O0/
      ISML=0
      IFT=0
      SGN=-1.
      DO 10 I=1,10
   10 IF(MEMBER(I).NE. BLANK .AND. MEMBER(I) .NE. ZERO)GO TO 20
      GO TO 30
   20 CALL FILES(3,MEMBER,10,0,1,0,LUNIT)
      READ(LUNIT,120,ERR=30) (NAMES(I),I=1,10)
      GO TO 40
   30 TYPE 110
      IERR=1
      GO TO 200
   40 READ(LUNIT,130,ERR=30) KOPT,LOPT,FRM,ISUM,IDGL
      IF(ISUM.LE.3) GO TO 60
      SGN=1.
      READ(LUNIT,140) FLINE,DSPACE,SNOFF,EEDGE,ESHFT,NX
      IF(NX.GE.1)GO TO 50
      IF(ISUM.EQ.4)GO TO 60
   50 ISML=1
      READ(LUNIT,150) XI,XJ
      READ(LUNIT,160) XNMIN,YNMIN,STPHGT,NEL1,NEL2,MCM,NR
      IF(ISUM.EQ.6) GO TO 60
      IFT=1
      READ(LUNIT,170) XKMIN,XKMAX,KNUM,SCPTS,IPRCNT,KPOWR
      IF(ISUM.LT.9) GO TO 60
      READ(LUNIT,180) XRMIN,XRMAX
   60 READ(LUNIT,190,END=70) (V(I),W(I),I=1,1024)
   70 NPTS=I-1
   80 IF(V(NPTS).NE.0.)GO TO 90
      NPTS=NPTS-1
      XKMIN=V(1)
      XKMAX=V(NPTS)
      GO TO 80
   90 DO 100 I=1,NPTS
  100 V(I)=SGN*V(I)
C      FORM=0.
C      IF(FORM.EQ.IOLD)FORM=1.
C      KLOG=0
C      ICAPX=1
C      IF(ISUM.EQ.3)KLOG=1
C      IF(ISUM.EQ.4)ICAPX=2
C      IF(ISUM.EQ.5)ICAPX=5
C      IF(ISUM.EQ.6)ICAPX=3
C      IF(ISUM.EQ.8)ICAPX=4
C      IF(ISUM.EQ.9)ICAPX=3
C      IF(ISUM.EQ.9)KLOG=3
C      ICAPY=ISUM-ICAPY-KLOG
      CALL FILES(3,MEMBER,10,0,0,0,LUNIT)
  110 FORMAT(' FILE NOT FOUND OR ERROR READING FILE')
  120 FORMAT(/,21X,10A1)
  130 FORMAT(40X,/,17X,/,10X,A1,1X,A1,2X,A2,2(2X,I1))
  140 FORMAT(29X,F10.0,/,26X,F10.5,/,21X,F10.2,/,30X,F10.2,
     C    /,18X,F10.2,22X,I2)
  150 FORMAT(37X,2(2PE12.4,5X))
  160 FORMAT(2(8X,2PE12.4),10X,1PE12.4,2(2X,I3),24X,I2,11X,I2)
  170 FORMAT(7X,1PE12.4,9X,2PE12.4,/,22X,I4,24X,0PF6.0,
     C    /,20X,I2,19X,I1)
  180 FORMAT(25X,/,7X,1PE12.4,9X,1PE12.4,10X)
  190 FORMAT(4(5X,1P2E13.5))
  200 RETURN
      END
C----------------------SUBROUTINE--CHIGEN-----------------------------
C
      SUBROUTINE CHIGEN(N,Z,A,B,C,R,SIGMA,P,G,H,M,
     C    KPOWER,CHI,CAY,Y,CAYO,NPTS)
C * THIS SUBROUTINE GENERATES ALL OF THE CHI(L)'S AND THE OTHER FUNTIONS
C * F(L),T(L),PSI(L) AND LAMBDA(L) FOR SUBSEQUENT USE IN EVALUATING THE
C * PARTIAL DERIVATIVES
      DIMENSION CHI(1024),CAY(1024),Y(1024),CAYO(1024)
      DIMENSION Z(10),A(10),B(10),C(10),R(10),SIGMA(10),
     C    P(10,4),G(10),H(10)
      COMMON /CARD5/ NSHELL(10),NENV(10),MMBR(10,10),QV(10),POWERV(10)
      COMMON /CARD7/ ENVLOP(100,2,10)
      COMMON /CARD8/ NPS(10),XKAY(20,10),PSI(20,10)
      COMMON /TERMS/  DF(10),DPSI(10),DLAMDA(10),DT(10),DCHIL(10),
     C    DPWR(10),DRK(10),DCOST(10),DSINT(10),DNOSIN(10),DCAYMC(10),
     C    DELAST(10),DC(10),D2RONL(10)
      REAL*8 CHI,CHIN
      LOGICAL*1 MMBR
      CHIN=0.0D0
      CAYN=CAY(N)
      CAYON=CAYO(N)
      CAYN2=CAYN*CAYN
      CAYN3=CAYN*CAYN2
      CAYNKP=1.
      IF(KPOWER.NE.0)CAYNKP=CAYN**KPOWER
      DO 10 L=1,M
      DC(L)=2.0*SIGMA(L)*CAYN2
      DPWR(L)=0.
      IF(DC(L).LT.100.)DPWR(L)=EXP(-DC(L))
      DCAYMC(L)=CAYN-C(L)
      IF(A(L).NE.1.0.OR.B(L).NE.1.0.OR.C(L).NE.1.0)
     C    DF(L)=A(L)/(1.0+B(L)*B(L)*DCAYMC(L)*DCAYMC(L))
      ENVSAV=ENVLOP(1,1,L)
      NSAV=1
C * IN THE FOLLOWING STATEMENT THE "ORIGINAL" VALUE OF K IS BEING
C * USED TO LOOK UP THE ENVELOPE VALUE RATHER THAN THE CURRENT K,
C * WHICH MAY HAVE BEEN MODIFIED IF THE EDGE ENERGY IS VARYING.
      IF(A(L).EQ.1.0.AND.B(L).EQ.1.0.AND.C(L).EQ.1.0)
     C    DF(L)=TBLP(ENVLOP(1,1,L),ENVLOP(1,2,L),CAYON,ENVSAV,
     C    NENV(L),NSAV,0)
      XKSAV=XKAY(1,L)
      MSAV=1
      IF(P(L,1).NE.1..OR.P(L,2).NE.1..OR.P(L,3).NE.1..OR.P(L,4).NE.1.)
     C    DPSI(L)=P(L,1)+CAYN*(P(L,2)+CAYN*P(L,3))+P(L,4)/CAYN3
      IF(P(L,1).EQ.1.AND.P(L,2).EQ.1.AND.P(L,3).EQ.1.AND.P(L,4).EQ.1.)
     C    DPSI(L)=TBLP(XKAY(1,L),PSI(1,L),CAYN,XKSAV,NPS(L),MSAV,0)
      DLAMDA(L)=G(L)*CAYN+H(L)
      DTL=2.0*CAYN*R(L)+DPSI(L)
      DT(L)=DTL
      DSINT(L)=SIN(DTL)
      DCOST(L)=COS(DTL)
      D2RONL(L)=2.0*R(L)/DLAMDA(L)
      DELAST(L)=EXP(-D2RONL(L))
      DNOSIN(L)=Z(L)*DPWR(L)*DF(L)/R(L)/R(L)*DELAST(L)*CAYNKP
      DCHIL(L)=Z(L)*DPWR(L)*DF(L)/R(L)/R(L)*DSINT(L)*DELAST(L)*CAYNKP
      CHIN=CHIN+DCHIL(L)
   10 CONTINUE
      CHI(N)=CHIN
      RETURN
      END
C----------------------SUBROUTINE--COEGEN------------------------------
C
      SUBROUTINE COEGEN(CHI,CAY,Y,CAYO,NPTS,IP,ISH,
     C    ITERNO,COEFF,RHS,MP)
C * SUBROUTINE TO GENERATE THE COEFFICIENT MATRIX WITH ZERO DAMPING
      LOGICAL*1 CHKOUT,CHKPRT,TFLAG,TSTNAM,MEMBER
      DIMENSION CHI(1024),CAY(1024),Y(1024),CAYO(1024),COEFF(MP,MP),
     C    RHS(111),IP(111),ISH(111),PARSHL(111)
      REAL*8 CHI,PARSHL,RHS
      COMMON /CARD1/ MEMBER(10),M,NDS,NITER,NSTEP,DAMPER,DESIRE,Q,
     C    POWER,KPOWER,TFLAG,CHKPRT,CHKOUT
      COMMON /CARD4/ ZA(10),Z(10),A(10),R(10),SIGMA(10),B(10),
     C    C(10),P(10,4),G(10),H(10),EZ
      COMMON /CARD6/ TSTNAM(10)
C * CARD TO CAUSE GENERATION OF ANALYTICAL DATA FOR TEST PURPOSES
C ***************************************************************
      REAL*4 F2(3),F3(3),F4(3),F5(3)
      DATA F2/'(F7.','2,1X',')   '/
      DATA F3/'(F7.','3,1X',')   '/
      DATA F4/'(F7.','4,1X',')   '/
      DATA F5/'(F7.','5,1X',')   '/
C **************************************************************
      DO 20 L=1,MP
      DO 10 I=L,MP
   10 COEFF(I,L)=0.0
      RHS(L)=0.0D0
   20 CONTINUE
      DO 50 N=1,NPTS
      NAN=N
      CALL CHIGEN(NAN,Z,A,B,C,R,SIGMA,P,G,H,M,KPOWER,CHI,
     C    CAY,Y,CAYO,NPTS)
      CALL DCHIDP(NAN,CHI,CAY,Y,CAYO,NPTS,IP,ISH,MP,PARSHL)
      YDEL=Y(N)-CHI(N)
      DO 40 L=1,MP
      DO 30 I=L,MP
      COEFF(I,L)=COEFF(I,L)+PARSHL(I)*PARSHL(L)
   30 CONTINUE
      RHS(L)=RHS(L)+YDEL*PARSHL(L)
   40 CONTINUE
   50 CONTINUE
C * FILL IN THE UPPER TRIANGLE OF THE MATRIX
      DO 60 J=2,MP
      L=J-1
      DO 60 I=1,L
   60 COEFF(I,J)=COEFF(J,I)
C * CARD TO CAUSE GENERATION OF ANALYTICAL DATA FOR TEST PURPOSES
C ***************************************************************
      IF(.NOT.TFLAG) GO TO 100
      MM=M
      IF(MM.GT.4)MM=4
C *      ENCODE(EZ,1,LABEL(7),8,F4,IERRTN)
      DO 70 J=1,MM
C *      ENCODE( Z(J)      ,1,LABEL(15 +96*(J-1)),8,F2,IERRTN)
C *      ENCODE( A(J)      ,1,LABEL(23 +96*(J-1)),8,F2,IERRTN)
C *      ENCODE( R(J)      ,1,LABEL(31 +96*(J-1)),8,F3,IERRTN)
C *      ENCODE( SIGMA(J)  ,1,LABEL(39 +96*(J-1)),8,F2,IERRTN)
C *      ENCODE( B(J)      ,1,LABEL(47 +96*(J-1)),8,F3,IERRTN)
C *      ENCODE( C(J)      ,1,LABEL(55 +96*(J-1)),8,F4,IERRTN)
C *      ENCODE( P(J,1)    ,1,LABEL(63 +96*(J-1)),8,F5,IERRTN)
C *      ENCODE( P(J,2)    ,1,LABEL(71 +96*(J-1)),8,F4,IERRTN)
C *      ENCODE( P(J,3)    ,1,LABEL(79 +96*(J-1)),8,F4,IERRTN)
C *      ENCODE( P(J,4)    ,1,LABEL(87 +96*(J-1)),8,F3,IERRTN)
C *      ENCODE( G(J)      ,1,LABEL(95 +96*(J-1)),8,F2,IERRTN)
C *      ENCODE( H(J)      ,1,LABEL(103+96*(J-1)),8,F4,IERRTN)
C ***************************************************************
   70 CONTINUE
C *      WRITE(5,6001) (TSTNAM(J),J=1,10),(LABEL(K),K=1,120)
      NPTS2=2*NPTS
   80 CONTINUE
   90 FORMAT('0TEST DATA GENERATED AS FOLLOWS:'/,
     C    1X,10A1,120A1)
  100 RETURN
      END
C-------------------SUBROUTINE--DCHIDP-----------------------------------
C
      SUBROUTINE DCHIDP(N,CHI,CAY,Y,CAYO,NPTS,IP,ISH,MP,PARSHL)
C * SUBROUTINE TO GENERATE THE COMPLETE SET OF PARTIAL DERIVATIVES FOR
C * THE N'TH DATA POINT
      COMMON /CARD1/ MEMBER(10),M,NDS,NITER,NSTEP,DAMPER,DESIRE,Q,
     C    POWER,KPOWER,TFLAG,CHKPRT,CHKOUT
      COMMON /CARD4/ ZA(10),Z(10),A(10),R(10),SIGMA(10),B(10),
     C    GAMMA(10),P(10,4),G(10),H(10),EZ
      COMMON /TERMS/ DF(10),DPSI(10),DLAMDA(10),DT(10),DCHIL(10),
     C    DPWR(10),DRK(10),DCOST(10),DSINT(10),DNOSIN(10),DCAYMC(10),
     C    DELAST(10),DC(10),D2RONL(10)
      DIMENSION CHI(1024),CAY(1024),CAYO(1024),Y(1024),IP(111),ISH(111),
     C    PARSHL(111)
      REAL*8 CHI,PARSHL,VALUE
      LOGICAL*1 MEMBER,TFLAG,CHKPRT,CHKOUT
      DO 110 L=1,MP
      IL=IP(L)
      IF(IL.EQ.0)GO TO 120
      ISHL=ISH(L)
      GO TO (10,20,30,40,50,60,60,60,60,70,80,90),IL
C * THE PARAMETER IS AN N-A PRODUCT
   10 CONTINUE
      PARSHL(L)=DCHIL(ISHL)/ZA(ISHL)
      GO TO 110
C * THE PARAMETER IS AN R
   20 CALL DDR(N,ISHL,VALUE,CHI,CAY,Y)
      GO TO 100
C * THE PARAMETER IS A SIGMA
   30 CALL DDSIG(N,ISHL,VALUE,CHI,CAY,Y)
      GO TO 100
C * THE PARAMETER IS A B
   40 CALL DDB(N,ISHL,VALUE)
      GO TO 100
C * THE PARAMETER IS A GAMMA
   50 CALL DDGAM(N,ISHL,VALUE)
      GO TO 100
C * THE PARAMETER IS A PSI COEFFICIENT
   60 CALL DDPSI(N,ISHL,IL-5,VALUE,CHI,CAY,Y)
      GO TO 100
C * THE PARAMETER IS A LAMBDA SLOPE
   70 CALL DDLSL(N,ISHL,VALUE,CHI,CAY,Y)
      GO TO 100
C * THE PARAMETER IS A LAMBDA CONSTANT
   80 CALL DDLC(N,ISHL,VALUE)
      GO TO 100
C * THE PARAMETER IS THE EDGE ENERGY
   90 CALL DDEDGE(N,VALUE,CHI,CAY,Y,CAYO)
  100 CONTINUE
      PARSHL(L)=VALUE
  110 CONTINUE
  120 CONTINUE
      RETURN
      END
C--------------------SUBROUTINE--DDB-----------------------------------
C
      SUBROUTINE DDB(N,L,VALUE)
C * SUBROUTINE TO EVALUATE THE PARTIAL DERIVATIVE OF CHI WITH RESPECT TO
C * THE PARAMETER B FOR SHELL L, USING CAY(N)
      COMMON /CARD4/ ZA(10),Z(10),A(10),R(10),SIGMA(10),B(10),
     C    GAMMA(10),P(10,4),G(10),H(10),EZ
      COMMON /TERMS/ DF(10),DPSI(10),DLAMDA(10),DT(10),DCHIL(10),
     C    DPWR(10),DRK(10),DCOST(10),DSINT(10),DNOSIN(10),DCAYMC(10),
     C    DELAST(10),DC(10),D2RONL(10)
      REAL*8 VALUE
   10 VALUE=-DCHIL(L)*DF(L)/A(L)*DCAYMC(L)*2.0*B(L)*DCAYMC(L)
      RETURN
      END
C-------------------SUBROUTINE--DDEDGE---------------------------------
C
      SUBROUTINE DDEDGE(N,VALUE,CHI,CAY,Y,CAYO)
C * SUBROUTINE TO EVALUATE THE PARTIAL DERIVATIVE OF CHI WITH RESPECT 
C * TO THE EDGE ENERGY, EVALUATED AT CAY(N)
      COMMON /CARD1/ MEMBER(10),M,NDS,NITER,NSTEP,DAMPER,DESIRE,
     C    Q,POWER,KPOWER,TFLAG,CHKPRT,CHKOUT
      COMMON /CARD4/ ZA(10),Z(10),A(10),R(10),SIGMA(10),B(10),
     C    GAMMA(10),P(10,4),G(10),H(10),EZ
      COMMON /CARD5/ NSHELL(10),NENV(10),MMBR(10,10),QV(10),POWERV(10)
      COMMON /CARD7/ ENVLOP(100,2,10)
      COMMON /CARD8/ NPS(10),XKAY(20,10),PSI(20,10)
      COMMON /TERMS/ DF(10),DPSI(10),DLAMDA(10),DT(10),DCHIL(10),
     C    DPWR(10),DRK(10),DCOST(10),DSINT(10),DNOSIN(10),DCAYMC(10),
     C    DELAST(10),DC(10),D2RONL(10)
      DIMENSION CHI(1024),CAY(1024),CAYO(1024),Y(1024)
      LOGICAL*1 MEMBER,TFLAG,CHKPRT,CHKOUT,MMBR
      REAL*8 CHI,VALUE
      CAYN=CAY(N)
      CAYON=CAYO(N)
      VALUE=0.0D0
      DO 10 L=1,M
      DU=-2.0*DC(L)/CAYN
      DV=D2RONL(L)/DLAMDA(L)*G(L)
      DW=-2.0*DF(L)*B(L)*B(L)*DCAYMC(L)/A(L)
      IF(A(L).EQ.1..AND.B(L).EQ.1..AND.GAMMA(L).EQ.1.)
     C    DW=TBLSL(ENVLOP(1,1,L),ENVLOP(1,2,L),NENV(L),CAYON)/DF(L)
      DX=DCOST(L)
      DY=P(L,2)+2.0*(R(L)+CAYN*P(L,3))-3.0*P(L,4)/(CAYN**4)
      IF(P(L,1).EQ.1.AND.P(L,2).EQ.1.AND.P(L,3).EQ.1.AND.P(L,4).EQ.1)
     C    DY=2.0*R(L)+TBLSL(XKAY(1,L),PSI(1,L),NPS(L),CAYN)
      CAYPWR=KPOWER/CAYN
      DZ=(DU+DV+DW+CAYPWR)*DSINT(L)
      VALUE=VALUE-DNOSIN(L)*(DZ+DX*DY)
   10 CONTINUE
      VALUE=.131234*VALUE/CAYN
      RETURN
      END
C---------------------SUBROUTINE--DDGAM-----------------------------------
C
      SUBROUTINE DDGAM(N,L,VALUE)
C * SUBROUTINE TO EVALUATE THE PARTIAL DERIVATIVE OF CHI WITH RESPECT TO
C * THE PARAMETER GAMMA
      COMMON /CARD4/ ZA(10),Z(10),A(10),R(10),SIGMA(10),B(10),
     C    GAMMA(10),P(10,4),G(10),H(10),EZ
      COMMON /TERMS/ DF(10),DPSI(10),DLAMDA(10),DT(10),DCHIL(10),
     C    DPWR(10),DRK(10),DCOST(10),DSINT(10),DNOSIN(10),DCAYMC(10),
     C    DELAST(10),DC(10),D2RONL(10)
      REAL*8 VALUE
   10 VALUE=2.0*DCHIL(L)*B(L)*B(L)*DF(L)/A(L)*DCAYMC(L)
      RETURN
      END
C------------------------SUBROUTINE--DDLC----------------------------------
C
      SUBROUTINE DDLC(N,L,VALUE)
C * SUBROUTINE TO EVALUATE THE PARTIAL DERIVATIVE OF CHI WITH RESPECT TO 
C * THE INTERCEPT LAMBDA, FOR SHELL L, USING CAY(N)
      COMMON /CARD4/ ZA(10),Z(10),A(10),R(10),SIGMA(10),B(10),
     C    GAMMA(10),P(10,4),G(10),H(10),EZ
      COMMON /TERMS/ DF(10),DPSI(10),DLAMDA(10),DT(10),DCHIL(10),
     C    DPWR(10),DRK(10),DCOST(10),DSINT(10),DNOSIN(10),DCAYMC(10),
     C    DELAST(10),DC(10),D2RONL(10)
      REAL*8 VALUE
   10 VALUE=2.0*DCHIL(L)*R(L)/DLAMDA(L)/DLAMDA(L)
      RETURN
      END
C------------------------SUBROUTINE--DDLSL---------------------------------
C
      SUBROUTINE DDLSL(N,L,VALUE,CHI,CAY,Y)
C * SUBROUTINE TO EVALUATE THE PARTIAL DERIVATIVE OF CHI WITH RESPECT TO 
C * THE SLOPE OF LAMBDA, FOR SHELL L, USING CAY(N)
      COMMON /CARD4/ ZA(10),Z(10),A(10),R(10),SIGMA(10),B(10),
     C    GAMMA(10),P(10,4),G(10),H(10),EZ
      COMMON /TERMS/ DF(10),DPSI(10),DLAMDA(10),DT(10),DCHIL(10),
     C    DPWR(10),DRK(10),DCOST(10),DSINT(10),DNOSIN(10),DCAYMC(10),
     C    DELAST(10),DC(10),D2RONL(10)
      REAL*8 VALUE
      DIMENSION CHI(1024),CAY(1024),Y(1024)
      REAL*8 CHI
   10 VALUE=CAY(N)*DCHIL(L)*R(L)*2.0/DLAMDA(L)/DLAMDA(L)
      RETURN 
      END
C---------------------SUBROUTINE--DDPSI-----------------------------------
C
      SUBROUTINE DDPSI(N,L,J,VALUE,CHI,CAY,Y)
C * SUBROUTINE TO EVALUATE THE PARTIAL DERIVATIVE OF CHI WITH RESPECT TO
C * THE J'TH COEFFICIENT OF THE PARAMETER PSI, FOR SHELL L, USING CAY(N)
      COMMON /TERMS/ DF(10),DPSI(10),DLAMDA(10),DT(10),DCHIL(10),
     C    DPWR(10),DRK(10),DCOST(10),DSINT(10),DNOSIN(10),DCAYMC(10),
     C    DELAST(10),DC(10),D2RONL(10)
      DIMENSION CHI(1024),CAY(1024),Y(1024)
      REAL*8 CHI,VALUE
      IF(J.EQ.1) VALUE= DNOSIN(L)*DCOST(L)
      IF(J.EQ.2) VALUE= DNOSIN(L)*DCOST(L)*CAY(N)
      IF(J.EQ.3) VALUE= DNOSIN(L)*DCOST(L)*CAY(N)*CAY(N)
      IF(J.EQ.4) VALUE= DNOSIN(L)*DCOST(L)/CAY(N)/CAY(N)/CAY(N)
      RETURN 
      END
C---------------------SUBROUTINE--DDR-----------------------------------
C
      SUBROUTINE DDR(N,L,VALUE,CHI,CAY,Y)
C * SUBROUTINE TO EVALUATE THE PARTIAL DERIVATIVE OF CHI WITH RESPECT TO
C * THE PARAMETER R, FOR SHELL L, USING CAY(N)
      COMMON /CARD4/ ZA(10),Z(10),A(10),R(10),SIGMA(10),B(10),
     C    GAMMA(10),P(10,4),G(10),H(10),EZ
      COMMON /TERMS/ DF(10),DPSI(10),DLAMDA(10),DT(10),DCHIL(10),
     C    DPWR(10),DRK(10),DCOST(10),DSINT(10),DNOSIN(10),DCAYMC(10),
     C    DELAST(10),DC(10),D2RONL(10)
      DIMENSION CHI(1024),CAY(1024),Y(1024)
      REAL*8 CHI,VALUE
      DERIV=DSINT(L)*(1.0/R(L)+1.0/DLAMDA(L))
   10 VALUE=2.0*DNOSIN(L)*(CAY(N)*DCOST(L)-DERIV)
      RETURN
      END
C---------------------SUBROUTINE--DDSIG-----------------------------------
C
      SUBROUTINE DDSIG(N,L,VALUE,CHI,CAY,Y)
C * SUBROUTINE TO EVALUATE THE PARTIAL DERIVATIVE OF CHI WITH RESPECT TO
C * TO THE PARAMETER SIGMA, FOR SHELL L, USING CAY(N)
      COMMON /CARD4/ ZA(10),Z(10),A(10),R(10),SIGMA(10),B(10),
     C    GAMMA(10),P(10,4),G(10),H(10),EZ
      COMMON /TERMS/ DF(10),DPSI(10),DLAMDA(10),DT(10),DCHIL(10),
     C    DPWR(10),DRK(10),DCOST(10),DSINT(10),DNOSIN(10),DCAYMC(10),
     C    DELAST(10),DC(10),D2RONL(10)
      DIMENSION CHI(1024),CAY(1024),Y(1024)
      REAL*8 CHI,VALUE
   10 VALUE=-2.0*CAY(N)*CAY(N)*DCHIL(L)
      RETURN
      END
C----------------SUBROUTINE--EZGEN----------------------------------
C
      SUBROUTINE EZGEN(CAY,Y,CAYO,NPTS,EZERO)
      DIMENSION CHI(1024),CAY(1024),Y(1024),CAYO(1024),SNUM(2),
     C    SUM(2)
      REAL*8 SUM,CHI,DEL
      LOGICAL*1 MEMBER,TFLAG,CHKPRT,CHKOUT
      COMMON /CARD1/ MEMBER(10),M,NDS,NITER,NSTEP,DAMPER,DESIRE,
     C    Q,POWER,KPOWER,TFLAG,CHKPRT,CHKOUT
      COMMON /CARD4/ ZA(10),Z(10),A(10),R(10),SIGMA(10),B(10),
     C    C(10),P(10,4),G(10),H(10),EZ
      DO 10 N=1,2
      NAN=N
   10 CALL CHIGEN(NAN,Z,A,B,C,R,SIGMA,P,G,H,M,KPOWER,
     C    CHI,CAY,Y,CAYO,NPTS)
C * COMPARE PEAKS OF CHI WITH Y, SELECT EZERO TO MINIMIZE DELTA
      SNUM(1)=0.
      SNUM(2)=0.
      SUM(1)=0.0D0
      SUM(2)=0.0D0
      K=0
      IEXP=0
      NM1=NPTS-1
      DO 70 I=2,NM1
      CALL CHIGEN(I+1,Z,A,B,C,R,SIGMA,P,G,H,M,KPOWER,
     C    CHI,CAY,Y,CAYO,NPTS)
      IF(CHI(I).LT.CHI(I-1).OR.CHI(I).LT.CHI(I+1))GO TO 20
      CC=CAY(I)
      YC=CHI(I)
      K=2
      IF(Y(I).LT.Y(I-1).OR.Y(I).LT.Y(I+1))GO TO 40
      GO TO 30
   20 IF(Y(I).LT.Y(I-1).OR.Y(I).LT.Y(I+1))GO TO 70
   30 CY=CAY(I)
      YY=Y(I)
      K=K+1
   40 IF(K.EQ.IEXP)GO TO 50
      IF(K.NE.3)GO TO 60
C * COINCIDENT PEAKS, ADD DEL TO BOTH SUMS
      SNUM(1)=SNUM(1)+1.
      SNUM(2)=SNUM(2)+1.
      GO TO 60
C * SEPARATE PEAKS, ADD DEL TO APPROPRIATE SUM
   50 DEL=CY*CY-CC*CC
      SUM(K)=SUM(K)+DEL
      SNUM(K)=SNUM(K)+1.
   60 CONTINUE
      IEXP=3-K
      K=0
   70 CONTINUE
      IF(SNUM(1)*SNUM(2).LE.0.)RETURN
      DAMBST=0.
      CALL BSTPLT(0,0.0,EZERO,Y,NPTS,CHI,CAY,DAMBST)
      EZERO=SUM(1)/SNUM(1)
      IF(ABS(SNGL(SUM(2)/SNUM(2))).LT.ABS(EZERO))EZERO=SUM(2)/SNUM(2)
      EZERO=EZERO/.262467
      VSQDEL=.262467*(EZERO-EZ)
      DO 80 I=1,NPTS
      CAYSQ=CAY(I)*CAY(I)
   80 CAY(I)=SQRT(CAYSQ-VSQDEL)
      EZ=EZERO
      RETURN
      END
C---------------SUBROUTINE--AXES-----------------------------
C
      SUBROUTINE AXES(XMIN,XMAX,XDEL,NE,NXTX,TKT)
C * THIS SUBROUTINE COMPUTES THE NUMBER AND SCALING OF THE TICK
C * MARKS ALONG THE X-AXIS.
      INTEGER*4 IXMIN,IXMAX,IDELX
      DELX=XMAX-XMIN
      TKT=0.
      NE=0
   10 IF(DELX.GT. 5.)GO TO 20
      DELX=DELX*10.
      NE=NE+1
      GO TO 10
   20 IF(DELX.LE. 50.)GO TO 30
      DELX=DELX/10.
      NE=NE-1
      GO TO 20
   30 XMIN=XMIN*(10.**NE)
      IXMIN=XMIN-1.
      XMAX=XMAX*(10.**NE)
      IXMAX=XMAX+1.
      IDELX=IXMAX-IXMIN
      IF(IDELX.GT.10)GO TO 40
      XDEL=1.*(10.**(-NE))
      GO TO 60
   40 IF(IDELX.GT.25)GO TO 50
      XDEL=2.*(10.**(-NE))
      GO TO 60
   50 XDEL=5.*(10.**(-NE))
   60 XMIN=FLOAT(IXMIN)*(10.**(-NE))
      XMAX=FLOAT(IXMAX)*(10.**(-NE))
      NXTX=(XMAX-XMIN)/XDEL-.5
      IF(XMIN.LT.0. .AND. XMAX.GT.0.)TKT=1.
      IF(XMAX.LT.0. .AND. XMIN.GT.0.)TKT=1.
      RETURN
      END
C------------------SUBROUTINE--BOXES---------------------------------
C
      SUBROUTINE BOXES(V,W,N)
C * SUBROUTINE TO DRAW BOXES
      CALL MOVEA(V,W)
      CALL MOVREL(N,N)
      CALL DRWREL(-N,0)
      CALL DRWREL(0,-N)
      CALL DRWREL(N,0)
      CALL DRWREL(0,N)
      RETURN
      END
C----------------------SUBROUTINE--DECOMP----------------------------------
C
      SUBROUTINE DECOMP(NN,UL,IPS)
C * SUBROUTINE TO DECOMPOSE THE MATRIX
      DIMENSION UL(NN,NN),SCALES(111),IPS(NN)
      N=NN
C * INITIALIZE IPS AND SCALES
      DO 50 I=1,N
      IPS(I)=I
      ROWNRM=0.0
      DO 20 J=1,N
      IF(ROWNRM-ABS(UL(I,J))) 10,20,20
   10 ROWNRM=ABS(UL(I,J))
   20 CONTINUE
      IF(ROWNRM)30,40,30
   30 SCALES(I)=1.0/ROWNRM
      GO TO 50
   40 CALL SING(1)
      SCALES(I)=0.0
   50 CONTINUE
C * GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
      NM1=N-1
      DO 130 K=1,NM1
      BIG=0.0
      DO 70 I=K,N
      IP=IPS(I)
      SIZE=ABS(UL(IP,K))*SCALES(IP)
      IF(SIZE-BIG)70,70,60
   60 BIG=SIZE
      IDXPIV=I
   70 CONTINUE
      IF(BIG)90,80,90
   80 CALL SING(2)
      GO TO 130
   90 IF(IDXPIV-K)100,110,100
  100 J=IPS(K)
      IPS(K)=IPS(IDXPIV)
      IPS(IDXPIV)=J
  110 KP=IPS(K)
      PIVOT=UL(KP,K)
      KP1=K+1
      DO 120 I=KP1,N
      IP=IPS(I)
      EM=-UL(IP,K)/PIVOT
      UL(IP,K)=-EM
      DO 120 J=KP1,N
      UL(IP,J)=UL(IP,J)+EM*UL(KP,J)
C * INNER LOOP
  120 CONTINUE
  130 CONTINUE
      KP=IPS(N)
      IF(UL(KP,N)) 150,140,150
  140 CALL SING(2)
  150 RETURN
      END
C-------------SUBROUTINE--DIAMND-------------------------------------
C
      SUBROUTINE DIAMND(V,W,N)
C * SUBROUTINE TO DRAW DIAMONDS FOR DATA POINTS
      CALL MOVEA(V,W)
      CALL MOVREL(0,N)
      CALL DRWREL(N,-N)
      CALL DRWREL(-N,-N)
      CALL DRWREL(-N,N)
      CALL DRWREL(N,N)
      RETURN
      END
C-----------------SUBROUTINE--FILES------------------------------------
C
      SUBROUTINE FILES(L,FLNAME,N,J,K,NN,M)
C * SUBROUTINE TO HANDLE ALL OF THE FILE MANAGEMENT FOR THE EXAPLT PROGRAM
      LOGICAL*1 FLNAME(N),FILNAM(14),DEV1(3),DEV3(3),DEV10(3),ZERO
      DATA DEV1/' ',' ',' '/,DEV3/' ',' ',' '/,
     C	   DEV10/'T','I',':'/,
     C     ZERO/O0/,FILNAM/14*O0/
C * SYSTEM DEPENDENT CONVENTIONS FOR THIS COMPUTER SYSTEM ARE:
C
C *        L=1 ******** LOGICAL UNIT #1 ; SYSTEM DISKETTE
C *        L=2            "      "   #2 ; DEFAULT DISKETTE
C *        L=3            "      "   #3 ; DEFAULT DISKETTE
C *        L=4            "      "   #4 ; DEFAULT DISKETTE
C *        L=5            "      "   #5 ; READ FROM CONSOL TERMINAL
C *        L=6            "      "   #6 ; LINE PRINTER
C *        L=7            "      "   #7 ; WRITE TO CONSOL TERMINAL
C *        L=10           "      "   #10; DEVICE FOR GRAPHIC MODE OUTPUT
C
      IF(K) 10, 20, 30
   10 M=L
      GO TO 60
   20 M=L
      CLOSE(M)
      GO TO 60
   30 M=L
      DO 40 I=1,N
      IF(FLNAME(I).EQ.' ')FLNAME(I)=ZERO
   40 FILNAM(I+3)=FLNAME(I)
      DO 50 I=1,3
      IF(L.EQ.1) FILNAM(I)=DEV1(I)
      IF(L.GT.1.AND.L.LT.5) FILNAM(I)=DEV3(I)
      IF(L.GE.5) FILNAM(I)=DEV10(I)
   50 CONTINUE
      NUM=N+3
      IF(J.EQ.0) CALL ASSIGN(M,FILNAM,NUM,'OLD')
      IF(J.EQ.1) CALL ASSIGN(M,FILNAM,NUM,'NEW')
      IF(J.EQ.2) CALL ASSIGN(M,DEV10)
	IF (L.EQ.6) M=5
C *      IF(NN.EQ.1) DEFINE FILE M(10,10,U,NREC)
C *      IF(NN.EQ.3) DEFINE FILE M(100,25,U,NREC)
C *      IF(NN.EQ.4) DEFINE FILE M(1000,25,U,NREC)
   60 RETURN 
      END
C-------------------SUBROUTINE--MATINV-------------------------------------
C
      SUBROUTINE MATINV(A,B,C,N)
C * SUBROUTINE TO INVERT THE COEFFICIENT MATRIX
      DIMENSION A(N,N),B(111),C(111),NAR(4),IPS(111)
      REAL*8 B,C
      IF(N.EQ.1) GO TO 10
      NAR(1)=N
      NAR(2)=1
      NAR(3)=0
      NAR(4)=0
      CALL DECOMP(N,A,IPS)
      CALL SOLVE(N,B,C,A,IPS)
      GO TO 20
   10 C(1)=B(1)/A(1,1)
   20 CONTINUE
      RETURN 
      END
C------------------SUBROUTINE--SETUP---------------------------------
C
      SUBROUTINE SETUP(CAY,Y,NPTS,IREAD,IERR)
C * THIS PROGRAM SETS UP ALL THE CARD PARAMETERS FOR THE FITTING PROGRAM
      DIMENSION CAY(1024),Y(1024)
      COMMON /CARD1/ MEMBER(10),M,NDS,NITER,NSTEP,DAMPER,DESIRE,Q,
     C    POWER,KPOWER,TFLAG,CHKPRT,CHKOUT
      COMMON /CARD2/ CAYLO,CAYHI
      COMMON /CARD3/ FLAG(122),IFLGVR(13)
      COMMON /CARD4/ ZA(10),Z(10),A(10),R(10),SIGMA(10),B(10),
     C    GAMMA(10),P(10,4),G(10),H(10),EZ
      COMMON /CARD5/ NSHELL(10),NENV(10),MMBR(10,10),QV(10),POWERV(10)
      COMMON /CARD6/ TSTNAM(10)
      COMMON /CARD7/ ENVLOP(100,2,10)
      COMMON /CARD8/ NPS(10),XKAY(20,10),PSI(20,10)
      DATA IFLGVR /'N ','A ','R ','SI','B ','C ','P0','P1',
     C    'P2','P3','G ','H ','E0'/
      IFIRST=.FALSE.
   10 CALL NEWPAG
	CALL TOUTPT(24)		!SWITCH TO ANSI MODE
   20 TYPE 120 
      ACCEPT 130,ICARD
      IF(ICARD.EQ.0.AND.IREAD.NE.0) GO TO 40
      IF(ICARD.EQ.0.AND.IFIRST) GO TO 140
      IF(ICARD.EQ.1) GO TO 30
      IF(ICARD.EQ.2) GO TO 50
      IF(ICARD.EQ.3) GO TO 60
      IF(ICARD.EQ.4) GO TO 70
      IF(ICARD.EQ.5) GO TO 80
      IF(ICARD.EQ.6) GO TO 90
      IF(ICARD.EQ.7) GO TO 100
      IF(ICARD.EQ.8) GO TO 110
      IF(ICARD.EQ.0) TYPE *,'NO CASE NAME IN CARD1 !'
      GO TO 20
   30 CALL CARD1
   40 CALL CASEFL(CAY,Y,NPTS,MEMBER,XKMIN,XKMAX,IERR)
      IF(ICARD.EQ.0)GO TO 140
      IF(IERR.EQ.1)GO TO 20
      IFIRST=.TRUE.
      IREAD=0
      IERR=0
      GO TO 20
   50 CALL CARD2(XKMIN,XKMAX)
      GO TO 20
   60 CALL CARD3
      GO TO 20
   70 CALL CARD4
      GO TO 20
   80 CALL CARD5
      GO TO 20
   90 CALL CARD6
      GO TO 20
  100 CALL CARD7
      GO TO 20
  110 CALL CARD8
      GO TO 20
  120 FORMAT(/,'$TYPE INPUT CARD NUMBER (0-8) : ')
  130 FORMAT(I1)
  140 IF(IERR.EQ.1)GO TO 150
  150 RETURN
      END
C--------------------SUBROUTINE--SING------------------------------
C
      SUBROUTINE SING(IWHY)
C * SUBROUTINE TO CALL OUT THE MATRIX INVERSION ERRORS
   10 FORMAT('0MATRIX WITH ZERO ROW IN DECOMPOSE. ')
   20 FORMAT('0SINGULAR MATRIX IN DECOMPOSE.  ZERO DIVIDE IN SOLVE.')
   30 FORMAT('0NO CONVERGENCE IN IMPRUV. MATRIX IS NEARLY SINGULAR.')
      GO TO (40,50,60),IWHY
   40 TYPE 10
      GO TO 70
   50 TYPE 20
      GO TO 70
   60 TYPE 30
   70 RETURN
      END
C----------------------SUBROUTINE--SOLVE-------------------------------
C
      SUBROUTINE SOLVE(NN,B,X,UL,IPS)
C * SUBROUTINE TO SOLVE THE LINEAR SYSTEM OF EQUATIONS
      DIMENSION UL(NN,NN),B(NN),X(NN),IPS(NN)
      REAL*8 B,X
      N=NN
      NP1=N+1
      IP=IPS(1)
      X(1)=B(IP)
      DO 20 I=2,N
      IP=IPS(I)
      IM1=I-1
      SUM=0.0
      DO 10 J=1,IM1
   10 SUM=SUM+UL(IP,J)*X(J)
   20 X(I)=B(IP)-SUM
      IP=IPS(N)
      X(N)=X(N)/UL(IP,N)
      DO 40 IBACK=2,N
      I=NP1-IBACK
C *        I GOES (N-1),...,1
      IP=IPS(I)
      IP1=I+1
      SUM=0.0
      DO 30 J=IP1,N
   30 SUM=SUM+UL(IP,J)*X(J)
   40 X(I)=(X(I)-SUM)/UL(IP,I)
      RETURN
      END
C---------------SUBROUTINE--TICKS------------------------------------
C
      SUBROUTINE TICKS(XT,XV,YT,YV)
C * THIS SUBROUTINE DRAWS ALL OF THE TICK-MARKS ON EACH 
C * OF THE FOUR AXES OF THE FRAME SURROUNDING THE PLOT.
C * THIS EDITION HAS A PROVISION FOR A ZERO TICK-MARK SEARCH.
      COMMON /PLOT/ XDEL,NEX,NXTX,YDEL,NEY,NYTX,IMIN,IMAX,JMIN,JMAX
      COMMON /LIMITS/ XMIN,XMAX,YMIN,YMAX
      XV=XMIN
      YV=YMIN
      IF(XT.EQ.0.)GO TO 30
      DO 10 I=1,NXTX
      J=I-NXTX
      IF(J*XDEL.GT.XMIN)GO TO 20
   10 CONTINUE
   20 XV=J*XDEL
      CALL MOVEA(XV,YMIN)
      CALL DRWREL(0,15)
      CALL MOVEA(XV,YMAX)
      CALL DRWREL(0,-15)
      NXTX=NXTX-1
   30 DO 40 I=1,NXTX
      XVAL=XV+I*XDEL
      CALL MOVEA(XVAL,YMIN)
      CALL DRWREL(0,15)
      CALL MOVEA(XVAL,YMAX)
   40 CALL DRWREL(0,-15)
      IF(YT.EQ.0.)GO TO 70
      DO 50 I=1,NYTX
      K=I-NYTX
      IF(K*YDEL.GT.YMIN)GO TO 60
   50 CONTINUE
   60 YV=K*YDEL
      CALL MOVEA(XMIN,YV)
      CALL DRWREL(15,0)
      CALL MOVEA(XMAX,YV)
      CALL DRWREL(-15,0)
      NYTX=NYTX-1
   70 DO 80 I=1,NYTX
      YVAL=YV+I*YDEL
      CALL MOVEA(XMIN,YVAL)
      CALL DRWREL(15,0)
      CALL MOVEA(XMAX,YVAL)
   80 CALL DRWREL(-15,0)
      IF(XT.EQ.0.)GO TO 90
      NXTX=NXTX+1
   90 IF(YT.EQ.0.)GO TO 100
      NYTX=NYTX+1
  100 RETURN
      END
C-----------------SUBROUTINE--INPUT---------------------------------
C
      SUBROUTINE INPUT(CAY,Y,CAYO,NPTS,IP,ISH,MP,EZERO,SS,S,IREAD,
     C    IERR)
C * PROGRAM TO SETUP INPUT OF ALL NECESSARY PARAMETERS FOR THE 
C * FITTING ROUTINE
      DIMENSION CAY(1024),Y(1024),CCAY(1024),YY(1024),
     C    CAYO(1024),IP(111),ISH(111),SS(2,12,10),S(2)
      LOGICAL*1 YN,BLANK,ZERO,FLAG,TFLAG,CHKPRT,CHKOUT,
     C    MEMBER,MMBR,TSTNAM,SS,S,READIT,SAVE
      DATA BLANK,ZERO/' ',O0/
      COMMON /CARD1/ MEMBER(10),M,NDS,NITER,NSTEP,DAMPER,DESIRE,Q,
     C    POWER,KPOWER,TFLAG,CHKPRT,CHKOUT
      COMMON /CARD2/ CAYLO,CAYHI
      COMMON /CARD3/ FLAG(122),IFLGVR(13)
      COMMON /CARD4/ ZA(10),Z(10),A(10),R(10),SIGMA(10),B(10),
     C    GAMMA(10),P(10,4),G(10),H(10),EZ
      COMMON /CARD5/ NSHELL(10),NENV(10),MMBR(10,10),QV(10),POWERV(10)
      COMMON /CARD6/ TSTNAM(10)
      COMMON /CARD7/ ENVLOP(100,2,10)
      COMMON /CARD8/ NPS(10),XKAY(20,10),PSI(20,10)
      COMMON /LIMITS/ XMIN,XMAX,YMIN,YMAX
      COMMON /YSET/ YSMIN,YSMAX
      CALL INITT(10)
	CALL TOUTPT(24)		!SWITCH TO ANSI MODE
   10 IER=0
      READIT=.TRUE.
      SAVE=.FALSE.
      TYPE 15
15    FORMAT($' DO YOU WANT TO READ CARD SETUP PARAMETERS FROM DISK? ')
      ACCEPT 20, YN
   20 FORMAT(A1)
      IF(YN.EQ.'Y')CALL STORIT(READIT,SAVE,IREAD,IER)
      IF(IER.EQ.1)GO TO 10
      CALL SETUP(CAY,Y,NPTS,IREAD,IERR)
      IF(IERR.EQ.1)GO TO 130
      EZERO=EZ
      IF(NSTEP.EQ.0) NSTEP=1
      DO 40 J=1,M
      JIJ=J
      IF(A(J).NE.1.0.OR.B(J).NE.1.0.OR.GAMMA(J).NE.1.0) GO TO 40
      FLAG(10+J)=.FALSE.
      FLAG(40+J)=.FALSE.
      FLAG(50+J)=.FALSE.
      IF(MMBR(1,J).EQ.BLANK)GO TO 40
      CALL CASEFL (CCAY,YY,NPTV,MMBR(1,J),GARB1,GARB2,IERR)
      IF(IERR.EQ.1)GO TO 130
      DO 30 JJ=1,NPTV
   30 YY(JJ)=YY(JJ)*QV(J)*CCAY(JJ)**POWERV(J)
      CALL OUTLIN(YY,CCAY,NPTV,JIJ,ENVLOP,NENV)
   40 CONTINUE
      DO 50 J=1,M
      IF(P(J,1).NE.1.) GO TO 50
      IF(P(J,2).NE.1.) GO TO 50
      IF(P(J,3).NE.1.) GO TO 50
      IF(P(J,4).NE.1.) GO TO 50
C * IF ALL FOUR VALUES OF P FOR ANY SHELL ARE EQUAL TO 1.0,
C * THE ENVELOPE FOR THAT SHELL IS BEING FURNISHED AS A TABLE.
C * THE P'S CANNOT VARY.
      FLAG(60+J)=.FALSE.
      FLAG(70+J)=.FALSE.
      FLAG(80+J)=.FALSE.
      FLAG(90+J)=.FALSE.
   50 CONTINUE
      CALL SETIND(M,IP,ISH,MP,SS,S)
C * FORM THE N-A PRODUCT, WHICH IS TREATED AS A SINGLE VARIABLE
C * DURING THE FITTING PROCESS
      DO 60 L=1,M
   60 ZA(L)=Z(L)*A(L)
      VSQDEL=.262467*EZ
      DO 70 L=1,NPTS
      CAYO(L)=CAY(L)
      CAYSQ=CAY(L)*CAY(L)
   70 CAY(L)=SQRT(CAYSQ-VSQDEL)
C * MODIFY THE INPUT DATA SET BY TRUNCATING THE ENDS AT CAYLO
C * AND CAYHI AND PICK ONLY EVERY N"STEPTH" POINT.
   80 JLO=1
      DO 90 J=1,NPTS
      IF(CAYLO.GT.CAY(J)) GO TO 90
      JLO=J
      GO TO 100
   90 CONTINUE
C * FIND THE RANGE ON Y (THE GIVEN VALUES OF CHI) FOR PLOTTING
C * PURPOSES
  100 YSMIN=1.E30
      YSMAX=-YSMIN
      J=1
      DO 110 K=JLO,NPTS,NSTEP
      IF(CAY(K).GT.CAYHI)GO TO 120
      CAY(J)=CAY(K)
      CAYO(J)=CAYO(K)
      Y(J)=Q*Y(K)*CAY(K)**POWER
      IF(Y(J).GE.YSMAX)YSMAX=Y(J)
      IF(Y(J).LE.YSMIN) YSMIN=Y(J)
  110 J=J+1
  120 NPTS=J-1
C * IF NO FIRST GUESS FOR EZERO WHEN VARYING THE EDGE ENERGY
C * THEN CHOOSE ONE
  130 READIT=.FALSE.
      IF(IERR.EQ.1)SAVE=.TRUE.
      CALL STORIT(READIT,SAVE,IREAD,IER)
      YN='N'
C      IF(CHKPRT)GO TO 140
C      TYPE 135
C135	FORMAT($' WOULD YOU ALSO LIKE A SETUP LISTING?')
C      ACCEPT 20,YN
  140 IF(IERR.EQ.1)GO TO 160
      IF(CHKPRT.OR.YN.EQ.'Y')CALL PRNTIT
      DESIRE=DESIRE/100.
      IF(FLAG(121).AND.EZ.EQ.0.) GO TO 150
      GO TO 160
  150 CALL EZGEN(CAY,Y,CAYO,NPTS,EZERO)
  160 CONTINUE
      RETURN
      END
C---------------------------SUBROUTINE--ITER8-------------------------------
C
      SUBROUTINE ITER8(CHI,CAY,Y,CAYO,NPTS,ITERNO,RHS,MP,IP,ISH,
     C    SS,S,VPREV,CHANGE,EZERO,RTURN1)
C * SUBROUTINE TO GENERATE AN ITERATION TOWARD A CLOSER FIT TO THE
C * ORIGINAL DATA
      DIMENSION CHI(1024),CAY(1024),Y(1024),CAYO(1024),COEFF(111,111),
     C    RHS(111),IP(111),ISH(111),SOLN(121,2),SS(2,12,10),S(2)
      LOGICAL*1 SS,S,CHANGE
      REAL*8 CHI,RHS
C * SET MINIMUM SIGMA VALUE
      SIGMIN=-17.27/(CAY(NPTS)**2)
C * GENERATE THE COEFFICIENT MATRIX
      IF(MP.GT.0) CALL COEGEN(CHI,CAY,Y,CAYO,NPTS,IP,ISH,
     C	ITERNO,COEFF,RHS,MP)
C * GENERATE THE SOLUTIONS FOR THIS ITERATION UNTIL THE 
C * STANDARD DEVIATION INCREASES
      CALL SOLVER(ITERNO,SOLN,NBEST,DAMBST,SIGMIN,VPREV,EZERO,COEFF,
     C    RHS,MP,CHI,CAY,Y,CAYO,NPTS,IP,ISH,SS,S,CHANGE,RTURN1)
      CALL ADJUST(SOLN,NBEST,CHI,CAY,Y,NPTS,EZERO,SIGMIN)
      CALL RESULT(VPREV,SS,S,ITERNO)
      RETURN
      END
C------------------SUBROUTINE--OUTLIN-------------------------------------
C
      SUBROUTINE OUTLIN(Y,X,N,JSH,ENVLOP,NENV)
C * SUBROUTINE TO CALCULATE THE ENVELOPE SHAPE OF A REFERENCE SET OF
C * CHI VS K DATA
      DIMENSION X(1024),Y(1024)
      DIMENSION ENVLOP(100,2,10),NENV(10)
      K=1
      ENVLOP(K,1,JSH)=X(1)
      ENVLOP(K,2,JSH)=ABS(Y(1))
      K=2
      NM1=N-1
      DO 10 J=2,NM1
      IF(X(J).GT.X(J+1)) GO TO 20
      IF(Y(J).GT.Y(J-1).AND.Y(J).LT.Y(J+1)) GO TO 10
      IF(Y(J).LT.Y(J-1).AND.Y(J).GT.Y(J+1)) GO TO 10
      ENVLOP(K,1,JSH)=X(J)
      ENVLOP(K,2,JSH)=ABS(Y(J))
      K=K+1
      IF(K.GT.100) GO TO 20
   10 CONTINUE
   20 CONTINUE
      ENVLOP(K,1,JSH)=X(J)
      ENVLOP(K,2,JSH)=ABS(Y(J))
      NENV(JSH)=K
      RETURN
      END
C--------------------SUBROUTINE--OUTPUT---------------------------------
C
      SUBROUTINE OUTPUT(VPREV,S,SS)
C * SUBROUTINE TO GENERATE THE PRINTER OUTPUT OF ALL THE VARIABLES FOR
C * EACH COMPLETE ITERATION
C * THE PARAMETERS WHICH ARE VARIED ARE ENCLOSED IN PARENTHESES
      DIMENSION S(24,10),SS(2)
	LOGICAL*1 FLAG
      COMMON /CARD1/ MEMBER(10),M,NDS,NITER,NSTEP,DAMPER,DESIRE,Q,
     C    POWER,KPOWER,TFLAG,CHKPRT,CHKOUT
      COMMON /CARD3/ FLAG(122),IFLGVR(13)
      COMMON /CARD4/ ZA(10),Z(10),A(10),R(10),SIGMA(10),BETA(10),
     C    GAMMA(10),P(10,4),G(10),H(10),EZ
      LOGICAL*1 MEMBER,TFLAG,CHKPRT,CHKOUT,SS,S,AFTER
C	CALL NEWPAG
	CALL TOUTPT(24)		!SWITCH TO ANSI MODE
	KUNIT=5
      WRITE(KUNIT,30) (MEMBER(I),I=1,10)
   30 FORMAT('  EXAFIT PROGRAM TO FIT ANALYTICAL MODEL FUNCTION TO ',
     C    10A1,' EXAFS DATA')
      WRITE(KUNIT,40) M,SS(1),EZ,SS(2)
   40 FORMAT('0NOMINAL PARAMETERS      SHELL ',I2,10X,'EDGE ENERGY ',
     C    'SHIFT ',A1,1PE12.5,A1)
      WRITE(KUNIT,50) (L,L=1,M)
50	FORMAT(' SHELL #:',10(7X,I2,6X))
	WRITE(KUNIT,60)IFLGVR(1),(S(1,L),Z(L),S(2,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(2),(S(3,L),A(L),S(4,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(3),(S(5,L),R(L),S(6,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(4),(S(7,L),SIGMA(L),S(8,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(5),(S(9,L),BETA(L),S(10,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(6),(S(11,L),GAMMA(L),S(12,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(7),(S(13,L),P(L,1),S(14,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(8),(S(15,L),P(L,2),S(16,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(9),(S(17,L),P(L,3),S(18,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(10),(S(19,L),P(L,4),S(20,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(11),(S(21,L),G(L),S(22,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(12),(S(23,L),H(L),S(24,L),L=1,M)
60	FORMAT(1X,3X,A2,' : ',10(1X,A1,1PE12.5,A1))
      WRITE(KUNIT,70)
   70 FORMAT('0** PARENTHESES INDICATE ESTIMATED VALUES TO BE REFINED
     C BY THIS PROGRAM **',//,' ** IMPROVED VALUES ARE TABULATED
     C BELOW **')
	TYPE *, ' '
	TYPE *, 'HIT "ENTER" TO CONTINUE !'
	ACCEPT 25,IAA
25	FORMAT(I1)
   	RETURN
	END
C--------------------SUBROUTINE--PRNTIT--------------------------------
C
      SUBROUTINE PRNTIT
C * SUBROUTINE TO STORE OR READ THE CARD SETUP PARAMETERS
      LOGICAL*1 FLAG,TFLAG,CHKPRT,CHKOUT,MEMBER,MMBR,
     C    TSTNAM,NOMEM(10)
      DIMENSION IDAT(5)
      COMMON /CARD1/ MEMBER(10),M,NDS,NITER,NSTEP,DAMPER,DESIRE,
     C    Q,POWER,KPOWER,TFLAG,CHKPRT,CHKOUT
      COMMON /CARD2/ CAYLO,CAYHI
      COMMON /CARD3/ FLAG(122),IFLGVR(13)
      COMMON /CARD4/ VAL(131)
      COMMON /CARD5/ NSHELL(10),NENV(10),MMBR(10,10),QV(10),
     C    POWERV(10)
      COMMON /CARD6/ TSTNAM(10)
      COMMON /CARD7/ ENVLOP(100,2,10)
      COMMON /CARD8/ NPS(10),XKAY(20,10),PSI(20,10)
      DATA NOMEM/10*' '/,NONENV/0/,GARB/0./
	CALL NEWPAG
	CALL TOUTPT(24)		!SWITCH TO ANSI MODE
      CALL DATE(IDAT)
	KUNIT=5 
      WRITE(KUNIT,80) IDAT 
      WRITE(KUNIT,90) (MEMBER(I),I=1,10),M,NDS,NITER,NSTEP,DAMPER,
     C    DESIRE,Q,POWER,KPOWER,TFLAG,CHKPRT,CHKOUT
      WRITE(KUNIT,100) CAYLO,CAYHI
      WRITE(KUNIT,110)(LLL,LLL=1,M)
	DO 11 I=1,12
11      WRITE(KUNIT,120) IFLGVR(I),(FLAG((I-1)*10+II),II=1,M)
	WRITE(KUNIT,121) IFLGVR(13),FLAG(121)
	TYPE *,' '
	TYPE *,'HIT "ENTER" TO CONTINUE '
	ACCEPT 5,IA
5	FORMAT(I1)
	WRITE(KUNIT,125)
      DO 10 I=1,12
   10 WRITE(KUNIT,130) IFLGVR(I),(VAL(I*10+J),J=1,M)
      WRITE(KUNIT,140) IFLGVR(13),VAL(131)
      WRITE(KUNIT,150) 
      DO 30 I=1,M
      IF(VAL(20+I).NE.1..OR.VAL(50+I).NE.1..OR.VAL(60+I).NE.1.)
     C    GO TO 30
      IF(MMBR(1,I).NE.' ') GO TO 20
      WRITE(KUNIT,160) NSHELL(I),NONENV,(NOMEM(J),J=1,10),GARB,GARB
      GO TO 30
   20 WRITE(KUNIT,160) NSHELL(I),NENV(I),(MMBR(J,I),J=1,10),
     C    QV(I),POWERV(I)
   30 CONTINUE
      WRITE(KUNIT,170)
      IF(.NOT.TFLAG)GO TO 40
      WRITE(KUNIT,180) (TSTNAM(I),I=1,10)
40      WRITE(KUNIT,190)
      DO 50 J=1,M
      IF(VAL(20+J).NE.1..OR.VAL(50+J).NE.1..OR.VAL(60+J).NE.1.)
     C    GO TO 50
      IF(MMBR(1,J).NE.' ') GO TO 50
      WRITE(KUNIT,200) NSHELL(J),NENV(J)
      WRITE(KUNIT,210) (ENVLOP(K,1,J),ENVLOP(K,2,J),K=1,NENV(J))
   50 CONTINUE
      WRITE(KUNIT,220) 
      DO 70 J=1,M
      DO 60 I=10,40,10
   60 IF(VAL(60+I+J).NE.1.)GO TO 70
      WRITE(KUNIT,230) NSHELL(J),NPS(J)
      WRITE(KUNIT,240) (XKAY(N,J),PSI(N,J),N=1,NPS(J))
   70 CONTINUE
   80 FORMAT('    ***** EXAFIT PROGRAM  ***** ',
     C    5A2,' INPUT CARD PARAMETERS FROM SETUP')
   90 FORMAT(' CARD 1      ;MEMBER TO FIT :',1X,10A1,/,
     C    2X,4I5,1P4E10.3,I5,2X,3L1)
  100 FORMAT(' CARD 2      ;RANGE OF FIT  KMIN, KMAX',
     C    /,5X,1PE10.3,1X,2PE10.3)
110   FORMAT(' CARD 3      ;PARAMETERS BEING VARIED (T=TRUE   ',
     C    'F=FALSE)',/,' SHELL NO: ',10(3X,I2))
120   FORMAT(1X,3X,A2,4X,':',10(4X,L1))
121   FORMAT(1X,3X,A2,4X,':',4X,L1)
125   FORMAT(/,' CARD 4      ;VALUES OF THE FIT PARAMETERS')
  130 FORMAT(1X,A2,3X,1P10E12.5)
  140 FORMAT(1X,A2,3X,1PE12.5)
  150 FORMAT(' CARD 5      ;ENVELOPE MEMBER INPUT FROM DISK')
  160 FORMAT(1X,2I5,1X,10A1,1X,1P2E10.3)
  170 FORMAT(' CARD 6      ;TESTCASE MEMBER')
  180 FORMAT(1X,10A1)
  190 FORMAT(' CARD 7      ;ENVELOPE PAIRS ENTERED FROM TERMINAL')
  200 FORMAT(' SHELL',I5,'   NO. OF ENVELOPE PAIRS',I5,/)
  210 FORMAT(1X,1P8E12.4)
  220 FORMAT(' CARD 8      ;PHASE PAIRS ENTERED FROM TERMINAL')
  230 FORMAT(' SHELL',I5,'    NO. OF PHASE PAIRS',I5,/)
  240 FORMAT(1X,1P8E12.4)
  250 FORMAT(/)
      TYPE *,' '
	TYPE *,'HIT "ENTER" TO CONTINUE'
	ACCEPT 5,IA
	RETURN
      END
C--------------------SUBROUTINE--RESULT---------------------------------
C
      SUBROUTINE RESULT(VPREV,S,SS,ITERNO)
C * SUBROUTINE TO PRINT THE BEST FIT VALUES OF ALL THE VARIABLES FOR
C * EACH COMPLETE ITERATION
C * THE PARAMETERS WHICH ARE VARIED ARE ENCLOSED IN PARENTHESES
      DIMENSION S(24,10),SS(2)
      LOGICAL*1 SS,S
      LOGICAL*1 FLAG,TFLAG,CHKPRT,CHKOUT,MEMBER,MMBR,
     C    TSTNAM,READIT,SAVE,YN,FLN,EXT
      COMMON /CARD1/ MEMBER(10),M,NDS,NITER,NSTEP,DAMPER,DESIRE,
     C    Q,POWER,KPOWER,TFLAG,CHKPRT,CHKOUT
      COMMON /CARD2/ CAYLO,CAYHI
      COMMON /CARD3/ FLAG(122),IFLGVR(13)
      COMMON /CARD4/ ZA(10),Z(10),A(10),R(10),SIGMA(10),BETA(10),
     C    GAMMA(10),P(10,4),G(10),H(10),EZ
      COMMON /CARD5/ NSHELL(10),NENV(10),MMBR(10,10),QV(10),
     C    POWERV(10)
      COMMON /CARD6/ TSTNAM(10)
      COMMON /CARD7/ ENVLOP(100,2,10)
      COMMON /CARD8/ NPS(10),XKAY(20,10),PSI(20,10)
      COMMON /SUNAM/ FLN(10),EXT(6)
	KUNIT=5
C	CALL NEWPAG
	CALL TOUTPT(24)		!SWITCH TO ANSI MODE
      WRITE(KUNIT,10) ITERNO
   10 FORMAT(' THE BEST RESULTS FOR ITERATION ',I2,' ARE: ')
      WRITE(KUNIT,20) SS(1),EZ,SS(2)
   20 FORMAT(5X,' THE EDGE ENERGY SHIFT IS :',A1,1PE12.5,A1)
      WRITE(KUNIT,50) (L,L=1,M)
50	FORMAT(' SHELL #:',10(7X,I2,6X))
	WRITE(KUNIT,60)IFLGVR(1),(S(1,L),Z(L),S(2,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(2),(S(3,L),A(L),S(4,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(3),(S(5,L),R(L),S(6,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(4),(S(7,L),SIGMA(L),S(8,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(5),(S(9,L),BETA(L),S(10,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(6),(S(11,L),GAMMA(L),S(12,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(7),(S(13,L),P(L,1),S(14,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(8),(S(15,L),P(L,2),S(16,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(9),(S(17,L),P(L,3),S(18,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(10),(S(19,L),P(L,4),S(20,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(11),(S(21,L),G(L),S(22,L),L=1,M)
	WRITE(KUNIT,60)IFLGVR(12),(S(23,L),H(L),S(24,L),L=1,M)
60	FORMAT(1X,3X,A2,' : ',10(1X,A1,1PE12.5,A1))
C
      STDV=VPREV**(0.5)
      WRITE(KUNIT,40) STDV
   40 FORMAT('  THE STANDARD DEVIATION OF THIS FIT IS :'1PE12.4)
C
	OPEN(UNIT=4,FILE='SETUP.6',FORM='UNFORMATTED',
     C    STATUS='UNKNOWN')
      WRITE(4) MEMBER,M,NDS,NITER,NSTEP,NITER,DAMPER,
     C    DESIRE,Q,POWER,KPOWER,TFLAG,CHKPRT,CHKOUT,CAYLO,CAYHI
      WRITE(4) FLAG,IFLGVR,ZA,Z,A,R,SIGMA,BETA,GAMMA,P,G,H,EZ,
     C    NSHELL,NENV,MMBR,QV,POWERV,TSTNAM,ENVLOP,NPS,XKAY,PSI
	CLOSE(UNIT=4)
	TYPE *,'  HIT "ENTER" TO CONTINUE'
	ACCEPT 100, IIII
100	FORMAT(I1)
	RETURN
      END
C---------------SUBROUTINE--SCPLT----------------------------------
C
      SUBROUTINE SCPLT(I)
      COMMON /PLOT/ XDEL,NEX,NXTX,YDEL,NEY,NYTX,IMIN,IMAX,JMIN,JMAX
      COMMON /LIMITS/ XMIN,XMAX,YMIN,YMAX
      IF(I-2) 10,20,30
   10 CALL TWINDO (IMIN,IMAX,JMIN,JMAX)
      GO TO 50
   20 CALL DWINDO (XMIN,XMAX,YMIN,YMAX)
      GO TO 50
   30 IF(I.EQ.4)GO TO 40
      CALL HDCOPY
      GO TO 50
   40 CALL MOVEA(XMIN,YMIN)
      CALL DRAWA(XMAX,YMIN)
      CALL DRAWA(XMAX,YMAX)
      CALL DRAWA(XMIN,YMAX)
      CALL DRAWA(XMIN,YMIN)
      CONTINUE
   50 RETURN
      END
C-----------------SUBROUTINE--SETIND-----------------------------------
C
      SUBROUTINE SETIND(M,IP,ISH,MP,SS,S)
C * SUBROUTINE TO SET INDICES FOR EACH PARAMETER BEING VARIED
      DIMENSION IP(111),ISH(111),SS(2,12,10),S(2)
      COMMON /CARD3/ FLAG(122),IFLGVR(13)
      COMMON /CARD4/ ZA(10),Z(10),A(10),R(10),SIGMA(10),B(10),
     C    GAMMA(10),P(10,4),G(10),H(10),EZ
      LOGICAL*1 FLAG,LEFT,RIGHT,BLANK,SS,S
      DATA LEFT/'('/,RIGHT/')'/,BLANK/' '/
      DO 10 K=1,10
      DO 10 J=1,12
      SS(1,J,K)=BLANK
   10 SS(2,J,K)=BLANK
      S(1)=BLANK
      S(2)=BLANK
      J=0
      DO 30 K=1,M
      IF(.NOT.FLAG(K).AND..NOT.FLAG(K+10))GO TO 30
   20 J=J+1
      IP(J)=1
      ISH(J)=K
      IF(FLAG(K))SS(1,1,K)=LEFT
      IF(FLAG(K))SS(2,1,K)=RIGHT
      IF(FLAG(K+10))SS(1,2,K)=LEFT
      IF(FLAG(K+10))SS(2,2,K)=RIGHT
   30 CONTINUE
      DO 50 N=20,110,10
      DO 40 K=1,M
      IF(.NOT.FLAG(K+N))GO TO 40
      J=J+1
      IP(J)=N/10
      ISH(J)=K
      SS(1,IP(J)+1,K)=LEFT
      SS(2,IP(J)+1,K)=RIGHT
   40 CONTINUE
   50 CONTINUE
      IF(.NOT.FLAG(121))GO TO 60
      J=J+1
      IP(J)=12
      S(1)=LEFT
      S(2)=RIGHT
      ISH(J)=1
   60 CONTINUE
      MP=J
      RETURN
      END
C-----------------------------SUBROUTINE--SOLVER-----------------------------
C
      SUBROUTINE SOLVER(ITERNO,SOLN,NBEST,DAMBST,SIGMIN,VPREV,EZERO,
     C    COEFF,RHS,MP,CHI,CAY,Y,CAYO,NPTS,IP,ISH,SS,S,CHANGE,RTURN1)
C * SUBROUTINE TO OBTAIN TWO DIFFERENT SOLUTIONS,AND A STANDARD
C * DEVIATION ASSOCIATED WITH EACH
C * THERE WILL BE A BEST AND SECOND-BEST SOLUTION FOR EACH ITERATION
      COMMON /CARD1/ MEMBER(10),M,NDS,NITER,NSTEP,DAMPER,DESIRE,Q,
     C    POWER,KPOWER,TFLAG,CHKPRT,CHKOUT
      COMMON /CARD3/ FLAG(122),IFLGVR(13)
      COMMON /CARD4/ VALS(131)
      DIMENSION CHI(1024),CAY(1024),Y(1024),CAYO(1024),IP(111),ISH(111),
     C    ANSWER(111),SOLN(121,2),COEFF(MP,MP),RHS(111),SS(2,12,10),
     C    S(2)
      REAL*4 PVAL(121)
      LOGICAL*1 CHANGE,CHNG(2),MEMBER,TFLAG,CHKPRT,CHKOUT,FLAG,SS,S
      REAL*8 CHI,ANSWER,R,RHS
      R=1.0D0
	KUNIT=6
      R=1.8D0-R
      PI=3.14159
      TWOPI=6.28318
      NBEST=1
      DMPMIN=1.0
      JJ=99
      VN1=1.E30
      IF(ITERNO.GT.1) VN1=VPREV
      RTURN1=0.0
      CHNG(1)=.FALSE.
      J=1
      IF(MP.LE.0) GO TO 60
      JJ=0
      GO TO 20
   10 J=3-NBEST
C * APPLY DAMPING FACTOR TO THE COEFFICIENT MATRIX
   20 DAMP=DMPMIN
      IF(JJ.GT.0) DAMP=DMPMIN+(JJ*DAMPER/NDS)**2
      JJ=JJ+1
      DO 30 L=1,MP
   30 COEFF(L,L)=DAMP*COEFF(L,L)
      DO 40 L=1,10
   40 SOLN(L+10,J)=VALS(L)
      DO 50 L=11,111
   50 SOLN(L+10,J)=VALS(L+20)
      IF(JJ.GT.1.OR.ITERNO.GT.1)GO TO 70
C * CALCULATE THE STANDARD DEVIATION OF THE INITIAL GUESS
   60 CALL STDEV(0,0.0,0,0,VALS(11),VZERO,CHI,CAY,Y,
     C    CAYO,NPTS,S,SS,EZERO)
      VPREV=VZERO
      IF(MP.LE.0) GO TO 130
      VN1=VZERO
      CALL NEWPAG
      CALL BSTPLT(0,VZERO,EZERO,Y,NPTS,CHI,CAY,DAMBST)
      IF(NITER.EQ.0) GO TO 150
C * INVERT THE COEFFICIENT MATRIX
   70 MEQNS=MP
      CALL MATINV(COEFF,RHS,ANSWER,MEQNS)
      CHNG(J)=.FALSE.
      DO 80 L=1,MP
C * RESTORE COEFF MATRIX FOR NEXT DAMPING FACTOR APPLICATION
      COEFF(L,L)=COEFF(L,L)/DAMP     
      N=10*IP(L)+ISH(L)
      IF(DABS(ANSWER(L)).GT.ABS(DESIRE*SOLN(N,J)))CHNG(J)=.TRUE.
      SOLN(N,J)=SOLN(N,J)+ANSWER(L)
   80 CONTINUE
      DO 110 L=1,M
C * SORT OUT N AND A, GIVEN THEIR PRODUCT AND WHICH OF THEM IS BEING VARIED
C * IF N AND A ARE BOTH VARIABLE, SET N EQUAL TO 1.0
      IF(FLAG(L+10).AND.FLAG(L)) SOLN(L,J)=1.0
      IF(.NOT.FLAG(L+10)) SOLN(L,J)=SOLN(L+10,J)/VALS(L+20)
      IF(.NOT.FLAG(L)) SOLN(L,J)=VALS(L+10)
      SOLN(L+10,J)=SOLN(L+10,J)/SOLN(L,J)
C * MAKE "P0" POSITIVE AND LESS THAN TWOPI
      SOLN(L+60,J)=ABS(SOLN(L+60,J))
   90 IF(SOLN(L+60,J).LT.TWOPI)GO TO 100
      SOLN(L+60,J)=SOLN(L+60,J)-TWOPI
      GO TO 90
C * IF SIGMA IS NEGATIVE, LIMIT THE VALUE SUCH THAT
C *    EXP(-2.*SIGMA*KMAX**2) < 10.E15
  100 IF(SOLN(L+30,J).GE.SIGMIN)GO TO 110
      SOLN(L+30,J)= SIGMIN
  110 CONTINUE
      DO 120 I=1,121
  120 PVAL(I)=SOLN(I,J)
      CALL STDEV(ITERNO,DAMP,J,JJ,PVAL,VARN,CHI,CAY,Y,
     C    CAYO,NPTS,S,SS,EZERO)
      IF(VARN.GT.VPREV) GO TO 140
      VPREV=VARN
  130 NBEST=J
      DAMBST=DAMP
      IF(CHNG(J).AND.JJ.LT.NDS)GO TO 10
  140 IF(JJ.LT.NDS) GO TO 10
      CHANGE=CHNG(NBEST)
      EZERO=SOLN(121,NBEST)
      IF(VPREV.LT.VN1) GO TO 150
      WRITE(KUNIT,180)
C      REWIND KUNIT
      RTURN1=1.
      GO TO 160
  150 CONTINUE
      VN1=VPREV
      IF(.NOT.CHANGE)WRITE(KUNIT,190) DESIRE
 160    CALL NEWPAG
      CALL BSTPLT(ITERNO,VPREV,EZERO,Y,NPTS,CHI,CAY,DAMBST)
      IF(MP.LE.0)RTURN1=1.
  170 RETURN
  180 FORMAT( ' DIVERGING SOLUTION')
  190 FORMAT(' AT THIS ITERATION, ALL PARAMETERS CHANGED BY LESS THAN
     C   THE DESIRED PERCENTAGE (',2PF8.4,'%)')
      END
C----------------------SUBROUTINE--STDEV-------------------------------
C
      SUBROUTINE STDEV(ITERNO,DAMP,J,JJ,P,VARN,CHI,CAY,
     C    Y,CAYO,NPTS,SS,S,EZERO)
C * SUBROUTINE TO EVALUATE THE STANDARD DEVIATION FOR EACH SOLUTION
C * FROM SOLVER
      COMMON /CARD1/ MEMBER(10),M,NDS,NITER,NSTEP,DAMPER,DESIRE,Q,
     C    POWER,KPOWER,TFLAG,CHKPRT,CHKOUT
      COMMON /CARD3/ FLAG(122),IFLGVR(13)
      DIMENSION CHI(1024),CAY(1024),Y(1024),CAYO(1024),
     C    S(24,10),SS(2)
      REAL*8 CHI,SUM1,SUM2,QQ
      REAL*4 P(121)
      LOGICAL*1 MEMBER,TFLAG,CHKPRT,CHKOUT,FLAG,SS,S
      KK=0
	KUNIT=5
C * IF THE "INNER POTENTIAL ENERGY",(EZERO) IS BEING VARIED,
C * TRANFORM THE "K" ARRAY
   10 IF(.NOT.FLAG(121)) GO TO 30
      VSQ=.262467
      DELTA=P(121)-EZERO
      KK=KK+1
      VSQDEL=VSQ*DELTA
      DO 20 L=1,NPTS
      CAYSQ=CAY(L)*CAY(L)
      IF(CAYSQ.GT.VSQDEL)GO TO 20
      WRITE(KUNIT,80) L,CAY(L),DELTA,VSQ
      STOP
   20 CAY(L)=SQRT(CAYSQ-VSQDEL)
   30 SUM1=0.D0
      SUM2=0.D0
C *      IF(CHKPRT) WRITE(KUNIT,1105)
C * GENERATE THE "VARIANCE BY SUMMING OVER ALL THE POINTS IN THE ARRAY
      DO 40 N=1,NPTS
      NAN=N
      CALL CHIGEN(NAN,P( 1),P(11),P(41),P(51),P(21),P(31),P(61),
     C    P(101),P(111),M,KPOWER,CHI,CAY,Y,CAYO,NPTS)
      QQ=CHI(N)-Y(N)
      SUM1=SUM1+QQ*QQ
      SUM2=SUM2+1
C *      IF(MOD(N-1,5).EQ.0.AND.CHKPRT)
C *     C    WRITE(KUNIT,1106)N,CAY(N),CHI(N),Y(N),QQ
   40 CONTINUE
      MN=0
      VARN=SUM1/(SUM2-MN)
      IF(.NOT.CHKPRT) GO TO 50
C *      STDV=SQRT( ABS(VARN))
C *      WRITE(6,1107) STDV
C * IF THE "INNER POTENTIAL ENERGY" IS BEING VARIED, RESTORE THE
C * "K" ARRAY PROIR TO PERFORMING THE SOLUTION WITH THE NEXT 
C * DAMPING FACTOR
   50 IF(.NOT.FLAG(121))GO TO 70
      DO 60 L=1,NPTS
      CAYSQ=CAY(L)*CAY(L)
   60 CAY(L)=SQRT(CAYSQ+VSQDEL)
 70     RETURN
   80 FORMAT(' *** UNACCEPTABLE TRANSFORM *** ',/,
     C    ' L,CAY(L),DELTA,VSQ ',I5,1P3E12.4)
      END
C----------------SUBROUTINE--STORIT----------------------------------
C
      SUBROUTINE STORIT(READIT,SAVE,IREAD,IER)
C * SUBROUTINE TO STORE OR READ THE CARD SETUP PARAMETERS
      LOGICAL*1 FLAG,TFLAG,CHKPRT,CHKOUT,MEMBER,MMBR,
     C    TSTNAM,READIT,SAVE,YN,FLN,EXT
      COMMON /CARD1/ MEMBER(10),M,NDS,NITER,NSTEP,DAMPER,DESIRE,
     C    Q,POWER,KPOWER,TFLAG,CHKPRT,CHKOUT
      COMMON /CARD2/ CAYLO,CAYHI
      COMMON /CARD3/ FLAG(122),IFLGVR(13)
      COMMON /CARD4/ ZA(10),Z(10),A(10),R(10),SIGMA(10),B(10),
     C    GAMMA(10),P(10,4),G(10),H(10),EZ
      COMMON /CARD5/ NSHELL(10),NENV(10),MMBR(10,10),QV(10),
     C    POWERV(10)
      COMMON /CARD6/ TSTNAM(10)
      COMMON /CARD7/ ENVLOP(100,2,10)
      COMMON /CARD8/ NPS(10),XKAY(20,10),PSI(20,10)
      COMMON /SUNAM/ FLN(10),EXT(6)
      DATA FLN/'S','E','T','U','P','.','0',3*O0/
      DATA EXT/'1','2','3','4','5','6'/
      IF(READIT)GO TO 50
      IF(.NOT.SAVE) GO TO 10
      FLN(7)=EXT(6)
      TYPE 100
      CALL IOWAIT(5)
      GO TO 40
   10 IF(CHKOUT)GO TO 20
      TYPE *,'AFTER ALL YOUR WORK TO SETUP THIS FIT,'
      TYPE 15
15    FORMAT($'DO YOU WANT TO SAVE THE CARD SETUP PARAMETERS? ')
      ACCEPT 90,YN
      IF(YN.NE.'Y')GO TO 80
   20 TYPE 25
25    FORMAT($' DO YOU WANT TO LOOK AT A TABLE OF '
     1        'SETUP DATA FILES? ')
      ACCEPT 90,YN
      IF(YN.EQ.'Y')CALL TABLE(IER)
   30 TYPE 35
35    FORMAT($' TYPE NUMBER OF SETUP FILE FOR '
     1         'OUTPUT TO DISK (1-5) ')
      ACCEPT *, ISUNUM
      IF(ISUNUM.LT.1 .OR. ISUNUM.GT.5)GO TO 30
      FLN(7)=EXT(ISUNUM)
40    OPEN(UNIT=3,FILE=FLN,STATUS='UNKNOWN',FORM='UNFORMATTED')
      WRITE(3) MEMBER,M,NDS,NITER,NSTEP,NITER,DAMPER,
     C    DESIRE,Q,POWER,KPOWER,TFLAG,CHKPRT,CHKOUT,CAYLO,CAYHI
      WRITE(3) FLAG,IFLGVR,ZA,Z,A,R,SIGMA,B,GAMMA,P,G,H,EZ,
     C    NSHELL,NENV,MMBR,QV,POWERV,TSTNAM,ENVLOP,NPS,XKAY,PSI
	CLOSE(3)
      GO TO 80
   50 TYPE 55
55    FORMAT($' DO YOU WANT TO LOOK AT A TABLE OF '
     1    'SETUP DATA FILES? ')
      ACCEPT 90,YN
      IF(YN.EQ.'Y')CALL TABLE(IER)
      IF(IER.NE.1)GO TO 60
      TYPE *,'NO SETUP FILES ON DEFAULT DISK'
      GO TO 80
   60 TYPE 65
65	FORMAT($' TYPE NUMBER OF SETUP FILE FOR INPUT (1-5) ')
      ACCEPT *, ISUNUM
      IF(ISUNUM.LT.1 .OR. ISUNUM.GT.5)GO TO 60
      FLN(7)=EXT(ISUNUM)
	OPEN(UNIT=3,FILE=FLN,STATUS='OLD',FORM='UNFORMATTED')	
      READ(3,ERR=70) MEMBER,M,NDS,NITER,NSTEP,NITER,DAMPER,
     C    DESIRE,Q,POWER,KPOWER,TFLAG,CHKPRT,CHKOUT,CAYLO,CAYHI
      READ(3) FLAG,IFLGVR,ZA,Z,A,R,SIGMA,B,GAMMA,P,G,H,EZ,
     C    NSHELL,NENV,MMBR,QV,POWERV,TSTNAM,ENVLOP,NPS,XKAY,PSI
	CLOSE(3)
      IREAD=1
      GO TO 80
   70 TYPE *,'SETUP FILE',ISUNUM,' NOT FOUND ON DISK'
      GO TO 50
   80 RETURN
   90 FORMAT(A1)
  100 FORMAT(' ABNORMAL TERMINATION: SETUP DATA RECOVERED IN SETUP',
     C   ' FILE 6',/,' AFTER EXITING THE EXAFIT PROGRAM, RENAME THIS', 
     C   ' FILE OR DELETE IT')
      END
C-------------------SUBROUTINE--TABLE-----------------------------------
C
      SUBROUTINE TABLE(IER)
C * SUBROUTINE TO LIST A TABLE OF PRESENT SETUP DATA FILES
      DIMENSION NAME(10),NAMES(10,6)
      LOGICAL*1 UN(10),NAMES,NAME,FLN,EXT,NOMES,MES
      COMMON /SUNAM/ FLN(10),EXT(6)
      DATA NOMES/.FALSE./,MES/.TRUE./
      DATA UN/'<',' ','U','N','U','S','E','D',' ','>'/
      ICOUNT=0
      DO 50 I=1,5
      FLN(7)=EXT(I)
      FLN(8)="0
      FLN(9)="0
      FLN(10)="0
      OPEN(UNIT=3,FILE=FLN,STATUS='OLD',FORM='UNFORMATTED',ERR=20)
      READ(3,ERR=20)NAME
      DO 10 J=1,10
      IF(NAME(J).EQ.' ' .OR. NAME(J).EQ."0)NAME(J)=' '
   10 NAMES(J,I)=NAME(J)
      GO TO 40
   20 ICOUNT=ICOUNT+1
      DO 30 J=1,10
   30 NAMES(J,I)=UN(J)
   40 CLOSE(3)
   50 CONTINUE
      TYPE 70
      DO 60 I=1,5
   60 TYPE 80, I,(NAMES(J,I),J=1,10)
      I=6
      TYPE 90
      TYPE 80, I
      IF(ICOUNT.EQ.5)IER=1
      RETURN
   70 FORMAT(' TABLE OF CURRENT SETUP FILES ON DEFAULT DISK',/)
   80 FORMAT(' SETUP FILE',I2,8X,10A1)
   90 FORMAT(' SAVE SETUP FILE FOR ABNORMAL TERMINATIONS:')
      END
C-----------------------FUNCTION--TBLP-----------------------------------
C
      FUNCTION TBLP(X,Y,XN,XSAV,NV,NSAV,NDEX)
C * SUBROUTINE TO LINEARLY INTERPOLATE A TABLE OF X,Y PAIRS TO FIND THE
C * VALUE OF THE FUNCTION Y AT X=XN
      DIMENSION X(NV),Y(NV)
      DO 10 I=1,NV-1
      IF(X(I).LT.XN.AND.X(I+1).GE.XN)GO TO 20
   10 CONTINUE
   20 IJ=I
      SLOPE=(Y(IJ+1)-Y(IJ))/(X(IJ+1)-X(IJ))
      TBLP=(XN-X(IJ))*SLOPE+Y(IJ)
      RETURN
      END
C-----------------FUNCTION--TBLSL----------------------------------------
C
      FUNCTION TBLSL(X,Y,NV,CAYN)
C * SUBROUTINE FUNCTION TO RETURN SLOPE OF TABLE VALUES
      DIMENSION X(NV),Y(NV)
      DO 10 I=1,NV-1
      IF(X(I).LT.CAYN.AND.X(I+1).GE.CAYN)GO TO 20
   10 CONTINUE
   20 IJ=I
      TBLSL=(Y(IJ+1)-Y(IJ))/(X(IJ+1)-X(IJ))
      RETURN
      END

C-------------------SUBROUTINE--TITLE--------------------------------------
C
C * SUBROUTINE TO LABEL AND DATE THE PLOT
      SUBROUTINE TITLE
      LOGICAL*1 BLANK,ZERO,MEMBER,TFLAG,CHKPRT,CHKOUT
      COMMON /CARD1/ MEMBER(10),M,NDS,NITER,NSTEP,DAMPER,DESIRE,Q,
     C    POWER,KPOWER,TFLAG,CHKPRT,CHKOUT
      DIMENSION ICHAR(5)
      DATA ZERO,BLANK/O0,' '/
C * PUT FILE NAME AND DATE ON PLOT
	CALL ALFMOD
      CALL DATE(ICHAR)
       DO 20 I=1,10
      IF(MEMBER(I).EQ.ZERO) MEMBER(I)=BLANK
   20 CONTINUE
      CALL HOME
      TYPE 40, MEMBER,ICHAR
   40 FORMAT(1X,10A1,/,1X,5A2)
      RETURN
      END
C-------------------SUBROUTINE--XLABEL--------------------------------------
C
      SUBROUTINE XLABEL(XT,XV)
C * SUBROUTINE TO ANNOTATE X AXIS OF PLOT
      LOGICAL*1 CHAR1,CHAR2,CHAR3,CHAR,BLANK,ZERO,ZRO,KCHAR(3)
      COMMON /PLOT/XDEL,NEX,NXTX,YDEL,NEY,NYTX,
     1IMIN,IMAX,JMIN,JMAX
      COMMON /LIMITS/ XMIN,XMAX,YMIN,YMAX
      DATA IXCAP /'W','A','V','E',' ',
     1      'V','E','C','T','O','R',' ',18*' '/
      DATA NXCAP /12/
      DIMENSION IADE(10),IXCAP(30),ICHAR(5),CHAR1(60),
     1CHAR2(5),CHAR3(3),IDUM(30),CHAR(10),ZRO(2)
      EQUIVALENCE(IDUM(1),CHAR1(1))
      DATA ZRO/'0','.'/
      DATA CHAR2/'X','1','0','L','N'/
      DATA ZERO,BLANK/O0,' '/
      KPW=-KPOWER
      IONE=LINWDT(1)
      JONE=LINHGT(1)
      JJ=0
      I=1
C * PUT X-AXIS TICK VALUES ON PLOT
      NHGT=JMIN-2*JONE
      IF(XT.EQ.1.)GO TO 20
      XVAL=XMIN
      GO TO 80
   10 CALL MOVABS(IMIN-(IONE*NCHAR/2),NHGT)
      CALL ANSTR(NCHAR,IADE)
   20 XVAL=XMIN+I*XDEL
      IF(XT.EQ.1.)XVAL=XV+(I-1)*XDEL
      JJ=1
      GO TO 80
   30 CALL MOVEA(XVAL,YMIN)
      CALL MOVREL(-(IONE*NCHAR/2),NHGT-JMIN)
      CALL ANSTR(NCHAR,IADE)
      IF(I.EQ.NXTX)GO TO 40
      I=I+1
      GO TO 20
C * PUT X-AXIS LABEL ON PLOT
   40 IEX=6
      IF(NEX.EQ.0)IEX=0
      NWIDE=IONE*(NXCAP+IEX)
      NHGT=JMIN-4*JONE
      CALL MOVABS(((IMIN+IMAX-NWIDE)/2),NHGT)
      DO 50 I=1,30
      IDUM(I)=IXCAP(I)
   50 CONTINUE
      CALL ANSTR(NXCAP,IDUM)
      IF(NEX.EQ.0)GO TO 210
      DO 60 I=1,3
   60 ICHAR(I)=CHAR2(I)
      CALL ANSTR(3,ICHAR)
      N=-NEX
      GO TO 160
   70 CALL MOVREL(0,JONE/2)
      CALL ANSTR(NCHAR,IADE)
      GO TO 210
C * ENCODE THE TICKS LABELS
   80 XX=XVAL*(10.**NEX)
      NCHAR=0
   90 IF(ABS(XX).LT.1)GO TO 100
      XX=XX/10.
      NCHAR=NCHAR+1
      GO TO 90
  100 IF(XX.LT.0.)NCHAR=NCHAR+1
      IF(XX.EQ.0.)GO TO 130
      XX=XVAL*(10.**NEX)
      ENCODE(10,110,CHAR) XX
  110 FORMAT(F10.0)
      J=9-NCHAR
      NCHAR=NCHAR+1
      DO 120 II=1,NCHAR
  120 IADE(II)=CHAR(II+J)
      GO TO 150
  130 DO 140 II=1,2
  140 IADE(II)=ZRO(II)
      NCHAR=2
  150 IF(JJ.EQ.1)GO TO 30
      GO TO 10
C * SECTION TO ENCODE EXPONENTS
  160 XXX=N
      NCHAR=0
  170 IF(ABS(XXX).LT.0.999)GO TO 180
      XXX=XXX/10.
      NCHAR=NCHAR+1
      GO TO 170
  180 IF(N.LT.0)NCHAR=NCHAR+1
      ENCODE(3,190,KCHAR) N
  190 FORMAT(I3)
      J=3-NCHAR
      DO 200 L=1,NCHAR
  200 IADE(L)=KCHAR(L+J)
      GO TO 70
  210 RETURN
      END
C-------------------SUBROUTINE--YLABEL--------------------------------------
C
      SUBROUTINE YLABEL(YT,YV,KPOWER)
C * SUBROUTINE TO ANNOTATE Y AXIS OF PLOT
      LOGICAL*1 CHAR1,CHAR2,CHAR3,CHAR,BLANK,ZERO,ZRO,KCHAR(3)
	LOGICAL*1 KP(5)
      COMMON /PLOT/XDEL,NEX,NXTX,YDEL,NEY,NYTX,
     1IMIN,IMAX,JMIN,JMAX
      COMMON /LIMITS/ XMIN,XMAX,YMIN,YMAX
      DATA IYCAP /' K','CH','I ',27*'  '/
      DATA NYCAP /5/
      DIMENSION IADE(10),ICHAR(5),CHAR1(60),
     1CHAR2(5),CHAR3(3),IYCAP(30),IDUM(30),CHAR(10),ZRO(2)
      EQUIVALENCE(IDUM(1),CHAR1(1))
      DATA ZRO/'0','.'/
      DATA CHAR2/'X','1','0','L','N'/
      DATA ZERO,BLANK/O0,' '/
	DATA KP/'0','1','2','3','M'/
	KPW=KP(5)
	DO 3 II= 1,4
	   IF(KPOWER .EQ. II-2)KPW=KP(II)
3	CONTINUE
      IONE=LINWDT(1)
      JONE=LINHGT(1)
      JJ=0
      I=1
C * PUT Y-AXIS TICK VALUES ON PLOT
   10 IF(YT.EQ.1.)GO TO 30
      YVAL=YMIN
      GO TO 100
   20 CALL MOVABS(IMIN-IONE*(NCHAR+1),JMIN-(JONE/2))
      CALL ANSTR(NCHAR,IADE)
   30 YVAL=YMIN+I*YDEL
      IF(YT.EQ.1.)YVAL=YV+(I-1)*YDEL
      JJ=1
      GO TO 100
   40 CALL MOVEA(XMIN,YVAL)
      CALL MOVREL(-IONE*(NCHAR+1),-(JONE/2))
      CALL ANSTR(NCHAR,IADE)
      IF(I.EQ.NYTX)GO TO 50
      I=I+1
      GO TO 30
C * PUT Y AXIS LABEL ON PLOT
   50 ILOC=IMIN-IONE*(NCHAR+5)-(IONE/2)
      IF(ILOC.LE.0) GO TO 230
      DO 60 I=1,30
   60 IDUM(I)=IYCAP(I)
      IEY=3
      IF(NEY.EQ.0)IEY=0
      CALL MOVABS(ILOC,(JMIN+JMAX+JONE*(NYCAP+IEY))/2)
	CALL MOVREL(IONE,-JONE/2)
	CALL ANCHO(KPW)
	CALL MOVREL(-2*IONE,JONE/2)
      DO 70 I=1,NYCAP
      ICHAR(1)=CHAR1(I)
      CALL ANCHO(ICHAR(1))
      CALL MOVREL(-IONE,-JONE)
   70 CONTINUE
      IF(NEY.EQ.0)GO TO 230
      N=-NEY
      GO TO 180
   80 CALL MOVREL((-((NCHAR+2)*IONE)/2),-2*JONE)
      CALL SEELOC(IX,IY)
      IF(IY.LE.JMIN)GO TO 230
      IF(IX.LE.0)GO TO 230
      DO 90 I=1,3
   90 ICHAR(I)=CHAR2(I)
      CALL ANSTR(3,ICHAR)
      CALL MOVREL(0,JONE/2)
      CALL ANSTR(NCHAR,IADE)
      GO TO 230
C * SECTION TO ENCODE TICK LABELS
  100 YY=YVAL*(10.**NEY)
      NCHAR=0
  110 IF(ABS(YY).LT.1)GO TO 120
      YY=YY/10.
      NCHAR=NCHAR+1
      GO TO 110
  120 IF(YY.LT.0.)NCHAR=NCHAR+1
      IF(YY.EQ.0.)GO TO 150
      YY=YVAL*(10.**NEY)
      ENCODE(10,130,CHAR) YY
  130 FORMAT(F10.0)
      J=9-NCHAR
      NCHAR=NCHAR+1
      DO 140 II=1,NCHAR
  140 IADE(II)=CHAR(II+J)
      GO TO 170
  150 DO 160 II=1,2
  160 IADE(II)=ZRO(II)
      NCHAR=2
  170 IF(JJ.EQ.1.)GO TO 40
      GO TO 20
C * SECTION TO ENCODE EXPONENTS
  180 YYY=N
      NCHAR=0
  190 IF(ABS(YYY).LT.0.999)GO TO 200
      YYY=YYY/10.
      NCHAR=NCHAR+1
      GO TO 190
  200 IF(N.LT.0)NCHAR=NCHAR+1
      ENCODE(3,210,KCHAR) N
  210 FORMAT(I3)
      J=3-NCHAR
      DO 220 L=1,NCHAR
  220 IADE(L)=KCHAR(L+J)
      GO TO 80
  230 RETURN
      END
C--------------SUBROUTINE--GRAPH-----------------------
C
      SUBROUTINE GRAPH(U,V,W,NPTS,J)
C * THIS SUBROUTINE PLOTS THE DATA IN X-Y FORM
      DIMENSION U(1024),V(1024),W(1024)
      COMMON /LIMITS/ XMIN,XMAX,YMIN,YMAX
      REAL*8 U
      Q1=SNGL(U(1))
      R1=V(1)
      S1=W(1)
      CALL MOVEA(R1,S1)
      DO 10 I=1,NPTS
      R=V(I)
      S=W(I)
   10 CALL DRAWA(R,S)
      CALL MOVEA(R1,Q1)
      DO 30 I=1,NPTS
      R=V(I)
      Q=SNGL(U(I))
      IF(Q.GE.YMAX.OR.Q.LE.YMIN)GO TO 30
	CALL BOXES(R,Q,2)
C      IF(J.EQ.1)GO TO 20
C      CALL DIAMND(R,Q,2)
C      GO TO 30
C   20 CALL BOXES(R,Q,2)
   30 CONTINUE
      RETURN
      END
