C               ********************************
C               * ALEV AND SUPPORTING ROUTINES *
C               ********************************
C
C                  VERSION 1.0    30/04/1989
C
C     THESE SUBROUTINES SUPPORT THE "EV" COMMAND IN ALADDIN.  "ALEV"
C     MUST BE MODIFIED BY THE USER IF THIS COMMAND IS TO RECOGNIZE NEW
C     EVALUATION LABELS (FITTED DATA FORMS).
C
C     ORIGINAL VERSION DESIGNED AND WRITTEN BY R. A. HULSE, OF THE
C     PRINCETON UNIVERSITY, PLASMA PHYSICS LABORATORY.
C     [REFERENCE: "THE ALADDIN ATOMIC PHYSICS DATABASE SYSTEM,"
C      PP. 63-72 IN 'ATOMIC PROCESSES IN PLASMAS,'CONFERENCE
C      PROCEEDINGS 206, EDS. Y.K.KIM AND R.C. ELTON, AMERICAN
C      INSTITUTE OF PHYSICS, NEW YORK, 1990].
C-------------------------------------------------------------------------------
C
C              MODIFICATIONS MADE FOR VERSION 1.0
C
C      REV: 21/4/88 BY R. A. HULSE.
C           RE-ORGANIZE TABULAR OUTPUT AND ALLOW OUTPUT TO FILE
C
C      REV: 15/4/88 BY R. A. HULSE.
C            BREAKOUT CODE FOR EACH DATA TYPE TO SUBROUTINES;
C            GENERATE MULTI-POINT TABULAR OUTPUT;
C            ADD ORNL CHEBYSHEV POLYNOMIAL DATA TYPE
C
C      REV: 7/4/88 BY R. A. HULSE.
C            BASIC VERSION WITH #MEWE AND #TAB1D TYPES
C
C      REV: 1/03/89 BY J.J.SMITH IAEA ATOMIC AND MOLECULAR DATA UNIT
C            SUBROUTINE ALREVD INTRODUCED TO CONTROL THE READING AND
C            PASSING OF INPUT PARAMETERS TO BE USED IN GENERATING
C            TABLES FROM ALEV. SUBROUTINE ALTB1D ADDED TO GENERATE
C            LINEARLY INTERPOLABLE TABLES FOR FITTED DATA FORMS TO
C            REPRODUCE THE FUNCTION TO A INPUT SPECIFIED PERCENTAGE
C            ACCURACY. FUNCTION ALRND6 ADDED TO ROUND REAL NUMBERS TO
C            SIX DECIMAL PLACES FOR OUTPUT FROM THE SUBROUTINE ALTB1D.
C
C-------------------------------------------------------------------------------
C
C              MODIFICATIONS MADE TO VERSION 1.0
C
C      REV: 20/11/89 BY J.J.SMITH IAEA ATOMIC AND MOLECULAR DATA UNIT
C            FOR EACH EVALUATION FUNCTION A SEPERATE SUBROUTINE IS
C            DEFINED TO HANDLE THE INTERACTIVE DIALOG.
C
C#######################################################################
C
      SUBROUTINE ALEV
C
C-------------------------------------------------------------------
C
C          EV = EVALUATE DATA VALUE FOR CURRENT ENTRY
C
C-------------------------------------------------------------------
C
C          IN ORDER TO ADD A NEW DATA TYPE TO THE "EV" COMMAND, ADD AN
C          ADDITIONAL "ELSE IF(...)" STATEMENT TO THE MAIN IF...THEN...
C          ELSE CHAIN THAT CHECKS THE "#..." EVALUATION LABEL. AN
C          SUBROUTINE HAS TO BE INTRODUCED TO HANDLE THE INTERACTIVE
C          DIALOG.  FOR EACH DATA TYPE A
C          CALL MUST BE INCLUDED TO READ THE REQUIRED PERCENTAGE
C          ACCURACY AND LIMITS FOR GENERATION OF THE LINEAR DATA POINT
C          SERIES. THE SUBROUTINE ALTB1D REQUIRES THE NUMBER OF
C          INTERVALS AND THE EXTENTS OF EACH INTERVAL. FOR FITTED DATA
C          FORMS SPECIFIED OVER ONE INTERVAL THE EXTENTS CAN BE PASSED
C          DIRECTLY FROM THE VALUES FROM THE CALL TO ALREVD. FOR FITTED
C          DATA FORMS WHICH ARE DEFINED OVER MULTIPLE INTERVALS,
C          ADDITIONAL STATEMENTS HAVE TO BE ADDED TO CONSTRUCT THE
C          ARRAY OF INTERVALS FOR THE SUBROUTINE ALTB1D. A MAXIMUM OF
C          TEN INTERVALS FOR FITTED DATA FORMS HAS BEEN DEFINED.
C          (SEE RELEVANT CODING AND FURTHER COMMENTS BELOW).
C
C          SPECIAL NOTE: IF THE 1-DIMENSIONAL DATA EVALUATION LOOP
C          ROUTINE ALTB1D(AL----) IS TO BE USED, THE NEW EVALUATION
C          SUBROUTINE NAME **MUST** BE ADDED TO THE EXTERNAL STATEMENT
C          AT THE BEGINNING OF ALEV, AND THIS EVALUATION ROUTINE
C          **MUST** BE IN STANDARD FORM !!!
C          (AGAIN, SEE RELEVANT CODING BELOW FOR MORE DETAIL AND
C          EXAMPLE)
C
C-------------------------------------------------------------------
C
      INCLUDE 'ALPCOM.FOR'
C
      INCLUDE 'ALCOM.FOR'
C
C
C-------------------------------------------------------------------
C
      IF(LESEQN .EQ. 0) THEN
        WRITE(LUTOUT, '(/'' NO ENTRY CURRENTLY READ IN'')' )
        RETURN
      ENDIF
C
C-------------------------------------------------------------------
C
C      BEGIN SUBROUTINE SWITCH TO SELECT CORRECT SUBROUTINE TO CALL
C      IN ORDER TO EVALUATE THIS ENTRY, BASED ON THE EVALUATION
C      FUNCTION LABEL FROM THE ENTRY.  THIS LABEL IS STORED AS THE
C      LAST BOOLEAN LABEL.
C
C-------------------------------------------------------------------
C
      IF(BL(NBL) .EQ. '#MEWE') THEN
        CALL RMEWE
C------
      ELSE IF(BL(NBL) .EQ. '#BELI') THEN
        CALL RBELI
C------
      ELSE IF(BL(NBL) .EQ. '#CHEB') THEN
        CALL RCHEB
C------
      ELSE IF(BL(NBL) .EQ. '#KINGEX') THEN
        CALL RKINGEX
C------
      ELSE IF(BL(NBL) .EQ. '#PHACX') THEN
        CALL RPHACX
C------
      ELSE IF(BL(NBL) .EQ. '#JAN1') THEN
        CALL RALJAN1
C------
      ELSE IF(BL(NBL) .EQ. '#'.AND. BL(1) .EQ. '$TABSPY') THEN
        CALL RTABSPY
C------
C
C      OTHER EVALUATION LABELS GO HERE IN STRUCTURE OF THE FORM:
C
C      ELSE IF(BL(NBL) .EQ. '#{EVALUATION LABEL}') THEN
C             CALL {APPROPRIATE SUBROUTINE TO DO INTERACTIVE
C                    EVALUATION}
C
C      FOR 1-D DATA TYPES WITH SUBROUTINES IN PROPER FORM, ALTB1D CAN
C      BE USED TO CALL THE EVALUATION SUBROUTINE DIRECTLY BY NAME, WITH NO
C      NO SEPARATE ROUTINE NEEDED FOR ALEV TO INTERACT WITH THE
C      TERMINAL ETC.
C
C      IF ALTB1D CANNOT BE USED, THE TYPICAL STRUCTURE OF SUBROUTINE
C      CALLED BY ALEV FOR INTERACTIVE EVALUATION OF A DATA ENTRY CON
C      CONTAINS:
C
C          {CODING TO CONVERT COEFFICIENT FIELD AS REQUIRED}
C          {CODING TO OBTAIN REQUIRED INDEPENDENT VARIABLE(S) FOR THIS
C           DATA TYPE}
C          {CODING TO CALL EVALUATION FUNCTION SUBROUTINE}
C          {CODING TO OUTPUT RESULTING DATA POINT VALUE(S)}
C          {LOOP BACK FOR MORE EVALUATIONS OR RETURN}
C
C------
C
C          CATCH-ALL FOR ENTRIES WITH EVALUATION LABELS NOT CURRENTLY IN
C          THIS VERSION OF ALADDIN
C
      ELSE
        WRITE(LUTOUT,
     &        '(/'' THIS EVALUATION LABEL CURRENTLY UNAVAILABLE'')')
C
C------
C
      ENDIF
C
      RETURN
      END
      SUBROUTINE ALTB1D(PERUNC,EVSUB,NRANGE,XRANGE)
C
C     CONSTRUCT A LINEARLY INTERPOLABLE TABLE OF VALUES OF A FUNCTION
C     OF A SINGLE CONTINUOUS VARIABLE.
C
C     THE ROUTINE WILL CONSTRUCT A TABLE OF VALUES (X,Y) WHERE THE
C     TABULATED (X,Y) PAIRS WILL BE EXACTLY REPRODUCE THE FUNCTION
C     AND LINEAR INTERPOLATION BETWEEN THE TABULATED PAIRS WILL BE
C     WITHIN THE FRACTIONAL UNCERTAINTY (PERUNC).  THE FUNCTION CAN
C     BE DEFINED OVER A NUMBER OF X INTERVALS.
C
C     THE PASSED VARIABLE "EVSUB" IS THE NAME OF THE SUBROUTINE TO
C     TO EVALUATE THE DATA POINT.  THIS SUBROUTINE MUST HAVE AN
C     ARGUMENT LIST IN STANDARD FORM, AND USE TYPE REAL DOUBLE PRECISION
C     COEFFICIENTS, FOR ALTB1D TO BE USED.
C     OTHERWISE, A SEPARATE ROUTINE FOR THAT DATA TYPE WILL HAVE TO BE
C     TO BE WRITTEN.
C
C     THE PARAMETERS REQUIRED BY SUBROUTINE ALTB1D ARE:
C
C       PERCUNC = FRACTIONAL PERCENTAGE UNCERTAINTY FOR GENERATION OF
C                 THE TABLE.
C       EVSUB = NAME OF EVALUATION SUBROUTINE
C              THE REQUIRED FORM FOR AN EVALUATION SUBROUTINE CALLABLE
C              USING ALTB1D IS:
C
C           SUBROUTINE {EVSUB} (PX, PCF, KNCF, PY, KERMSG)
C
C             WHERE PX = INDEPENDENT VARIABLE
C             PCF = ARRAY OF REAL DOUBLE PRECISION COEFFICIENTS
C                   (E.G., ALEVID USES ALRECF)
C             KNCF = NUMBER OF COEFFICIENTS IN PCF
C             PY = RETURNED DOUBLE PRECISION VALUE CORRESPONDING TO PX
C             KERMSG = CHARACTER ARRAY RETURNING ERROR MESSAGE (BLANK =
C                      NO ERROR)
C       NRANGE  = NUMBER OF RANGES/INTERVALS USED IN REPRESENTING THE
C                 FUNCBTION
C       XRANGE  = ARRAY OF RANGES USED FOR FUNCTION WHERE THE RANGES ARE
C                 DEFINED AS
C
C                         XRANGE(1) TO XRANGE(2)   WITH NRANGE = 1
C                     AND XRANGE(2) TO XRANGE(3)   WITH NRANGE = 2 ETC
C
C      ON SUCCESFUL COMPLETION OF THE GENERATION OF THE LINEARLY
C      INTERPOLABLE TABLE THE NUMBER OF (X,Y) PAIRS REPRESENTING THE
C      FUNCTION IS OUTPUT. OTHERWISE EITHER CONVERGENCE HAS FAILED DUE
C      TO THE FUNCTIONAL  FORM OF THE EXPRESSION OR BECAUSE MORE THAN
C      THE MAXIMUM NUMBER OF POINTS IN A TABLE HAS BEEN EXCEDDED.
C
      INCLUDE 'ALCOM.FOR'
C
      CHARACTER*4 KSIGN
      DOUBLE PRECISION XRANGE,XNOW,YNOW,XLAST,YLAST,XEND,
     1 YEND,XMID,YMID,YLIN,XFINAL,YFINAL,PERUNC,XTAB,YTAB,
     2 DX,DXSTEP,RATLST,RATNOW,XNOWL,DYMID,XNOW1,XNOW2,ALRND6
C
      DOUBLE PRECISION ZERO,HUNDRD,HALF,TEN,SMALLR,SMALLX
C
      DIMENSION XRANGE(*)
      DIMENSION XTAB(5000),YTAB(5000),XNORM(6),KSIGN(6),KEXP(6),
     1 MARKER(5000)
C
C------DEFINE DIMENSION OF TABLES (TO PREVENT MEMORY OVERFLOW)
C
      DATA IPTMAX/5000/
      DATA MAXLP1/20/
      DATA MAXLP2/40/
      DATA ZERO/0.0D+00/
      DATA HALF/5.0D-01/
      DATA TEN/1.0D+01/
      DATA HUNDRD/1.0D+02/
      DATA SMALLR/9.99D-01/
      DATA SMALLX/1.0D-07/
C
C------INITIALIZE TOTAL POINT COUNT (ALL X REGIONS COMBINED).
C
      KPT=0
C
C          TREAT EACH X RANGE OF FUNCTION SEPARATELY.
C
      DO 340 IRANGE=1,NRANGE
C
C------- DETERMINE FUNCTION AT START AN END OF INTERVAL
C
      XNOW=XRANGE(IRANGE)
      IF (IRANGE .NE. 1) XNOW = XNOW + (SMALLX * XNOW)
      CALL EVSUB (XNOW, CF, NCF, YNOW, ERRMSG)
      IF (ERRMSG .NE. ' ')       GO TO 370
      XEND=XRANGE(IRANGE+1)
      CALL EVSUB (XEND, CF, NCF, YEND, ERRMSG)
      IF (ERRMSG .NE. ' ')       GO TO 370
C
C          SELECT LINEAR OR LOG INTERVALS.
C
      IF(XNOW.LE.ZERO) GO TO 10
      IF(XEND.LT.TEN*XNOW) GO TO 10
C
C------LOG INTERVALS.
C
      MYWAY=2
      DXSTEP=DLOG(XEND/XNOW)/HUNDRD
      DXSTEP=DEXP(DXSTEP)
      XNOW1=ALRND6(XNOW*DXSTEP)
      GO TO 20
C
C------LINEAR INTERVALS.
C
   10 MYWAY=1
      DXSTEP=(XEND-XNOW)/HUNDRD
      XNOW1=ALRND6(XNOW+DXSTEP)
C
C         START CONSTRUCTING TABLE.
C
C------INITIALIZE TABLE TO CONTAIN START OF INTERVAL.
C
   20 IPT=1
      XTAB(IPT)=XNOW
      YTAB(IPT)=YNOW
      MARKER(IPT)=0
C
C------NO LINEARIZING IF ALLOWABLE ERROR IS NOT POSITIVE OR CURRENT
C------INTERVAL IS OF ZERO LENGTH.
C
      IF(PERUNC.GT.ZERO.AND.XEND.GT.XNOW) GO TO 30
      IPT=2
      XTAB(IPT)=XEND
      YTAB(IPT)=YEND
      MARKER(IPT)=1
      GO TO 240
C
C          START NEXT X INTERVAL.
C
C------DEFINE BEGINNING OF SUB-INTERVAL.
C
   30 XLAST=XNOW
      YLAST=YNOW
C
C------INITIALIZE TABLE TO CONTAIN START OF INTERVAL.
C
      IF(XNOW.LT.XNOW1) GO TO 80
C
C------SAVE UNROUNDED VALUE TO SKIP VERY SHORT INTERVALS.
C
   40 XNOWL=XNOW1
C
C------DEFINE END OF SUB-INTERVAL.
C
      IF(MYWAY.EQ.2) GO TO 50
C
C------LINEAR INTERVALS.
C
      XNOW1=ALRND6(XNOW1+DXSTEP)
      XNOW2=XNOW1+DXSTEP
C
C------SPECIAL TREATMENT FOR 0.0.
C
      IF(XLAST.LT.ZERO.AND.XNOW1.GT.ZERO) XNOW1=ZERO
      GO TO 60
C
C------LOG INTERVALS.
C
   50 XNOW1=ALRND6(XNOW1*DXSTEP)
      XNOW2=XNOW1*DXSTEP
C
C------AVOID ROUNDOFF NEAR END OF INTERVAL.
C
   60 IF(XNOW2.GE.XEND) GO TO 90
C
C------STOP AT END OF INTERVAL.
C
      IF(XNOW1.GE.XEND) GO TO 90
C
C------IGNORE INTERVALS THAT CANNOT BE SUB-DIVIDED AT LEAST ONCE.
C
      XMID=HALF*(XNOW1+XNOW)
      XMID=ALRND6(XMID)
      IF(XMID.LE.XNOW.OR.XMID.GE.XNOW1) GO TO 70
      GO TO 80
   70 XNOW1=XNOWL+DXSTEP
      GO TO 40
   80 XNOW=XNOW1
C
C------INITIALIZE COUNT OF NUMBER OF TIMES INTERVAL HAS BEEN SUB-DIVIDED
C
      LOOP=0
      IF(XNOW.LT.XEND) GO TO 100
   90 XNOW=XEND
      YNOW=YEND
      XNOW1=XEND
      GO TO 110
  100 CALL EVSUB (XNOW, CF, NCF, YNOW, ERRMSG)
      IF (ERRMSG .NE. ' ')       GO TO 370
C
C          SUB-DIVIDE INTERVAL AND CHECK FOR CONVERGENCE.
C
C------INCREMENT COUNT OF NUMBER OF TIMES INTERVAL HAS BEEN SUB-DIVIDED
C------AND CHECK FOR MAXIMUM NUMBER OF SUB-DIVISIONS.
C
  110 LOOP=LOOP+1
      IF(LOOP.GT.MAXLP2) GO TO 150
C
C------COMPARE EXACT FUNCTION TO LINEAR INTERPOLATION OVER CURRENT
C------SUBINTERVAL.
C
      XMID=HALF*(XLAST+XNOW)
      XMID = ALRND6(XMID)
C
C------CHECK FOR MINIMUM X INTERVAL.
C
      IF(XMID.LE.XLAST.OR.XMID.GE.XNOW) GO TO 150
      CALL EVSUB (XMID, CF, NCF, YMID, ERRMSG)
      IF (ERRMSG .NE. ' ')  GO TO 370
      YLIN=((XMID-XLAST)*YNOW+(XNOW-XMID)*YLAST)/(XNOW-XLAST)
C
C------CHECK FOR CONVERGENCE.
C
      DYMID=YMID-YLIN
      IF(DABS(DYMID).LE.DABS(PERUNC*YMID)) GO TO 140
C
C------IF POSSIBLE SAVE CURRENT AND LAST RATIO AND CHECK TO SEE IF
C------METHOD IS CONVERGING.
C
      IF(YMID.EQ.ZERO) GO TO 130
      RATNOW=DABS((YMID-YLIN)/(PERUNC*YMID))
      IF(LOOP.LT.MAXLP1) GO TO 120
C
C------STOP SUB-DIVIDING IF METHOD IS NOT CONVERGING.
C
      IF(RATNOW.GE.SMALLR*RATLST) GO TO 150
  120 RATLST=RATNOW
C
C------NO CONVERERENCE...HALF SUBINTERVAL AND CONTINUE.
C
  130 XNOW=XMID
      YNOW=YMID
      GO TO 110
C
C------CONVERGENCE...SAVE MIDPOINT AND ENDPOINT FOR LATER THINNING.
C
  140 IF(IPT.GE.IPTMAX) GO TO 170
      IPT=IPT+1
      XTAB(IPT)=XMID
      YTAB(IPT)=YMID
      MARKER(IPT)=0
      IF(IPT.GE.IPTMAX) GO TO 170
      IPT=IPT+1
      XTAB(IPT)=XNOW
      YTAB(IPT)=YNOW
      MARKER(IPT)=0
      GO TO 160
C
C------MINIMUM X INTERVAL OR MAXIMUM NUMBER OF SUB-DIVISIONS...SAVE
C------ENDPOINT AND MARK AS UNCONVERGED.
C
  150 IF(IPT.GE.IPTMAX) GO TO 170
      IPT=IPT+1
      XTAB(IPT)=XNOW
      YTAB(IPT)=YNOW
      MARKER(IPT)=1
C
C------CONTINUE UNTIL END OF X INTERVAL.
C
  160 IF(XNOW.LT.XEND) GO TO 30
      GO TO 180
C
C          OVER 5000 POINTS REQUIRED. SAVE END OF X RANGE AND INDICATE N
C          CONVERGENCE...TOO MANY POINTS REQUIRED.
C
  170 XTAB(IPT)=XEND
      YTAB(IPT)=YEND
      MARKER(IPT)=2
C
C          THIN RESULTS.
C
C------INITIALIZE ORIGINAL POINT COUNT, FINAL POINT COUNT AND POINTE
C------BEGINNING OF THINNING INTERVAL.
C
  180 NPT=IPT
      IPT=1
      INOW=1
C
C------DEFINE X AND Y AT BEGINNING OF THINNING INTERVAL.
C
  190 XLAST=XTAB(INOW)
      YLAST=YTAB(INOW)
C
C------INITIALIZE THINNING INDICES TO TRY TO THIN 1 POINT ABOVE INOW
C
      JMIN=INOW+1
      JMAX=INOW+1
C
C------FINISHED IF NO MORE POINTS TO THIN.
C
  200 IF(JMAX.GE.NPT) GO TO 230
C
C------SAVE POINT IF NO CONVERGENCE OR X = 0.0.
C
      IF(MARKER(JMAX).NE.0.OR.XTAB(JMAX).EQ.ZERO) GO TO 220
C
C------DEFINE INDEX TO END OF THINNING INTERVAL.
C
      ITOP=JMAX+1
C
C------DEFINE X AND Y AT END OF THINNING INTERVAL AND X INTERVAL.
C
      XNOW=XTAB(ITOP)
      YNOW=YTAB(ITOP)
      DX=XNOW-XLAST
C
C------SEE IF POINTS JMIN TO JMAX CAN BE APPROXIMATED BY LINEAR
C------INTERPOLATION BETWEEN ENDS OF THINNING INTERVAL.
C
      DO 210 J=JMIN,JMAX
      XMID=XTAB(J)
      YMID=YTAB(J)
      YLIN=((XMID-XLAST)*YNOW+(XNOW-XMID)*YLAST)/DX
      IF(DABS(YMID-YLIN).GE.DABS(PERUNC*YMID)) GO TO 220
  210 CONTINUE
C
C------CAN ELIMINATE ALL POINTS JMIN TO JMAX. INCREMENT JMAX AND TRY
C------ELIMINATE MORE POINTS.
C
      JMAX=JMAX+1
      GO TO 200
C
C------MUST KEEP POINT JMAX.
C
  220 IPT=IPT+1
      XTAB(IPT)=XTAB(JMAX)
      YTAB(IPT)=YTAB(JMAX)
      MARKER(IPT)=MARKER(JMAX)
C
C------RESET INDEX TO RESTART THINNING AT JMAX
C
      INOW=JMAX
      GO TO 190
C
C------END OF THINNING SAVE LAST POINT.
C
  230 IPT=IPT+1
      XTAB(IPT)=XTAB(NPT)
      YTAB(IPT)=YTAB(NPT)
      MARKER(IPT)=MARKER(NPT)
C
C          PRINT RESULTS AFTER THINNING
C
C------IGNORE FIRST POINT IF AT BOUNDARY BETWEEN X RANGES AND FUNCTION
C------IS CONTINUOUS.
C
  240 NPT=1
      IF(IRANGE.EQ.1) GO TO 250
      IF(XTAB(1).EQ.XFINAL.AND.YTAB(1).EQ.YFINAL) NPT=2
C
C------OUTPUT TABULATED VALUES
C
  250 DO 330 I=NPT,IPT
C
C------INCREMENT TOTAL POINT COUNT AND CONVERT TO OUTPUT FORM.
C
      KPT=KPT+1
C
C     IF MARKER FOR THE INTERVAL LINEARISATION HAS FAILED A MESSAGE
C     IS DISPLAYED AMD LINEARISATION TERMINATED
C
      IF(MARKER(I).NE.0) GO TO 350
C
C------CONVERGENCE - WRITE DATA TO LUEF OR LUTOUT
C
      IF(EFNAME .NE. '**NONE**') THEN
        WRITE(LUEF,'(2(1PE11.4))') XTAB(I),YTAB(I)
      ELSE
        WRITE(LUTOUT,'(2(1PE11.4))') XTAB(I),YTAB(I)
      ENDIF
  330 CONTINUE
C
C------SAVE VALUE FROM END OF X REGION (TO CHECK FOR DISCONTINUITY).
C
      XFINAL=XEND
      YFINAL=YEND
C
C------END OF RANGE.
C
  340 CONTINUE
C
C------LINEARIZING IS COMPLETED.
C
      WRITE(LUTOUT, '('' LINEAR DATA POINT SERIES  '', I4,
     &                 '' VALUES IN TABLE'')') KPT
      IF(EFNAME .NE. '**NONE**')
     &  WRITE(LUEF, '(25X,'' LINEAR DATA POINT SERIES  '', I4,
     &                 '' VALUES IN TABLE'')') KPT
      RETURN
C
  350 IF(MARKER(I).EQ.2) GO TO 360
C
C------NO CONVERGENCE.
C
      WRITE(LUTOUT,'(I6,(2(1PE12.5)),'' NO CONVERGENCE OF ALTB1D'')')
     &      KPT,XTAB(J),YTAB(J)
C
      IF(EFNAME .NE. '**NONE**')
     &  WRITE(LUEF,'(I6,(2(1PE12.5)),'' NO CONVERGENCE OF ALTB1D'')')
     &      KPT,XTAB(J),YTAB(J)
      GO TO 380
C
C------TOO MANY POINTS.
C
  360 WRITE(LUTOUT,'(I6,2(1PE12.5),
     &  '' CONVERGENCE FAILED TOO MANY POINTS'')') KPT,XTAB(J),YTAB(J)
C
      IF(EFNAME .NE. '**NONE**')
     &  WRITE(LUEF,'(I6,2(1PE12.5),
     &  '' CONVERGENCE FAILED TOO MANY POINTS'')') KPT,XTAB(J),YTAB(J)
      GO TO 380
C
C------ PRINT FINAL ERROR MESSAGE ON LUTOUT
C
  370 WRITE(LUTOUT,'('' LINEARISING TERMINATED : '',A)' ) ERRMSG
  380 RETURN
      END
C
C####################################################################
C
      SUBROUTINE ALREVD(PERUNC,ZX1,ZX2)
C
      INCLUDE 'ALCOM.FOR'
C
      DOUBLE PRECISION PERUNC, ZX1, ZX2
C
      PERUNC=0.0D0
C
C      READ DATA FOR FUNCTION EVALUATION
C
  100 WRITE(LUTOUT,
     &    '(/'' EV OUTPUT DATA CURRENTLY DIRECTED TO FILE: '', A,
     &    /'' (TO CHANGE OUTPUT FILE, USE EF COMMAND AFTER'',
     &         '' EXITING AT NEXT PROMPT)'')') EFNAME
C
C      DETERMINE VALUES OF INDEPENDENT VARIABLE TO BE EVALUATED
C
      WRITE(LUTOUT,
     &       '(/''$LINEAR DATA POINT SERIES, INPUT  '')' )
C
      WRITE(LUTOUT,
     &'(/''$% UNCERTAINTY, FIRST, LAST VALUE (,,, TO EXIT) ? '')' )
      READ(LUTIN, *, IOSTAT=IOS) PERUNC, ZX1, ZX2
      IF (PERUNC .EQ. 0.0D0) RETURN
      IF(IOS .NE. 0 .OR. PERUNC .LT. 0.0D0
     &        .OR. (PERUNC .GT. 0.0D0 .AND. ZX2 .LE. ZX1) ) THEN
        WRITE(LUTOUT, '('' INVALID INPUT'')')
        GO TO 100
      ENDIF
      PERUNC = PERUNC*0.01D0
C
C          BEFORE DATA OUTPUT, WRITE APPROPRIATE HEADER LINE
C
      IF(EFNAME .NE. '**NONE**') THEN
        WRITE(LUTOUT, '('' ENTER HEADER LINE FOR THIS DATA ?'')')
        READ(LUTIN, '(A)') CMND
        WRITE(LUTOUT, '(1X,''DATA IS BEING WRITTEN TO   '',A)') EFNAME
        WRITE(LUEF, '(A)') CMND
      ENDIF
      RETURN
      END
C
C####################################################################
C
      DOUBLE PRECISION FUNCTION ALRND6(X)
C
C          ROUND NUMBER TO 6 DIGITS ACCURACY.
C
      DOUBLE PRECISION X,ABSX,XR,OFFSET,ZERO,HALF,ONE,TEN
      DATA ZERO/0.0D+00/
      DATA HALF/5.0D-01/
      DATA ONE/1.0D+00/
      DATA TEN/1.0D+01/
      DATA NMAX/1000000/
C
C------AVOID TRYING TO TAKE THE LOG OF 0.0 OR NEGATIVE NUMBERS.
C
      ABSX=DABS(X)
      IF(ABSX.GT.ZERO) GO TO 10
      ALRND6=ZERO
      RETURN
C
C------DEFINE EXPONENT TO ROUND TO NORMAL FORM (1.00.. TO 9.99....).
C
   10 I=DLOG10(ABSX)
      IF(ABSX.LT.ONE) I=I-1
C
C------ROUND TO 6 DIGIT INTEGER (TREAT -, 0 AND + EXPONENTS SEPARATE
C------TO AVOID ROUNDOFF WITH NEGATIVE POWERS OF 10..E.G., 10E-1 IS
C------A BINARY MULTIPLE).
C
      I=I-5
      IF(I) 20,40,60
C
C          NEGATIVE EXPONENT.
C
C------ROUND TO 6 DIGIT INTEGER.
C
   20 OFFSET=TEN**(-I)
      N=ABSX*OFFSET+HALF
C
C------CHECK FOR CHANGE IN EXPONENT DUE TO ROUNDING.
C
      IF(N.LT.NMAX) GO TO 30
      OFFSET=OFFSET/TEN
      N=ABSX*OFFSET+HALF
C
C------CONVERT TO 6 DIGIT ACCURACY FLOATING POINT NUMBER AND
C------RE-NORMALISE
C
   30 XR=N
      XR=XR/OFFSET
      GO TO 80
C
C          ZERO EXPONENT.
C
C------ROUND TO 6 DIGIT INTEGER.
C
   40 OFFSET=ONE
      N=ABSX+HALF
C
C------CHECK FOR CHANGE IN EXPONENT DUE TO ROUNDING.
C
      IF(N.LT.NMAX) GO TO 50
      I=-1
      OFFSET=TEN
      N=ABSX/OFFSET+HALF
C
C------CONVERT TO 6 DIGIT ACCURACY FLOATING POINT NUMBER AND
C------RE-NORMALISE
C
   50 XR=N
      XR=XR*OFFSET
      GO TO 80
C
C          POSITIVE EXPONENT.
C
C------ROUND TO 6 DIGIT INTEGER.
C
   60 OFFSET=TEN**I
      N=ABSX/OFFSET+HALF
C
C------CHECK FOR CHANGE IN EXPONENT DUE TO ROUNDING.
C
      IF(N.LT.NMAX) GO TO 70
      OFFSET=TEN*OFFSET
      N=ABSX/OFFSET+HALF
C
C------CONVERT TO 6 DIGIT ACCURACY FLOATING POINT NUMBER AND
C------RE-NORMALISE
C
   70 XR=N
      XR=XR*OFFSET
C
C------DEFINE SIGN OF RESULT.
C
   80 IF(X.LT.ZERO) XR=-XR
      ALRND6=XR
      RETURN
      END
C
C
C-----------------------------------------------------------------
C
C     ALADDIN SUBROUTINES TO SUPPORT VERSION 1.0
C
C-----------------------------------------------------------------
C
C
C####################################################################
C
      SUBROUTINE RMEWE
C
      EXTERNAL ALMEWE
C
      INCLUDE 'ALPCOM.FOR'
C
      INCLUDE 'ALCOM.FOR'
C
      DOUBLE PRECISION PERUNC, XRANGE
      DIMENSION XRANGE(10)
C
        WRITE(LUTOUT,
     &  '(/'' MEWE SPECTRAL LINE EXCITATION RATE COEFFICIENT (cm3/s),'',
     &    /'' AS FUNCTION OF ELECTRON TEMPERATURE (keV)'')')
C
C--------------  RETREIVE NUMERIC COEFFICIENTS
C
        CALL ALRECF(CF, NCF, NCFMX, ERRMSG)
        IF(NCF .NE. 8 .AND. ERRMSG .EQ. ' ')
     &         ERRMSG = 'INCORRECT NUMBER OF COEFFICIENTS'
        IF (ERRMSG .NE.  ' ') THEN
          WRITE(LUTOUT, '(/1X,A)') ERRMSG
          RETURN
        ENDIF
C
        CALL ALREVD(PERUNC,XRANGE(1),XRANGE(2))
        IF (PERUNC .EQ. 0.0D0) RETURN
        NRANGE=1
C
        CALL ALTB1D(PERUNC,ALMEWE,NRANGE,XRANGE)
C
        RETURN
        END
C
C####################################################################
C
      SUBROUTINE RBELI
C
      EXTERNAL ALBELI
C
      INCLUDE 'ALPCOM.FOR'
C
      INCLUDE 'ALCOM.FOR'
C
      DOUBLE PRECISION PERUNC, XRANGE
      DOUBLE PRECISION ZX1, ZX2
      DIMENSION XRANGE(10)
C
        WRITE(LUTOUT,
     &       '(/'' EVALUATION FUNCTION ALBELI'')')
        WRITE(LUTOUT,
     &       '(/'' ELECTRON IMPACT IONIZATION CROSS SECTION cm[2], '',
     &         /'' AS FUNCTION OF ELECTRON ENERGY (eV)'')')
C
C--------------  RETREIVE NUMERIC COEFFICIENTS
C
        CALL ALRECF(CF, NCF, NCFMX, ERRMSG)
        IF (NCF. EQ. 0) THEN
          WRITE(LUTOUT,
     &     '(/'' NO FITTING PARAMETERS   EV FUNCTION TERMINATED'')' )
          RETURN
        ENDIF
        IF(ERRMSG .EQ. ' ') THEN
          WRITE(LUTOUT, '(/'' (INPUT VALID RANGE FROM THRESHOLD OF'',
     &          1PE12.4, '' (eV) )'')') CF(1)
        ELSE
          WRITE(LUTOUT, '(/1X,A)') ERRMSG
          RETURN
        ENDIF
C
C--------------  READ PERCENTAGE UNCERTAINTY AND RANGE FOR EVALUATION
C
        CALL ALREVD(PERUNC,ZX1,ZX2)
        IF (PERUNC .EQ. 0.0D0) RETURN
C
C------      DEFINE EXTENT AND NUMBER OF INTERVALS FOR TABLE
C------ (SEE ROUTINE ALBELI FOR DETAILS OF INTERPRETATION OF PARAMETERS)
C
        XRANGE(1) = ZX1
        IF (NCF .GT. 7 ) THEN
          IF (ZX2 .LT. CF(8) .OR. ZX1 .GT. CF(8)) THEN
C
C------      RANGE ZX1 YTO ZX2 OUTSIDE OF AUTOIONIZATION THRESHOLD
C
            XRANGE(2) = ZX2
            NRANGE=1
          ELSE
C
C------      RANGE ZX1 YTO ZX2 BOUNDS THE AUTOIONIZATION THRESHOLD
C
            XRANGE(2) = CF(8)
            XRANGE(3) = ZX2
            NRANGE=2
          ENDIF
        ELSE
C
C------      NO AUTOIONIZATION STRUCTURE
C
          XRANGE(2) = ZX2
          NRANGE = 1
        ENDIF
C
C--------------  GENERATE TABLE OF VALUES OF FUNCTION
C
        CALL ALTB1D(PERUNC,ALBELI,NRANGE,XRANGE)
     
        RETURN
        END
C
C####################################################################
C
      SUBROUTINE RCHEB
C
      EXTERNAL ALCHEB
C
      INCLUDE 'ALPCOM.FOR'
C
      INCLUDE 'ALCOM.FOR'
C
      DOUBLE PRECISION PERUNC, XRANGE
      DIMENSION XRANGE(10)
C
        WRITE(LUTOUT,
     &       '(/'' ORNL:CFADC CHEBYSHEV POLYNOMIAL FIT'')' )
C
        CALL ALRECF(CF, NCF, NCFMX, ERRMSG)
        IF(NCF .NE. 11 .AND. ERRMSG .EQ. ' ')
     &         ERRMSG = 'INCORRECT NUMBER OF COEFFICIENTS'
        IF(ERRMSG .EQ. ' ') THEN
          WRITE(LUTOUT, '(/'' (INPUT VALID RANGE FROM'',
     &          1PE10.2, '' TO'', 1PE10.2, '')'')' ) CF(10), CF(11)
        ELSE
          WRITE(LUTOUT, '(/1X,A)') ERRMSG
          RETURN
        ENDIF
C
        CALL ALREVD(PERUNC,XRANGE(1),XRANGE(2))
        IF (PERUNC .EQ. 0.0D0) RETURN
        NRANGE=1
C
        CALL ALTB1D(PERUNC,ALCHEB,NRANGE,XRANGE)
C
        RETURN
        END
C
C####################################################################
C
      SUBROUTINE RKINGEX
C
      EXTERNAL ALKING
C
      INCLUDE 'ALPCOM.FOR'
C
      INCLUDE 'ALCOM.FOR'
C
      DOUBLE PRECISION PERUNC, XRANGE
C
      DIMENSION XRANGE(10)
        WRITE(LUTOUT,
     &       '(/'' EVALUATION FUNCTION ALKING'')' )
        WRITE(LUTOUT,
     &  '(/'' POLYNOMIAL FIT OF EXCITATION RATE COEFFICIENT (cm3/s),'',
     &    /'' AS FUNCTION OF ELECTRON TEMPERATURE (eV)'')')
C
        CALL ALRECF(CF, NCF, NCFMX, ERRMSG)
        IF(ERRMSG .EQ. ' ') THEN
          WRITE(LUTOUT, '(/'' (INPUT VALID RANGE FROM'',
     &          1PE12.4, '' TO'', 1PE12.4, '')'')' ) CF(2), CF(3)
        ELSE
          WRITE(LUTOUT, '(/1X,A)') ERRMSG
          RETURN
        ENDIF
C
        CALL ALREVD(PERUNC,XRANGE(1),XRANGE(2))
        IF (PERUNC .EQ. 0.0D0) RETURN
        NRANGE=1
C
        CALL ALTB1D(PERUNC,ALKING,NRANGE,XRANGE)
C
        RETURN
        END
C
C####################################################################
C
      SUBROUTINE RPHACX
C
      EXTERNAL ALPHACX
C
      INCLUDE 'ALPCOM.FOR'
C
      INCLUDE 'ALCOM.FOR'
C
      DOUBLE PRECISION PERUNC, XRANGE
      DOUBLE PRECISION Q
      DIMENSION XRANGE(10)
C
        WRITE(LUTOUT,
     &       '(/'' EVALUATION FUNCTION ALPHACX'')' )
        WRITE(LUTOUT,
     &       '(/'' CHARGE EXCHANGE CROSS SECTION (cm[2]),'',
     &         /'' AS FUNCTION OF PROJECTILE ENERGY (eV/amu)'')')
C
        CALL ALRECF(CF, NCF, NCFMX, ERRMSG)
        IF(ERRMSG .NE. ' ') THEN
          WRITE(LUTOUT, '(/1X,A)') ERRMSG
          RETURN
        ENDIF
C
C-------   READ THE CHARGE STATE OF THE ION AND THEN DETERMINE THE
C-------   ENERGY RANGE OF VALIDITY FROM THE SCALED ENERGY
C
  10    WRITE(LUTOUT, '(/'' INPUT THE CHARGE STATE OF THE ION,'',
     &  '' (TO EXIT INPUT ,) ? '')' )
C
        INPTS = 0
        READ(LUTIN, *, IOSTAT=IOS) Q
        IF (Q .EQ. 0.0D0 ) RETURN
        IF (IOS .NE. 0) THEN
          WRITE(LUTOUT, '('' INVALID INPUT'')')
          GO TO 10
        ENDIF
C
C-------   THE SCALED CROSS SECTION APPLIES TO CHARGE STATES FROM
C-------   5 TO 26.
C
        IF(Q .LT. 5 . OR. Q .GT. 26) THEN
          WRITE(LUTOUT,
     &     '(/'' INVALID INPUT, THE CHARGE STATE (Q) MUST '',
     &       ''BE A VALUE IN THE RANGE 5 <= Q <= 26 '')')
          GO TO 10
        ENDIF
C
C-------   PASS Q IN THE COEFFICIENT ARRAY TO EVALUATION FUNCTION
C-------   ALPHACX. THERE ARE SIX COEFFICIENTS FOR THESE ENTRIES
C
       CF(7) = Q
C
       ELOW = CF(1)*DSQRT(Q)
       EHIGH = CF(2)*DSQRT(Q)
       WRITE(LUTOUT, '(/'' (INPUT VALID RANGE FROM'',
     &       1PE11.3, '' TO'', 1PE11.3, '' (eV/amu))'')' ) ELOW, EHIGH
C
      CALL ALREVD(PERUNC,XRANGE(1),XRANGE(2))
      IF (PERUNC .EQ. 0.0D0) RETURN
      NRANGE=1
C
      CALL ALTB1D(PERUNC,ALPHACX,NRANGE,XRANGE)
C
      RETURN
      END
C
C####################################################################
C
      SUBROUTINE RALJAN1
C
      EXTERNAL ALJAN1
C
      INCLUDE 'ALPCOM.FOR'
C
      INCLUDE 'ALCOM.FOR'
C
      DOUBLE PRECISION PERUNC, XRANGE
      DIMENSION XRANGE(10)
C
        WRITE(LUTOUT,
     &      '(/'' EVALUATION FUNCTION JAN1 - LOG(N) POLYNOMIAL FIT'')')
C
        CALL ALRECF(CF, NCF, NCFMX, ERRMSG)
        IF(NCF .NE. 11 .AND. ERRMSG .EQ. ' ')
     &         ERRMSG = 'INCORRECT NUMBER OF COEFFICIENTS'
        IF(ERRMSG .EQ. ' ') THEN
          WRITE(LUTOUT, '(/'' (INPUT VALID RANGE FROM'',
     &          1PE10.2, '' TO'', 1PE10.2, '')'')' ) CF(1), CF(2)
        ELSE
          WRITE(LUTOUT, '(/1X,A)') ERRMSG
          RETURN
        ENDIF
C
        CALL ALREVD(PERUNC,XRANGE(1),XRANGE(2))
        IF (PERUNC .EQ. 0.0D0) RETURN
        NRANGE=1
C
        CALL ALTB1D(PERUNC,ALJAN1,NRANGE,XRANGE)
C
        RETURN
        END
C
C####################################################################
C
      SUBROUTINE RTABSPY
C
      EXTERNAL ETABSQ
C
      INCLUDE 'ALPCOM.FOR'
C
      INCLUDE 'ALCOM.FOR'
C
      CHARACTER*2  CTARGT, CION
      DOUBLE PRECISION PERUNC, XRANGE
      DIMENSION XRANGE(10)
C
        WRITE(LUTOUT,
     &       '(/'' EVALUATION FUNCTION ETABSQ'')')
        WRITE(LUTOUT,
     &   '(/'' SPUTTERING YIELDS (ATOMS/ION) OF MONATOMIC SOLIDS),'',
     &     /'' AS FUNCTION OF INCIDENT ION ENERGY (eV)'')')
C
C------  THE INCIDENT ION AND TARGET PAIR MUST BE SPECIFIED
C
 100    WRITE(LUTOUT,
     &'(/'' ENTER THE INCIDENT ION, TARGET PAIR '',
     &   ''(WITH 2 CHARACTERS FOR EACH REACTANT)'',
     & /'' EXAMPLE H ,Cu  . (TO EXIT INPUT ,,) ? '')' )
C
      READ(LUTIN, '(A2,1X,A2)', IOSTAT=IOS) CION,CTARGT
      INPTS = 0
C
C------  CHECK RETURNED VALUES FOR THE INCIDENT ION AND TARGET
C
      IF (CION .EQ. ',,') RETURN
      IF (CION .EQ. ' ' .AND. CTARGT .EQ. ' ') RETURN
      IF(IOS .NE. 0 .OR. CION .EQ. ' ' .OR. CTARGT .EQ. ' ')  THEN
        WRITE(LUTOUT, '('' INVALID INPUT'')')
        GO TO 100
      ENDIF
C
C--  CALL RTABSQ TO READ DATA FOR CALCULATION OF SPUTTERING YIELDS
C--  CHECK THE TARGET, ION COMBINATION AND PASS BACK ARRAY CF FOR
C--  ROUTINE ETABSQ TO PRODUCE THE SPUTTERING YIELD.
C
        CALL RTABSQ(CION,CTARGT,CF,ERRMSG)
        IF(ERRMSG .NE. ' ') THEN
          WRITE(LUTOUT, '(/1X,A)') ERRMSG
          RETURN
        ENDIF
C-----
C-----  DETERMINE AND PRINT THE THRESHOLD ENERGY FOR THE REACTION
C-----  (ROUNDED UP TO THE NEXT .1 EV).
C-----  CF(2) = ATOMIC MASS OF ION   CF(4) = ATOMIC MASS OF TARGET
C-----  CF(5) = SUBLIMATION ENERGY.
C
        GAMMA= 4.0D0 * CF(2)* CF(4) / (CF(2) + CF(4))**2.0D0
        IF (CF(2) .GE. CF(4))    THEN
          ETH = (4.0D0 / 3.0D0)**6.0D0 * CF(5) / GAMMA
        ELSE
          ETH = (CF(5) / GAMMA) *
     &     ((2.0D0 * (CF(2) + CF(4))) / (CF(2) + 2.0D0 * CF(4)))**6.0D0
        ENDIF
        ETH = ETH + 0.0999D0
        WRITE(LUTOUT, '(/'' (INPUT VALID RANGE FROM THRESHOLD OF'',
     &          1PE10.2, '' (eV) )'')') ETH
C
        CALL ALREVD(PERUNC,XRANGE(1),XRANGE(2))
        IF (PERUNC .EQ. 0.0D0) RETURN
        NRANGE=1
C
        CALL ALTB1D(PERUNC,ETABSQ,NRANGE,XRANGE)
C
C
        RETURN
        END