      SUBROUTINE DINTA (ANSWER,WORK,IFLAG)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
c>> 2009-11/03 DINTA Krogh  Set initial value for a few variables.
c>> 2009-10-17 DINTA Krogh  AACUM now set with ACUM on coalescing nodes.
c>> 2008-02-23 DINTA Krogh/Snyder Misc. Changes
c>> 2001-09-11 DINTA Krogh  Fixed so really small intervals not picked.
c>> 1996-04-27 DINTA Krogh  Changes to use .C. and C%%.
C>> 1996-03-31 DINTA Krogh  Removed unused variable in common.
c>> 1996-03-30 DINTA Krogh  Change specific intrinsics to generics.
c>> 1995-11-28 DINTA Krogh  Converted from SFTRAN to Fortran 77.
C>> 1995-10-24 DINTA Krogh  Removed blanks in numbers for C conversion.
C>> 1994-11-14 DINTA Krogh  Declared all vars.
c>> 1994-10-19 DINTA Krogh  Changes to use M77CON
c>> 1994-09-01 DINTA Snyder don't give round-off msg if err << eps
C>> 1994-08-22 DINTA Krogh -- Modified data for C conversion.
c>> 1994-08-19 DINTA Snyder fix divide by 0 if ERR=0 and don't believe.
c>> 1994-08-15 DINTA Snyder convert MAX(A,B,C) to MAX(A,MAX(B,C))
c>> 1994-07-07 DINTA Snyder set up for CHGTYP.
C>> 1994-06-29 DINTA Snyder suppress discontinuity msg if one jump.
C>> 1994-06-22 DINTA Snyder don't use RE when reducing stepsize.
C>> 1993-05-26 DINTA Krogh -- Added data for C conversion.  Needs edit.
C>> 1993-05-18 DINTA Krogh -- Changed "END" to "END PROGRAM"
C>> 1992-03-03 DINTA Added call to DINTO for an error message.
C>> 1991-09-20 DINTA Krogh converted '(1)' dimensioning to '(*)'.
C>> 1989-11-08 DINTA Snyder  Revise coalesced abscissa test for CRAY
C>> 1988-05-17 DINTA Snyder  Initial code.
c--D replaces "?": ?INT, ?INTA, ?INTC, ?INTDL, ?INTDU, ?INTEC
c--&        ?INTF, ?INTM, ?INTMA, ?INTNC, ?INTNS, ?INTO, ?INTSM, ?INT1
C
C     THIS IS THE INTEGRATION PROGRAM FOR DINT.
C     THE RESULT IS OBTAINED USING QUADRATURE FORMULAE DUE TO
C     T. N. L. PATTERSON, MATHEMATICS OF COMPUTATION, VOLUME 22,
C     PAGES 847-856, 1968.
C
C     *****    WARNING   ***********************************************
C
C     THE RELIABILITY AND EFFICIENCY OF THIS PROGRAM ARE STRONGLY
C     INFLUENCED BY DISCONTINUITIES IN THE FUNCTION OR IT'S
C     DERIVATIVES, INTEGRABLE SINGULARITIES ON THE REAL AXIS,
C     NON-INTEGRABLE SINGULARITIES NEAR THE REGION, AND POLES
C     NOT ON THE REAL AXIS.  THE EFFICIENCY AND RELIABILITY
C     OF INTEGRATING SUCH FUNCTIONS MAY BE GREATLY IMPROVED BY MANUALLY
C     SUBDIVIDING THE INTERVAL OF INTEGRATION AT THE DISCONTINUITY,
C     SINGULARITY OR REAL PART OF THE POLE, AND SUMMING THE ANSWERS.
C     A CHANGE OF VARIABLE TO ELIMINATE OR REDUCE THE STRENGTH OF THE
C     SINGULARITY WILL SIGNIFICANTLY IMPROVE PERFORMANCE.
C
C     *****    FORMAL ARGUMENTS    *************************************
C
C ANSWER  IS THE COMPUTED INTEGRAL IF IFLAG .LT. 0.  IF IFLAG = 6,
C         ANSWER IS THE APPROXIMATE LOCATION OF A SINGULARITY WHICH
C         APPEARS NOT TO BE INTEGRABLE.  IF IFLAG = 0, ANSWER IS
C         USED TO PASS A FUNCTION VALUE FROM THE USER PROGRAM TO DINTA.
      DOUBLE PRECISION ANSWER
C WORK    CONTAINS THE ESTIMATED ABSOLUTE ERROR IF IFLAG .LT. 0.
C         IF IFLAG = 0, WORK IS USED TO PASS AN ABSCISSA TO THE
C         USER PROGRAM (FOR REVERSE COMMUNICATION).
      DOUBLE PRECISION WORK(*)
C IFLAG   IS USED TO CONTROL REVERSE COMMUNICATION AND TO RETURN
C         STATUS INFORMATION TO THE USER.  VALUES ARE
C          0 - INTEGRAND VALUE NEEDED.  COMPUTE ANSWER=INTEGRAND
C              EVALUATED AT WORK(1), ..., WORK(NDIM).
C         -1 - NORMAL TERMINATION, EITHER THE ABSOLUTE OR THE RELATIVE
C              ERROR CRITERIA ARE SATISFIED.
C         -2 - NORMAL TERMINATION, NEITHER THE ABSOLUTE NOR RELATIVE
C              ERROR CRITERIA ARE SATISFIED, BUT THE RELATIVE TO THE
C              OBTAINABLE ERROR CRITERION IS SATISFIED.
C         -3 - NORMAL TERMINATION, BUT NONE OF THE ERROR CRITERIA ARE
C              SATISFIED.
C          4 - BAD IOPT VALUE (SEE DINT1 OR DINTM).
C          5 - TOO MANY FUNCTION VALUES NEEDED.
C          KDIM+5 - APPARENT NON-INTEGRABLE SINGULARITY NEAR ANSWER.
C                   KDIM IS CURRENT DIMENSION - 1 IN 1 DIMENSIONAL CASE.
      INTEGER IFLAG(*)
C
C     *****    EXTERNAL REFERENCES   ***********************************
C
C DINTDL  TO FORM AND ANALYZE DIFFERENCE LINES.
C
C DINTDU  TO UPDATE DIFFERENCE LINES.
C
C DINTF   PROVIDE THE INTEGRAND.  DINTF IS EXPLAINED IN
C         EXPLAINED IN DINT1 AND DINTM.
C
C DINTNS  TO INCREASE AND REDUCE NSUB, THE DEGREE OF THE INDEPENDENT
C         VARIABLE TRANSFORMATION.
C
C DINTO   USED FOR PRINTING OUTPUT.
C
C DINTSM  USED TO CALCULATE THE MINIMUM STEPSIZE.
      EXTERNAL DINTSM
      DOUBLE PRECISION DINTSM
      INTEGER NSRB,NSRA,NSIB,NSIA,MAXK,NSUBMX,KONVRG,IC1,IC2,JUMPS
C
C     *****     PARAMETERS     *****************************************
C
C NSIA    AS FOR NSRA, BUT NSUB IS TO BE INCREASED.
C NSIB    AS FOR NSRB, BUT NSUB IS TO BE INCREASED.
C NSRA    IS AN ARGUMENT OF DINTNS THAT MEANS REDUCE NSUB, BUT ONLY
C         ALOCAL NEEDS TO BE RECALCULATED.
C NSRB    AS FOR NSRA, BUT BOTH ALOCAL AND BLOCAL MUST BE RECALCULATED.
C
      PARAMETER (NSRB=1, NSRA=NSRB+1, NSIB=NSRA+1, NSIA=NSIB+1)
C
C     *****     INTERNAL AND COMMON VARIABLES     **********************
C
C AACUM   IS THE INTEGRAL OF THE ABSOLUTE VALUE OF THE INTEGRAND OVER
C         CURRENT PANEL.
C ABSCIS  IS THE ABSCISSA TO BE TRANSFORMED TO USER COORDINATES.
C         SEE NSUB.
C ABSDIF  IS THE ABSOLUTE VALUE OF ONE HALF OF THE LENGTH OF THE PANEL
C         (SMALLEST SUBDIVISION) CURRENTLY BEING INTEGRATED.
C ACUM    IS THE ESTIMATE OF THE INTEGRAL OVER THE CURRENT PANEL.
C AIMFOR  IS AN INDEX OF THE FORMULA TO AIM FOR AFTER A PANEL HAS BEEN
C         SUCCESSFULLY INTEGRATED.  THE CURRENT VALUE OF AIMFOR
C         SPECIFIES A 63-POINT FORMULA.
      INTEGER AIMFOR
C AINIT   IS THE INITIAL LOWER REGION BOUNDARY.
C ALOCAL  IS THE 'LEFT' BOUNDARY OF THE PANEL BEING INTEGRATED.
      DOUBLE PRECISION ALOCAL
C ALPHA   IS USED DURING DETERMINATION OF THE NEXT STEP SIZE AS AN
C         INDICATOR OF THE PERFORMANCE OF THE INTEGRATION OF THE
C         PREVIOUS STEP, WITH RESPECT TO THE COMMITTED AND PERMITTED
C         ERROR.
      DOUBLE PRECISION ALPHA
C BETA    IS USED DURING COMPUTATION OF ABSCISSAE USED BY THE
C         QUADRATURE FORMULA AND DURING COMPUTATION OF THE NEXT STEP.
      DOUBLE PRECISION BETA
C BINIT   IS THE INITIAL UPPER REGION BOUNDARY.
C BLOCAL  IS THE 'RIGHT' BOUNDARY OF THE PANEL CURRENTLY BEING
C         INTEGRATED.
      DOUBLE PRECISION BLOCAL
C COUNT   IS USED DURING A SEARCH TO DISTINGUISH INTEGRABLE FROM
C         NON-INTEGRABLE SINGULARITIES.
C DELMIN  IS THE MINIMUM STEPSIZE ALLOWED, LEST THE ABSCISSAE COALESCE.
C DELTA   IS THE ABSOLUTE VALUE OF THE CURRENT STEP SIZE.
C DID1    IS INITIALLY FALSE, AND SET TRUE WHEN ONE PANEL HAS
C         BEEN SUCCESSFULLY INTEGRATED.
C DIFF    IS ONE HALF OF THE SIGNED CURRENT STEP SIZE.
C DISCF   IS NON-ZERO WHEN A DISCONTINUITY IS DETECTED.  IF DISCF IS
C         POSITIVE AND IPRINT IS NOT ZERO AFTER SUCCESSFULLY
C         INTEGRATING A PANEL IN PART = 1, THE BOUNDARIES FOR THE
C         PANEL ARE PRINTED.  IN ADDITION, IF ABS(DISCF) IS 2 AFTER A
C         PANEL IS INTEGRATED, THE REST OF THE PART IS USED FOR THE
C         NEXT PANEL.
C DISCHK  IS SET TO MAX(DISCHK,1) WHEN A JUMP IS SEEN IN THE DIFFERENCE
C         LINE, BUT THE NEXT FORMULA IS TRIED, OR WHEN A DISCONTINUITY
C         IS DETECTED.  IF THERE ARE JUMPS IN THE DIFFERENCE LINE WHICH
C         SUBSEQUENTLY VANISH, DISCHK IS SET TO 2.  WHEN A PANEL IS
C         ACCEPTED, DISCHK IS DECREMENTED AS FAR AS ZERO.  DISCHK IS
C         SET TO -1 WHEN NSUB IS CHANGED.  DISCHK IS SET TO -2 WHEN A
C         CAUTION POINT IS ENCOUNTERED.  IF DISCHK IS POSITIVE, THE
C         DIFFERENCE LINES MUST BE FORMED EVERY QUADRATURE STEP, AND
C         THE END POINTS MUST BE ADDED BEFORE THE ANSWER CAN BE
C         ACCEPTED.  IF DISCHK IS -1, THE DIFFERENCE LINE MUST BE
C         FORMED BEFORE THE ANSWER CAN BE ACCEPTED.  IF DISCHK IS -2
C         THE DIFFERENCE LINES MUST BE FORMED AND THE END POINTS ADDED
C         BEFORE THE ANSWER CAN BE ACCEPTED.
C DISCX   IS THE APPROXIMATE LOCATION OF A DISCONTINUITY.
C EDUE2A, EDUE2B  ARE THE ERRORS IN THE INTEGRAL DUE TO IMPRECISE LIMITS
C         A AND B RESPECTIVELY.
C END     STORES THE 'RIGHT' BOUNDARIES OF MAJOR SUBDIVISIONS.
C ENDPHI  CONTAINS THE END OF THE DIFFERENCE LINE ON WHICH THERE CAN BE
C         A JUMP.  ENDPHI(1) = 1 BECAUSE THERE CAN ONLY BE A JUMP IN PHI
C         AT 1.  ENDPHI(2) = LENDT BECAUSE THERE CAN ONLY BE A JUMP IN
C         PHIT AT LENDT.  ENDPHI IS USED ONLY WHEN CHECKING WHETHER THE
C         ERROR IN THE FUNCTION AT THE END OF A PANEL IS SO LARGE THAT
C         A JUMP IS NOT NECESSARILY BELIEVABLE.
      INTEGER ENDPHI(2)
C ENDPTS  IS USED TO CONTROL ADDITION OF THE END POINTS TO THE
C         DIFFERENCE LINES:
C           1.  FUNCTION NEED NOT BE EVALUATED AT A.
C           2.  FUNCTION NEED NOT BE EVALUATED AT B.
C           3.  FUNCTION NEED NOT BE EVALUATED AT EITHER A AND B.
C EP      IS THE PREVIOUS VALUE OF ERRC.
C EPNOIZ  IS THE MAXIMUM ERROR THRESHOLD FOR NOISE DETECTION.
C EPS     IS THE ERROR ALLOWED ON THE CURRENT PANEL.
C EPSMAX  IS THE MOST STRINGENT TOLERANCE ALLOWED ON THE CURRENT PANEL.
C         SEE XEPS.
C EPSMIN  IS AN ESTIMATE OF THE SMALLEST ERROR POSSIBLE ON THE CURRENT
C         PANEL.
C EPSO    IS THE ORIGINAL ERROR TOLERANCE.
C EPSR    IS THE REMAINING ERROR TOLERANCE.
C EPSS    IS A STORAGE CELL FOR THE MAXIMUM OF EPS AND EPSMAX.  EPS IS
C         CHANGED AFTER EPSS IS COMPUTED.
C ERR     IS THE EXTRAPOLATED ERROR COMMITTED DURING THE CURRENT
C         INTEGRATION STEP.
C ERRAT   ERRAT(I) CONTAINS THE ERROR IN FAT(I) IF FATAS(I), I=1 OR 2.
C ERRC    IS THE ERROR COMMITTED DURING THE CURRENT INTEGRATION STEP,
C         FOR PURPOSES OF DETERMINING THE NEXT STEPSIZE.
C ERRCF   IS USED TO COMPUTE A MINIMUM VALUE FOR THE ERROR
C         WHEN K.LE.4, USING THE DIFFERENCE LINES.
      DOUBLE PRECISION ERRCF(3)
C ERRF    IS THE ESTIMATED ERROR IN THE INTEGRAND.  WHEN KDIM=1, ERRF
C         IS OBTAINED FROM WORK(FEA) (IF FEA .NE. 0).  WHEN KDIM.GT.1
C         ERRF IS THE ESTIMATED ERROR IN THE INNER INTEGRAL.
C ERRI    IS THE ERROR ESTIMATED BY THE DIFFERENCE BETWEEN TWO
C         SUCCESSIVE INTEGRATION FORMULAE APPLIED TO THE CURRENT PANEL.
C ERRINA, ERRINB  ARE THE ERRORS THE USER DECLARES BY WAY OF AN OPTION
C         TO BE PRESENT IN THE LIMITS A AND B RESPECTIVELY.
C ERRT    IS THE TOTAL ERROR COMMITTED ON A MAJOR SUBDIVISION.
C ESOLD   IS THE PREVIOUS ERROR ESTIMATED ACROSS A JUMP BY THE SEARCH.
C EXTRA   IS ANY EXTRA AMOUNT ADDED TO THE ERROR BECAUSE OF OUR
C         SUSPICIOUS NATURE.  EXTRA IS BASED UPON HIGH-ORDER
C         DIFFERENCES IF THERE ARE NO JUMPS.
C FAIL    INDICATES TO THE STEPSIZE SELECTION ROUTINE THAT INTEGRATION
C         WAS NOT SUCCESSFUL.  FAIL IS ALSO SET IF A JUMP IN THE
C         DIFFERENCE LINE VANISHES WHEN ACCURATE LOCATION IS ATTEMPTED.
C FAT     IS A VECTOR EQUIVALENT TO FATA,FATB.
C FATA    SAVES THE FUNCTION VALUE AT ALOCAL.
      DOUBLE PRECISION FATA
C FATAS   IS TRUE WHEN FATA AND ERRAT(1) CONTAIN VALID DATA.
      LOGICAL FATAS
C FATB    SAVES THE FUNCTION VALUE AT BLOCAL.
      DOUBLE PRECISION FATB
C FATBS   IS TRUE WHEN FATB AND ERRAT(2) CONTAIN VALID DATA.
      LOGICAL FATBS
C FATS    IS A VECTOR EQUIVALENT TO FATAS,FATBS.
C FEA     IS THE INDEX IN WORK OF THE ABSOLUTE ERROR COMMITTED
C         COMPUTING F.  WORK(FEA) MAY BE CHANGED DURING THE INTEGRATION.
C         FEA APPLIES ONLY TO THE INNERMOST INTEGRATION OF A MULTIPLE
C         QUADRATURE.
C FER     IS THE RELATIVE ERROR COMMITTED COMPUTING F.
C FNCVAL  IS THE COMPUTED FUNCTION VALUE.
C FSAVE   IS USED TO STORE THE FUNCTION VALUE AT THE ABSCISSA WHERE
C         THE INTERVAL IS SUBDIVIDED.
C FSAVED  IS SET IF FSAVE CONTAINS USEFUL DATA.
C FSTORE  IS USED TO DETERMINE WHICH FUNCTION VALUES ARE TO BE STORED,
C         AND WHERE THEY ARE TO BE STORED IN THE FUNCTION TABLE.
C         FSTORE IS READ-ONLY.  THE CONTENT OF EACH ELEMENT OF FSTORE
C         IS A TWO DIGIT BASE-100 NUMBER, WHERE THE FIRST DIGIT
C         INDICATES WHERE THE FIRST FUNCTION VALUE COMPUTED DURING THE
C         CURRENT INTEGRATION STEP IS TO BE STORED, AND THE SECOND DIGIT
C         INDICATES THE LAST CELL OF THE FUNCTION TABLE TO BE USED FOR
C         THE CURRENT INTEGRATION STEP.
      INTEGER FSTORE(7)
C FT      USED TO HOLD, TEMPORARILY, THE FUNCTION VALUES USED.
C         THESE FUNCTION VALUES ARE USED TO COMPUTE THE DIFFERENCE
C         LINE.  SEE XT.
C FUDGE   IS USED DURING THE CONVERGENCE TEST TO RELATE THE ESTIMATED
C         ERROR TO THE ROUND-OFF LEVEL.  FUDGE(K-KAIMT) GIVES THE
C         RELATIVE PRECISION THAT MUST BE OBTAINED IF THE ERROR
C         ESTIMATE IS NOT TO BE INCREASED.
      DOUBLE PRECISION FUDGE(5)
C FUNCT   IS USED TO STORE FUNCTION VALUES WHICH MAY BE RE-USED BY
C         THE PATTERSON FORMULAE, AND ADDITIONAL FUNCTION VALUES TO
C         COMPUTE THE DIFFERENCE LINE.
C         FUNCT(1-17)  IS STORAGE NEEDED FOR THE PATTERSON METHOD
C         FUNCT(18-24) IS STORAGE TO RECOVER F(X1) THROUGH F(X15).
C F1, F2  CONTAIN NEWLY COMPUTED FUNCTION VALUES DURING THE QUADRATURE
C         STEP.  THEY ARE ALSO USED ELSEWHERE FOR TEMPORARY STORAGE.
C GAMMA   IS USED DURING COMPUTATION OF THE NEXT STEPSIZE.
      DOUBLE PRECISION GAMMA(7)
C HAVDIF  IS SET TO INDICATE THAT THE DIFFERENCE LINE HAS BEEN COMPUTED
C         USING 15 FUNCTION VALUES.
C I       IS USED FREELY AS AN INDEX.
C IC1 IC2 ARE THE INDICES OF CONVERGENCE IN THE FORWARD AND BACKWARD
C         DIFFERENCE LINES.  IC1 AND IC2 ARE EQUIVALENT TO ISTOP(1,1)
C         AND ISTOP(1,2). (SEE BELOW).
C IEND    IS SET WHEN THE END OF A MAJOR SUBDIVISION IS IN THE CURRENT
C         PANEL.
C IH      ARE INDICES USED TO LOCATE ABSCISSAE IN THE TABLE OF
C         QUADRATURE RULES.  IH IS USED DURING COMPUTATION OF THE
C         DIFFERENCE LINE.
      INTEGER IH(7)
C INC INC2 ARE USED DURING ANALYSIS OF THE DIFFERENCE LINE AS INDEX
C         INCREMENTORS.  THEY ARE ALSO USED FREELY AS INDICES.
C INEW    IS USED DURING THE QUADRATURE STEP TO CONTROL A LOOP.  INEW
C         MAY NOT BE CHANGED BETWEEN QUADRATURE STEPS.
C INIT    IS SET TO INDICATE THAT THE ORIGINAL INTERVAL IS BEING
C         INTEGRATED.
C INSTOP  IS USED TO CONTROL A LOOP DURING ANALYSIS OF THE DIFFERENCE
C         LINE.
C IOLD    IS USED DURING THE QUADRATURE STEP TO CONTROL A LOOP.  IOLD
C         MAY NOT BE CHANGED BETWEEN QUADRATURE STEPS.
C IP      IS USED DURING THE QUADRATURE STEP TO INDEX THE TABLE OF
C         ABSCISSAE AND WEIGHTS.  IP MAY NOT BE CHANGED BETWEEN
C         QUADRATURE STEPS.
C IPRINT  CONTROLS THE LEVEL OF DIAGNOSTIC PRINTING.
C         0 - NO PRINTING
C         1 - MINIMUM PRINTING - ERROR MESSAGES (DEFAULT)
C         2 - PANEL BOUNDARIES AND ANSWERS
C         3 - ERROR ESTIMATES FOR EACH QUADRATURE FORMULA
C         4 - DETAILED OUTPUT (DIFFERENCE LINES, ETC).
C ISTOP   HOLDS RESULTS OF ANALYSIS OF DIFFERENCE LINE.
C         ENTRIES IN ISTOP MEAN
C         ISTOP(1,1) - INDEX OF XT TO WHICH CONVERGENCE EXTENDS FROM
C              THE ALOCAL END OF THE INTERVAL,
C         ISTOP(1,2) - INDEX OF XT TO WHICH CONVERGENCE EXTENDS FROM
C              THE BLOCAL END OF THE INTERVAL
C         ISTOP(2,1) AND ISTOP(2,2) MAY INDICATE INDICES OF XT WHERE
C              THE INTERVAL SHOULD BE SUBDIVIDED.
C IXKDIM  INDEXES THE CELL IN IFLAG IN WHICH THE CURRENT DIMENSION OF A
C         MULTIDIMENSIONAL INTEGRATION IS STORED.  THE DEFAULT IS 1,
C         CORRESPONDING WITH THE USUAL CELL USED FOR COMMUNICATION WITH
C         DINTF.  THE USER MAY NEED TO CHANGE THIS BY USING OPTION 12
C         IF NON-STANDARD DIMENSION CHANGES ARE NECESSARY, AND A CHANGE
C         OF VARIABLE IS SIMULTANEOUSLY NECESSARY.  SEE DINTM AND
C         DINTMA.
C J       IS FREELY USED AS AN INDEX.
C JUMPS   STORES THE QUALITATIVE RESULT OF THE DIFFERENCE LINE ANALYSIS.
C         1 - NO JUMPS
C         2 - ONE JUMP, ON THE END OF THE LINE
C         3 - ONE JUMP, NOT ON THE END
C         4 - TWO JUMPS.
C J1 J2   ARE USED DURING THE QUADRATURE STEP TO CONTROL STORING INTO
C         THE FUNCT ARRAY.  THEY ARE RECOMPUTED AT EACH STEP.  THEY
C         ARE ALSO USED DURING THE ANALYSIS OF THE DIFFERENCE LINE TO
C         STORE 'JUMP' INDICES.
C J1OLD   ARE USED TO REMEMBER J1 AND J2 DURING A SEARCH
C J2OLD   FOR A JUMP.
C K       IS THE INDEX FOR THE QUADRATURE FORMULA LAST USED.  THE
C         NUMBER OF POINTS USED IN THE KTH FORMULA IS 2**K - 1.
C KAIMT   IS THE INDEX OF THE QUADRATURE FORMULA EXPECTED TO CONVERGE
C         ON THE NEXT STEP.
C KDIM    IS THE CURRENT DIMENSION - ALWAYS 1 IN 1 DIMENSIONAL CASE.
C KK      IS A TEMPORARY INTEGER STORAGE.
C KMAX    IS THE MAXIMUM VALUE ALLOWED FOR K ON THE CURRENT INTEGRATION
C         STEP.
C KMAXF   IS THE VALUE TO USE FOR KMAX INITIALLY AND AFTER A FAILURE.
C KMIN    IS THE MINIMUM QUADRATURE FORMULA INDEX TO ACCEPT IF THERE
C         ARE NO JUMPS.  KMIN IS USUALLY ZERO, BUT IS SET TO K IF THE
C         STEP SIZE INCREASES.
C KONVRG  KONVRG IS THE TOTAL NUMBER OF NODES OF CONVERGENCE OBSERVED
C         IN THE DIFFERENCE LINES.  KONVRG=IC1+LENDT-IC2.
C KORECT  IS USED TO DETERMINE WHICH OF THE STORED FUNCTION VALUES ARE
C         TO BE USED TO CORRECT ONE-HALF OF THE PREVIOUSLY COMPUTED
C         INTEGRAL ON THE CURRENT INTERVAL.  THE INDICES ARE STORED IN
C         KORECT AS BASE-20 DIGITS.  KORECT IS READ-ONLY, INDEXED BY K.
      INTEGER KORECT(6)
C L       IS USED DURING ANALYSIS OF THE DIFFERENCE LINE.
C LENDT   IS THE CURRENT LENGTH OF THE DIFFERENCE LINE.
C LOCAL   IS A VECTOR EQUIVALENT TO ALOCAL,BLOCAL.
C         LOCAL(3),LOCAL(4) ARE ALOCAL,BLOCAL IN USER COORDINATES.
C MAXK    IS THE MAXIMUM VALUE ALLOWED FOR K.
C NEEDH   IS SET IF THE HEADER FOR DIAGNOSTIC INFORMATION IS NEEDED.
C NEWEPS  MEANS EPS MUST BE RECOMPUTED TO ACCOUNT FOR ERROR IN LIMITS
      LOGICAL NEWEPS
C NFEVAL  IS THE CURRENT NUMBER OF FUNCTION EVALUATIONS.
C NFJUMP  IS THE VALUE OF NFEVAL WHEN A JUMP IS FIRST NOTICED.
C NFMAX   IS THE INDEX OF THE ENTRY IN THE OPTION VECTOR (IOPT) TO
C         USE TO CONTROL THE MAXIMUM NUMBER OF FUNCTION EVALUATIONS.
C NOMOUT  IS THE NOMINAL OUTPUT UNIT.
C NSUB    IS THE POWER OF THE TRANSFORMATION OF THE FORM
C         X=TA+(T-TA)*((T-TA)/(TB-TA))**(NSUB-1)
C         WHERE X IS THE USER'S INDEPENDENT VARIABLE AND T IS THE
C         INDEPENDENT VARIABLE TO USE FOR INTEGRATION.  NSUB MAY
C         ONLY BE A POWER OF 2.  NSUB=0 REALLY MEANS NSUB=1, BUT
C         IS SIMPLER TO USE.
C NSUBMX  IS THE MAXIMUM VALUE FOR NSUB.
C NSUBSV  IS THE VALUE TO USE FOR NSUB WHEN RESTARTING
C         AFTER A SUBDIVISION.
C NXKDIM  IS THE VALUE OF KDIM BEFORE AN INNER INTEGRAL WAS REQUESTED.
C         THE INNER INTEGRAL IS NOT NECESSARILY OF DIMENSION KDIM-1.
C         SEE IXKDIM ABOVE, DINTM AND DINTMA.
C P       IS THE ARRAY OF ABSCISSAE AND WEIGHTS.  DATA FOR P IS INCLUDED
C         AFTER THE DESCRIPTION OF INTERNAL VARIABLES.
      DOUBLE PRECISION P(305)
C PAACUM  IS THE PREVIOUS VALUE OF AACUM.
C PACUM   IS THE PREVIOUS VALUE OF ACUM.  ABS(ACUM-PACUM) IS USED AS AN
C         INITIAL ERROR ESTIMATE.
C PART    IS THE INDEX OF THE CURRENT MAJOR SUBDIVISION.  START, END,
C         STEP, ERRT AND RESULT ARE SUBSCRIPTED BY PART.  PART BEGINS
C         SET TO 1, IS CHANGED TO 2 WHEN THE INTERVAL IS SUBDIVIDED,
C         AND SET BACK TO 1 WHEN PART 2 IS INTEGRATED.  WHEN PART 1
C         IS FINISHED, THE ENTIRE PROBLEM IS FINISHED.
C PEPSMN  IS THE PREVIOUS VALUE OF EPSMIN.  PEPSMN MAY NOT BE CHANGED
C         BETWEEN QUADRATURE STEPS.
C PERR    IS THE ERROR COMMITTED ON THE LAST SUCCESSFULLY INTEGRATED
C         PANEL.  PERR IS USED DURING THE NOISE TEST.
C PF1 PF2 ARE THE PREVIOUS VALUES OF F1 AND F2.  PF1 AND PF2 ARE USED TO
C         COMPUTE EPSMIN.
C PHI     CONTAINS THE DIFFERENCE LINES.  PHI IS DIMENSIONED AND
C         EQUIVALENCED WITH PHIT AFTER THE DESCRIPTION OF INTERNAL
C         VARIABLES.
C PHIS    IS A TWO-DIMENSIONAL ARRAY EQUIVALENCED TO PHI AND PHIT.  IT
C         IS USED TO ACCESS THE DIFFERENCE LINES USING AN INDEX.
      DOUBLE PRECISION PHIS(17,2)
C PHISUM  IS THE SUM OF 3 MIDDLE ENTRIES OF PHI.
C PHIT    IS THE BACKWARD DIFFERENCE LINE.  SEE PHI.
      DOUBLE PRECISION PHIT(17)
C PHTSUM  IS THE SUM OF 3 MIDDLE ENTRIES OF PHIT.
C PX      IS THE PREVIOUS VALUE OF THE OFFSET FROM THE ENDS OF THE
C         CURRENT INTERVAL.  PX IS USED TO COMPUTE EPSMIN.
C RE      IS THE RATIO OF THE ERROR COMMITTED WITH THE CURRENT
C         QUADRATURE FORMULA TO THAT COMMITTED WITH THE PREVIOUS
C         QUADRATURE FORMULA.  RE IS USED TO DETECT CONVERGENCE AND
C         TO AVOID FALSE CONVERGENCE, AND TO DETECT NOISE.
C RECORR  IS USED TO CORRECT RE WHEN NSUB IS NOT ZERO.  ONLY THE FIRST
C         AND THIRD ELEMENTS OF RECORR ARE USED EXPLICITLY.
      DOUBLE PRECISION RECORR(3)
C RELEPS  IS THE ABSOLUTE ERROR TOLERANCE, BASED UPON THE RELATIVE
C         RELATIVE TOLERANCE REQUESTED AND THE BEST AVAILABLE ESTIMATE
C         OF THE ANSWER.
C RELOBT  IS THE PRECISION RELATIVE TO THE OBTAINABLE PRECISION.
C         THE ESTIMATED ERROR IS COMPARED TO
C         (EPSMIN/AACUM)**RELOBT*AACUM.
C RELOOK  IS USED TO FORCE EXAMINATION OF THE DIFFERENCE LINES, IF
C         THE METHOD IS TO ACCEPT AN ANSWER BUT RE*REP .GT. RELOOK.
      DOUBLE PRECISION RELOOK
C RELTOL  IS THE INDEX INTO WORK OF THE PARAMETERS FOR ABSOLUTE, LOCAL
C         RELATIVE, AND GLOBAL RELATIVE TOLERANCE PARAMETERS.
C REPROD  IS THE PRODUCT OF ALL VALUES OF RE FOR THE CURRENT PANEL.
C         REPROD IS NOT ALLOWED TO BE GREATER THAN 1.0.
C REP     IS THE PREVIOUS VALUE OF RE.  IT IS USED SIMILARLY TO RE.
C RESULT  IS THE CUMULATIVE INTEGRAL FOR THE CURRENT MAJOR SUBDIVISION.
C REVERS  IS THE INDEX OF THE ENTRY IN THE OPTION VECTOR (IOPT) TO
C         USE TO CONTROL REVERSE COMMUNICATION.
C RNDC    IS THE NOISE OR ROUND-OFF LEVEL, INITIALLY EQUAL TO FER.
C         RNDC IS CHANGED IF NOISE IS DETECTED.
C ROUNDF  IS SET TO INDICATE THAT ROUND-OFF IS LIMITING THE ERROR, OR
C         THAT THE INTEGRAND IS BEING TREATED AS PURE NOISE WITH RESPECT
C         TO THE REQUESTED ERROR TOLERANCE.
C S       IS A  SCRATCH VARIABLE.
C SEARCH  IS USED TO INDICATE WHETHER POINTS NOT USED BY THE INTEGRATION
C         FORMULAE ARE BEING ADDED TO THE DIFFERENCE LINES.
C         1 - POINTS NOT BEING ADDED
C         2 - SEARCHING FOR INITIAL INTERVAL SIZE
C         3 - MAKING SURE AN END-POINT JUMP IS NOT INSIDE THE INTERVAL
C         4 - SEARCHING FOR A DISCONTINUITY
C         5 - CHANGING FROM SEARCH=2 OR 3 TO SEARCH=4
C         6 - NSUB WAS JUST REDUCED BECAUSE WE SAW JUMPS.
C             NSUB CANNOT BE INCREASED WITHOUT DOING A TYPE-2 SEARCH.
C         NEGATIVE - INITIALIZE THE SEARCH (TRANSIENT STATE).
C SMIN    IS THE MINIMUM STEPSIZE TO USE FOR A PANEL BEGINNING AT SX.
      DOUBLE PRECISION SMIN
C SPACE   RESERVES SOME SPACE IN /DINTC/ IN CASE IT IS NEEDED.  THIS
C         IS DONE SO USERS' CODE WON'T NEED TO CHANGE.
C START   STORES THE 'LEFT' BOUNDARY OF A PART.
C STEP    STORES THE SIZE OF THE FIRST STEP TO USE ON A MAJOR
C         SUBDIVISION.  THE SIGN OF STEP IS THE SIGN TO USE FOR
C         ACCUMULATION OF THE ANSWER INTO RESULT(PART).
C SUM     IS THE CENTER OF THE CURRENT INTERVAL.
C T       IS THE TRANSFORMED INDEPENDENT VARIABLE.  SEE NSUB.
C TA TB   ARE THE INTERVAL ON WHICH A VARIABLE TRANSFORMATION
C         IS DEFINED.  SEE NSUB.  THE SIGN OF TB IS ALSO USED TO
C         INDICATE THE DIRECTION OF STEP ACCUMULATION.
C TALOC   IF NON-ZERO IS THE INDEX IN WORK OF THE LOCATION OF A
C         SINGULARITY.
C TASAVE  IS THE VALUE TO USE FOR TA WHEN RESTARTING PART 1
C         AFTER SUBDIVIDING.
C TDECR   IS AN ARITHMETIC STATEMENT FUNCTION DEFINED BELOW.
      DOUBLE PRECISION TDECR
C TEND    IS THE POINT AT WHICH USAGE OF A T**N SUBSTITUTION SHOULD BE
C         DISCONTINUED.
C THROW   IS USED DURING THE SEARCH FOR A DISCONTINUITY TO INDICATE
C         WHICH END OF THE DIFFERENCE LINES SHOULD BE DISCARDED.
      INTEGER THROW
C TINCR   IS AN ARITHMETIC STATEMENT FUNCTION DEFINED BELOW.
      DOUBLE PRECISION TINCR
C TLEN    IS THE TOTAL LENGTH OF THE UN-INTEGRATED PART.  TLEN IS USED
C         TO PARTION ERROR ALLOTTMENTS.
C TP TP1  ARE TEMPORARY VARIABLES.
      DOUBLE PRECISION TP1
C TPS     IS A TEMPORARY VARIABLE.
C WHERE   IS USED AS A COMPUTED GO TO INDEX.
C WHERE2  IS USED AS A COMPUTED GO TO INDEX.
C WORRY   IS USED TO STORE THE CAUTION POINTS FOR EACH PART.
C X       DURING THE QUADRATURE STEP, X IS THE DISTANCE FROM THE ENDS OF
C         THE INTERVAL, AT WHICH ABSCISSAE THE INTEGRAND IS EVALUATED.
C         X IS ALSO USED ELSEWHERE FOR TEMPORARY STORAGE.
C XEPS    IS THE MINIMUM PART OF THE REMAINING ERROR TOLERANCE TO
C         REQUIRE ON THE CURRENT INTERVAL.
      DOUBLE PRECISION XEPS
C XJ      IS THE MAGNITUDE OF THE LARGEST JUMP IN THE BACKWARD
C         DIFFERENCE LINE.
C XJP     IS THE MAGNITUDE OF THE LARGEST JUMP IN THE FORWARD
C         DIFFERENCE LINE.
C XJUMP   IS A THRESHOLD FOR JUMPS IN THE DIFFERENCE LINE.
C XJUMPS  IS THE VALUE OF XJUMP TO USE DURING A SEARCH.
      DOUBLE PRECISION XJUMPS
C XSTEP   IS USED TO INCREASE THE CURRENT EFFECTIVE STEPSIZE AFTER
C         SUCCESSFULLY INTEGRATING AN INTERVAL.  THE CURRENT EFFECTIVE
C         STEPSIZE IS USED TO COMPUTE THE NEXT STEPSIZE.  XSTEP IS
C         READ-ONLY, INDEXED BY MAX(AIMFOR-K+3,1).
      DOUBLE PRECISION XSTEP(7)
C XT      THE ABSCISSAE OF THE QUADRATURE FORMULA.  XT ARE USED TO
C         COMPUTE THE DIFFERENCE LINE.
C X1 X2   ARE TEMPORARY STORAGE.
C Z       ARE ADJUSTABLE PARAMETERS USED FOR DEVELOPMENT.
      DOUBLE PRECISION Z(4)
C ZL1     IS THE ARGUMENT OF ARITHMETIC STATEMENT FUNCTIONS.
      DOUBLE PRECISION ZL1
C
C     *****    COMMON STORAGE ******************************************
C
C     COMMON /DINTNC/ CONTAINS VARIABLES NOT SEPARATELY SAVED FOR
C     EACH DIMENSION OF A MULTIPLE QUADRATURE.  COMMON /DINTC/
C     CONTAINS VARIABLES THAT MUST BE SAVED FOR EACH DIMENSION OF THE
C     QUADRATURE.  THE VARIABLES IN EACH COMMON BLOCK ARE STORED IN THE
C     ORDER - ALWAYS DOUBLE, DOUBLE IF DOUBLE PRECISION PROGRAM, DOUBLE
C     IF DOUBLE PRECISION PROGRAM AND EXPONENT RANGE OF DOUBLE AND
C     SINGLE VERY DIFFERENT, SINGLE, INTEGER, LOGICAL.  A PAD OF LOGICAL
C     VARIABLES IS INCLUDED AT THE END OF /DINTC/.  THE DIMENSION OF
C     THE PAD MAY NEED TO BE VARIED SO THAT NO VARIABLES BEYOND THE END
C     OF THE COMMON BLOCK ARE ALTERED.
C
C     DECLARATIONS OF COMMON /DINTNC/ VARIABLES.
C
      DOUBLE PRECISION AINIT, BINIT, FNCVAL, S, TP
      DOUBLE PRECISION FER, FER1, RELOBT, TPS, XJ, XJP
      INTEGER     FEA,       FEA1,      INC,       INC2,      IPRINT,
     1 ISTOP(2,2),JPRINT,    KDIM,      KK,        KMAXF,     NDIM,
     2 NFINDX,    NFMAX,     NFMAXM,    RELTOL,    REVERM,    REVERS,
     3 WHEREM
      LOGICAL NEEDH
C
C     DECLARATIONS OF COMMON /DINTC/ VARIABLES.
C
c--D Next line special: S => D, X => Q, D => D, P => D
      DOUBLE PRECISION ACUM, PACUM, RESULT(2)
C     139 $.TYPE.$ VARIABLES
      DOUBLE PRECISION
     1 AACUM,     ABSCIS,    DELMIN,    DELTA,     DIFF,      DISCX(2),
     2 END(2),    ERRINA,    ERRINB,    FAT(2),    FSAVE,
     3 FUNCT(24), F1,        F2,        LOCAL(4),  PAACUM,    PF1,
     4 PF2,       PHISUM,    PHTSUM,    PX,        SPACE(6),
     5 STEP(2),   START(2),  SUM,       T,         TA,        TASAVE,
     6 TB,        TEND,      WORRY(2),  X,         X1,
     7 X2,        XT(17),    FT(17),    PHI(34)
c Note XT, FT, and PHI above are last, because they must be in adjacent
c locations in DINTC.
C     30 $DSTYP$ VARIABLES
      DOUBLE PRECISION
     1 ABSDIF,    COUNT,     EDUE2A,    EDUE2B,    EP,        EPNOIZ,
     2 EPS,       EPSMAX,    EPSMIN,    EPSO,      EPSR,      EPSS,
     3 ERR,       ERRAT(2),  ERRC,      ERRF,      ERRI,      ERRT(2),
     4 ESOLD,     EXTRA,     PEPSMN,    RE,        RELEPS,    REP,
     5 REPROD,    RNDC,      TLEN,      XJUMP
C     29 INTEGER VARIABLES
      INTEGER     DISCF,     DISCHK,    ENDPTS,    I,         INEW,
     1 IOLD,      IP,        IXKDIM,    J,         J1,        J1OLD,
     2 J2,        J2OLD,     K,         KAIMT,     KMAX,      KMIN,
     3 L,         LENDT,     NFEVAL,    NFJUMP,    NSUB,      NSUBSV,
     4 NXKDIM,    PART,      SEARCH,    TALOC,     WHERE,     WHERE2
C     11 TO 18 LOGICALS (7 ARE PADDING).
      LOGICAL     DID1,      FAIL,      FATS(2),   FSAVED,    HAVDIF,
     1 IEND,      INIT,      ROUNDF,    XCDOBT(2), PAD(7)
C
C     THE COMMON BLOCKS.
C
      COMMON /DINTNC/
c        1       2       3     4        5       6       7        8
     W AINIT,  BINIT,  FNCVAL, S,      TP,     FER,    FER1,   RELOBT,
c       9      10       11      12      13       1       2        3
     X TPS,    XJ,     XJP,    FEA,    FEA1,   KDIM,    INC,    INC2,
c     4 (2,2)    8       9     10       11      12       13      14
     Y ISTOP,  JPRINT, IPRINT, KK,     KMAXF,  NDIM,   NFINDX, NFMAX,
c        15     16       17      18      19      20
     Z NFMAXM, RELTOL, REVERM, REVERS, WHEREM, NEEDH
      COMMON /DINTC/
     1 ACUM,   PACUM,  RESULT
      COMMON /DINTC/
c        1     2 (4)     6      7        8       9      10     11 (2)
     1 AACUM,  LOCAL,  ABSCIS, TA,     DELTA,  DELMIN, DIFF,   DISCX,
c     13 (2)     15      16    17 (2)   19     20 (24) 44
     2 END,    ERRINA, ERRINB, FAT,    FSAVE,  FUNCT,  F2,
c       45      46     47       48      49     50      51 (6)
     3 PAACUM, PF1,    PF2,    PHISUM, PHTSUM, PX,     SPACE,
c      57 (2)  59 (2)   61     62        63    64       65
     4 STEP,   START,  SUM,    T,      TASAVE, TB,     TEND,
c      66 (2)  68      69      70      71       72
     5 WORRY,  X1,     X2,     X,      F1,     COUNT,
c      73 (17) 90 (17) 107 (34)
     6 XT,     FT,     PHI
      COMMON /DINTC/
c       141     142    143     144      145     146
     1 ABSDIF, EDUE2A, EDUE2B, EP,     EPNOIZ, EPSMAX,
c       147     148     149    150 (2)  152     153
     2 EPSO,   EPSR,   EPSS,   ERRAT,  ERRC,   ERRF,
c     154 (2)   156     157     158     159    160
     3 ERRT,   ESOLD,  EXTRA,  PEPSMN, RELEPS, REP,
c       161     162     163
     4 RNDC,   TLEN,   XJUMP,
c       164    165      166    167    168       169
     5 ERRI,   ERR,    EPSMIN, EPS,    RE,     REPROD
      COMMON /DINTC/
c       170     171     172
     1 DISCF,  DISCHK, ENDPTS, INEW,   IOLD,   IP,     IXKDIM,
     2 J,      J1,     J1OLD,  J2,     J2OLD,  KMAX,   KMIN,
     3 L,      LENDT,  NFJUMP, NSUBSV, NXKDIM, TALOC,  WHERE2,
c      1       2          3      4       5         6      7       8
     4 I,      K,      KAIMT,  NSUB,   PART,   SEARCH, WHERE, NFEVAL
      COMMON /DINTC/
     1 DID1,   FAIL,   FATS,   FSAVED, HAVDIF, IEND,   INIT,   ROUNDF,
     2 XCDOBT, PAD
      SAVE /DINTNC/, /DINTC/
C
C     THE VARIABLES HERE DEFINE THE MACHINE ENVIRONMENT.  ALL ARE SET
C     IN DINTOP.  THE MEANING ATTACHED TO THESE VARIABLES CAN BE
C     FOUND BY LOOKING AT THE DEFINITIONS IN DINTOP.
      DOUBLE PRECISION
     1  EMEPS,  EEPSM8, EDELM2, EDELM3, ESQEPS, ERSQEP, ERSQE6, EMINF,
     2  ESMALL, ENZER,  EDELM1, ENINF
      COMMON /DINTEC/
     1  EMEPS,  EEPSM8, EDELM2, EDELM3, ESQEPS, ERSQEP, ERSQE6, EMINF,
     2  ESMALL, ENZER,  EDELM1, ENINF
      SAVE /DINTEC/
C
C     *****    EQUIVALENCE STATEMENTS    *******************************
C
      EQUIVALENCE (PHI(18),PHIT),(PHI(1),PHIS(1,1))
      EQUIVALENCE (FAT(1),FATA), (FAT(2),FATB)
      EQUIVALENCE (FATS(1),FATAS), (FATS(2),FATBS)
      EQUIVALENCE (ISTOP(1,1),IC1), (ISTOP(1,2),IC2)
      EQUIVALENCE (LOCAL(1),ALOCAL), (LOCAL(2),BLOCAL)
C
C     *****    DATA STATEMENTS   ***************************************
C
      DATA AIMFOR /6/
      DATA ENDPHI(1) /1/
      DATA ERRCF /1.579D0,5.233D-2,8.501D-3/
      DATA FSTORE /202,304,508,916,1217,1417,100/
      DATA FUDGE /1.0D-8,1.0D-6,1.0D-5,1.0D-4,3.0D-4/
C     K-KAIMT =     -4      -3     -2     -1      0
      DATA GAMMA /1.42d0,2.28d0,4.62d0,19.00d0,75.00d0,75.00d0,75.00d0/
      DATA IH /14,6,16,2,18,8,20/
      DATA KORECT /41,81,161,321,55890063,56395705/
      DATA MAXK /8/
      DATA NSUBMX /4/
      DATA RECORR /1.0D-3, 0.0D0, 1.0D-1/
      DATA RELOOK /5.0D-3/
      DATA XEPS / 0.02D0 /
      DATA XJUMPS /0.5D0/
      DATA XSTEP /0.5d0,0.875d0,1.0d0,2.25d0,5.0d0,15.0d0,50.0d0/
      DATA Z /1.0d0, 1.0d0, 1.75d0, 6.0d0/
C     MULTIPLY DELTA BY Z(1) IN INITIAL INTERVAL SEARCH
C     DELTA.LE.Z(2)*ABSDIF IN INITIAL INTERVAL SEARCH MEANS DONT GO ON
C     DELTA.LT.Z(3)*ABSDIF IF K.LT.6 AS FOR Z(2)
C     DELTA.GT.Z(4)*ABSDIF MEANS KMAX=5 AFTER INITIAL INTERVAL SEARCH
C
C     IN THE COMMENTS BELOW, F(K,I) REFERS TO THE FUNCTION VALUE
C     COMPUTED FOR THE I'TH NODE OF THE K'TH FORMULA.
C    P   Points
C Indices      Use  (C = Correction Coefficient N&W = Nodes and Weights)
C    1      3  C for F(1,1)
C   2-  3      N&W for F(2,1)
C   4-  5   7  C for F(1,1) and F(2,1)
C   6-  9      N&W for F(3, 1-2)
C  10- 13  15  C for F(1,1), F(2,1), F(3,1-2)
C  14- 21      N&W for F(4, 1-4)
C  22- 29  31  C for F(1,1), F(2,1), F(3,1-2), F(4,1-4).
C  30- 45      N&W FOR F(5,1-8).
C  46- 61  63  C for F(1,1), F(2,1), F(3,1-2), F(4,1-4), F(5,1-8).
C  62- 93      N&W FOR F(6,1-16).
C  94-105 127  C for F(3,1), F(4,1-2), F(5,1-3), F(6,1-6).
C 106-169      N&W FOR F(7,1-32).
C 170-177 256  C for F(4,1), F(5,1), F(6,1-2), F(7,1-4).
C 178:305      N&W FOR F(8,1-64).
C
C     CORRECTIONS, NODES, AND WEIGHTS FOR THE  3-POINT FORMULA.
C
c++ Save data by elements if ~.C.
      DATA P(1) / -.11111111111111111111D+00 /
      DATA P(2) / +.22540333075851662296D+00 /
      DATA P(3) / +.55555555555555555556D+00 /
      DATA P(4) / +.647209421402969791D-02 /
      DATA P(5) / -.928968790944433705D-02 /
      DATA P(6) / +.39508731291979716579D-01 /
      DATA P(7) / +.10465622602646726519D+00 /
      DATA P(8) / +.56575625065319744200D+00 /
      DATA P(9) / +.40139741477596222291D+00 /
      DATA P(10) /  +.5223046896961622D-04 /
      DATA P(11) /  +.17121030961750000D-03 /
      DATA P(12) /  -.724830016153892898D-03 /
      DATA P(13) /  -.7017801099209042D-04 /
      DATA P(14) /  +.61680367872449777899D-02 /
      DATA P(15) /  +.17001719629940260339D-01 /
      DATA P(16) /  +.11154076712774300110D+00 /
      DATA P(17) /  +.92927195315124537686D-01 /
      DATA P(18) /  +.37889705326277359705D+00 /
      DATA P(19) /  +.17151190913639138079D+00 /
      DATA P(20) /  +.77661331357103311837D+00 /
      DATA P(21) /  +.21915685840158749640D+00 /
      DATA P(22) /  +.682166534792D-08 /
      DATA P(23) /  +.12667409859336D-06 /
      DATA P(24) /  +.59565976367837165D-05 /
      DATA P(25) /  +.1392330106826D-07 /
      DATA P(26) /  -.6629407564902392D-04 /
      DATA P(27) /  -.704395804282302D-06 /
      DATA P(28) /  -.34518205339241D-07 /
      DATA P(29) /  -.814486910996D-08 /
      DATA P(30) /  +.90187503233240234038D-03 /
      DATA P(31) /  +.25447807915618744154D-02 /
      DATA P(32) /  +.18468850446259893130D-01 /
      DATA P(33) /  +.16446049854387810934D-01 /
      DATA P(34) /  +.70345142570259943330D-01 /
      DATA P(35) /  +.35957103307129322097D-01 /
      DATA P(36) /  +.16327406183113126449D+00 /
      DATA P(37) /  +.56979509494123357412D-01 /
      DATA P(38) /  +.29750379350847292139D+00 /
      DATA P(39) /  +.76879620499003531043D-01 /
      DATA P(40) /  +.46868025635562437602D+00 /
      DATA P(41) /  +.93627109981264473617D-01 /
      DATA P(42) /  +.66886460674202316691D+00 /
      DATA P(43) /  +.10566989358023480974D+00 /
      DATA P(44) /  +.88751105686681337425D+00 /
      DATA P(45) /  +.11195687302095345688D+00 /
      DATA P(46) /  +.371583D-15 /
      DATA P(47) /  +.21237877D-12 /
      DATA P(48) /  +.10522629388435D-08 /
      DATA P(49) /  +.1748029D-14 /
      DATA P(50) /  +.3475718983017160D-06 /
      DATA P(51) /  +.90312761725D-11 /
      DATA P(52) /  +.12558916D-13 /
      DATA P(53) /  +.54591D-15 /
      DATA P(54) /  -.72338395508691963D-05 /
      DATA P(55) /  -.169699579757977D-07 /
      DATA P(56) /  -.854363907155D-10 /
      DATA P(57) /  -.12281300930D-11 /
      DATA P(58) /  -.462334825D-13 /
      DATA P(59) /  -.42244055D-14 /
      DATA P(60) /  -.88501D-15 /
      DATA P(61) /  -.40904D-15 /
      DATA P(62) /  +.12711187964238806027D-03 /
      DATA P(63) /  +.36322148184553065969D-03 /
      DATA P(64) /  +.27937406277780409196D-02 /
      DATA P(65) /  +.25790497946856882724D-02 /
      DATA P(66) /  +.11315242452570520059D-01 /
      DATA P(67) /  +.61155068221172463397D-02 /
      DATA P(68) /  +.27817125251418203419D-01 /
      DATA P(69) /  +.10498246909621321898D-01 /
      DATA P(70) /  +.53657141626597094849D-01 /
      DATA P(71) /  +.15406750466559497802D-01 /
      DATA P(72) /  +.89628843042995707499D-01 /
      DATA P(73) /  +.20594233915912711149D-01 /
      DATA P(74) /  +.13609206180630952284D+00 /
      DATA P(75) /  +.25869679327214746911D-01 /
      DATA P(76) /  +.19305946804978238813D+00 /
      DATA P(77) /  +.31073551111687964880D-01 /
      DATA P(78) /  +.26024395564730524132D+00 /
      DATA P(79) /  +.36064432780782572640D-01 /
      DATA P(80) /  +.33709033997521940454D+00 /
      DATA P(81) /  +.40715510116944318934D-01 /
      DATA P(82) /  +.42280428994795418516D+00 /
      DATA P(83) /  +.44914531653632197414D-01 /
      DATA P(84) /  +.51638197305415897244D+00 /
      DATA P(85) /  +.48564330406673198716D-01 /
      DATA P(86) /  +.61664067580126965307D+00 /
      DATA P(87) /  +.51583253952048458777D-01 /
      DATA P(88) /  +.72225017797817568492D+00 /
      DATA P(89) /  +.53905499335266063927D-01 /
      DATA P(90) /  +.83176474844779253501D+00 /
      DATA P(91) /  +.55481404356559363988D-01 /
      DATA P(92) /  +.94365568695340721002D+00 /
      DATA P(93) /  +.56277699831254301273D-01 /
      DATA P(94) /  +.1041098D-15 /
      DATA P(95) /  +.249472054598D-10 /
      DATA P(96) /  +.55D-20 /
      DATA P(97) /  +.290412475995385D-07 /
      DATA P(98) /  +.367282126D-13 /
      DATA P(99) /  +.5568D-18 /
      DATA P(100) /  -.871176477376972025D-06 /
      DATA P(101) /  -.8147324267441D-09 /
      DATA P(102) /  -.8830920337D-12 /
      DATA P(103) /  -.18018239D-14 /
      DATA P(104) /  -.70528D-17 /
      DATA P(105) /  -.506D-19 /
      DATA P(106) /  +.17569645108401419961D-04 /
      DATA P(107) /  +.50536095207862517625D-04 /
      DATA P(108) /  +.40120032808931675009D-03 /
      DATA P(109) /  +.37774664632698466027D-03 /
      DATA P(110) /  +.16833646815926074696D-02 /
      DATA P(111) /  +.93836984854238150079D-03 /
      DATA P(112) /  +.42758953015928114900D-02 /
      DATA P(113) /  +.16811428654214699063D-02 /
      DATA P(114) /  +.85042788218938676006D-02 /
      DATA P(115) /  +.25687649437940203731D-02 /
      DATA P(116) /  +.14628500401479628890D-01 /
      DATA P(117) /  +.35728927835172996494D-02 /
      DATA P(118) /  +.22858485360294285840D-01 /
      DATA P(119) /  +.46710503721143217474D-02 /
      DATA P(120) /  +.33362148441583432910D-01 /
      DATA P(121) /  +.58434498758356395076D-02 /
      DATA P(122) /  +.46269993574238863589D-01 /
      DATA P(123) /  +.70724899954335554680D-02 /
      DATA P(124) /  +.61679602220407116350D-01 /
      DATA P(125) /  +.83428387539681577056D-02 /
      DATA P(126) /  +.79659974529987579270D-01 /
      DATA P(127) /  +.96411777297025366953D-02 /
      DATA P(128) /  +.10025510022305996335D+00 /
      DATA P(129) /  +.10955733387837901648D-01 /
      DATA P(130) /  +.12348658551529473026D+00 /
      DATA P(131) /  +.12275830560082770087D-01 /
      DATA P(132) /  +.14935550523164972024D+00 /
      DATA P(133) /  +.13591571009765546790D-01 /
      DATA P(134) /  +.17784374563501959262D+00 /
      DATA P(135) /  +.14893641664815182035D-01 /
      DATA P(136) /  +.20891506620015163857D+00 /
      DATA P(137) /  +.16173218729577719942D-01 /
      DATA P(138) /  +.24251603361948636206D+00 /
      DATA P(139) /  +.17421930159464173747D-01 /
      DATA P(140) /  +.27857691462990108452D+00 /
      DATA P(141) /  +.18631848256138790186D-01 /
      DATA P(142) /  +.31701256890892077191D+00 /
      DATA P(143) /  +.19795495048097499488D-01 /
      DATA P(144) /  +.35772335749024048622D+00 /
      DATA P(145) /  +.20905851445812023852D-01 /
      DATA P(146) /  +.40059606975775710702D+00 /
      DATA P(147) /  +.21956366305317824939D-01 /
      DATA P(148) /  +.44550486736806745112D+00 /
      DATA P(149) /  +.22940964229387748761D-01 /
      DATA P(150) /  +.49231224246628339785D+00 /
      DATA P(151) /  +.23854052106038540080D-01 /
      DATA P(152) /  +.54086998801016766712D+00 /
      DATA P(153) /  +.24690524744487676909D-01 /
      DATA P(154) /  +.59102017877011132759D+00 /
      DATA P(155) /  +.25445769965464765813D-01 /
      DATA P(156) /  +.64259616216846784762D+00 /
      DATA P(157) /  +.26115673376706097680D-01 /
      DATA P(158) /  +.69542355844328595666D+00 /
      DATA P(159) /  +.26696622927450359906D-01 /
      DATA P(160) /  +.74932126969651682339D+00 /
      DATA P(161) /  +.27185513229624791819D-01 /
      DATA P(162) /  +.80410249728889984607D+00 /
      DATA P(163) /  +.27579749566481873035D-01 /
      DATA P(164) /  +.85957576684743982540D+00 /
      DATA P(165) /  +.27877251476613701609D-01 /
      DATA P(166) /  +.91554595991628911629D+00 /
      DATA P(167) /  +.28076455793817246607D-01 /
      DATA P(168) /  +.97181535105025430566D+00 /
      DATA P(169) /  +.28176319033016602131D-01 /
      DATA P(170) /  +.3326D-18 /
      DATA P(171) /  +.114094770478D-11 /
      DATA P(172) /  +.2952436056970351D-08 /
      DATA P(173) /  +.51608328D-15 /
      DATA P(174) /  -.110177219650597323D-06 /
      DATA P(175) /  -.58656987416475D-10 /
      DATA P(176) /  -.23340340645D-13 /
      DATA P(177) /  -.1248950D-16 /
      DATA P(178) /  +.24036202515353807630D-05 /
      DATA P(179) /  +.69379364324108267170D-05 /
      DATA P(180) /  +.56003792945624240417D-04 /
      DATA P(181) /  +.53275293669780613125D-04 /
      DATA P(182) /  +.23950907556795267013D-03 /
      DATA P(183) /  +.13575491094922871973D-03 /
      DATA P(184) /  +.61966197497641806982D-03 /
      DATA P(185) /  +.24921240048299729402D-03 /
      DATA P(186) /  +.12543855319048853002D-02 /
      DATA P(187) /  +.38974528447328229322D-03 /
      DATA P(188) /  +.21946455040427254399D-02 /
      DATA P(189) /  +.55429531493037471492D-03 /
      DATA P(190) /  +.34858540851097261500D-02 /
      DATA P(191) /  +.74028280424450333046D-03 /
      DATA P(192) /  +.51684971993789994803D-02 /
      DATA P(193) /  +.94536151685852538246D-03 /
      DATA P(194) /  +.72786557172113846706D-02 /
      DATA P(195) /  +.11674841174299594077D-02 /
      DATA P(196) /  +.98486295992298408193D-02 /
      DATA P(197) /  +.14049079956551446427D-02 /
      DATA P(198) /  +.12907472045965932809D-01 /
      DATA P(199) /  +.16561127281544526052D-02 /
      DATA P(200) /  +.16481342421367271240D-01 /
      DATA P(201) /  +.19197129710138724125D-02 /
      DATA P(202) /  +.20593718329137316189D-01 /
      DATA P(203) /  +.21944069253638388388D-02 /
      DATA P(204) /  +.25265540247597332240D-01 /
      DATA P(205) /  +.24789582266575679307D-02 /
      DATA P(206) /  +.30515340497540768229D-01 /
      DATA P(207) /  +.27721957645934509940D-02 /
      DATA P(208) /  +.36359378430187867480D-01 /
      DATA P(209) /  +.30730184347025783234D-02 /
      DATA P(210) /  +.42811783890139037259D-01 /
      DATA P(211) /  +.33803979910869203823D-02 /
      DATA P(212) /  +.49884702478705123440D-01 /
      DATA P(213) /  +.36933779170256508183D-02 /
      DATA P(214) /  +.57588434808916940190D-01 /
      DATA P(215) /  +.40110687240750233989D-02 /
      DATA P(216) /  +.65931563842274211999D-01 /
      DATA P(217) /  +.43326409680929828545D-02 /
      DATA P(218) /  +.74921067092924347640D-01 /
      DATA P(219) /  +.46573172997568547773D-02 /
      DATA P(220) /  +.84562412844234959360D-01 /
      DATA P(221) /  +.49843645647655386012D-02 /
      DATA P(222) /  +.94859641186738404810D-01 /
      DATA P(223) /  +.53130866051870565663D-02 /
      DATA P(224) /  +.10581543166444097714D+00 /
      DATA P(225) /  +.56428181013844441585D-02 /
      DATA P(226) /  +.11743115975265809315D+00 /
      DATA P(227) /  +.59729195655081658049D-02 /
      DATA P(228) /  +.12970694445188609414D+00 /
      DATA P(229) /  +.63027734490857587172D-02 /
      DATA P(230) /  +.14264168911376784347D+00 /
      DATA P(231) /  +.66317812429018878941D-02 /
      DATA P(232) /  +.15623311732729139895D+00 /
      DATA P(233) /  +.69593614093904229394D-02 /
      DATA P(234) /  +.17047780536259859981D+00 /
      DATA P(235) /  +.72849479805538070639D-02 /
      DATA P(236) /  +.18537121234486258656D+00 /
      DATA P(237) /  +.76079896657190565832D-02 /
      DATA P(238) /  +.20090770903915859819D+00 /
      DATA P(239) /  +.79279493342948491103D-02 /
      DATA P(240) /  +.21708060588171698360D+00 /
      DATA P(241) /  +.82443037630328680306D-02 /
      DATA P(242) /  +.23388218069623990928D+00 /
      DATA P(243) /  +.85565435613076896192D-02 /
      DATA P(244) /  +.25130370638306339718D+00 /
      DATA P(245) /  +.88641732094824942641D-02 /
      DATA P(246) /  +.26933547875781873867D+00 /
      DATA P(247) /  +.91667111635607884067D-02 /
      DATA P(248) /  +.28796684463774796540D+00 /
      DATA P(249) /  +.94636899938300652943D-02 /
      DATA P(250) /  +.30718623022088529711D+00 /
      DATA P(251) /  +.97546565363174114611D-02 /
      DATA P(252) /  +.32698116976958152079D+00 /
      DATA P(253) /  +.10039172044056840798D-01 /
      DATA P(254) /  +.34733833458998250389D+00 /
      DATA P(255) /  +.10316812330947621682D-01 /
      DATA P(256) /  +.36824356228880576959D+00 /
      DATA P(257) /  +.10587167904885197931D-01 /
      DATA P(258) /  +.38968188628481359983D+00 /
      DATA P(259) /  +.10849844089337314099D-01 /
      DATA P(260) /  +.41163756555233745857D+00 /
      DATA P(261) /  +.11104461134006926537D-01 /
      DATA P(262) /  +.43409411457634557737D+00 /
      DATA P(263) /  +.11350654315980596602D-01 /
      DATA P(264) /  +.45703433350168850951D+00 /
      DATA P(265) /  +.11588074033043952568D-01 /
      DATA P(266) /  +.48044033846254297801D+00 /
      DATA P(267) /  +.11816385890830235763D-01 /
      DATA P(268) /  +.50429359208123853983D+00 /
      DATA P(269) /  +.12035270785279562630D-01 /
      DATA P(270) /  +.52857493412834112307D+00 /
      DATA P(271) /  +.12244424981611985899D-01 /
      DATA P(272) /  +.55326461233797152625D+00 /
      DATA P(273) /  +.12443560190714035263D-01 /
      DATA P(274) /  +.57834231337383669993D+00 /
      DATA P(275) /  +.12632403643542078765D-01 /
      DATA P(276) /  +.60378719394238406082D+00 /
      DATA P(277) /  +.12810698163877361967D-01 /
      DATA P(278) /  +.62957791204992176986D+00 /
      DATA P(279) /  +.12978202239537399286D-01 /
      DATA P(280) /  +.65569265840056197721D+00 /
      DATA P(281) /  +.13134690091960152836D-01 /
      DATA P(282) /  +.68210918793152331682D+00 /
      DATA P(283) /  +.13279951743930530650D-01 /
      DATA P(284) /  +.70880485148175331803D+00 /
      DATA P(285) /  +.13413793085110098513D-01 /
      DATA P(286) /  +.73575662758907323806D+00 /
      DATA P(287) /  +.13536035934956213614D-01 /
      DATA P(288) /  +.76294115441017027278D+00 /
      DATA P(289) /  +.13646518102571291428D-01 /
      DATA P(290) /  +.79033476175681880523D+00 /
      DATA P(291) /  +.13745093443001896632D-01 /
      DATA P(292) /  +.81791350324074780175D+00 /
      DATA P(293) /  +.13831631909506428676D-01 /
      DATA P(294) /  +.84565318851862189130D+00 /
      DATA P(295) /  +.13906019601325461264D-01 /
      DATA P(296) /  +.87352941562769803314D+00 /
      DATA P(297) /  +.13968158806516938516D-01 /
      DATA P(298) /  +.90151760340188079791D+00 /
      DATA P(299) /  +.14017968039456608810D-01 /
      DATA P(300) /  +.92959302395714482093D+00 /
      DATA P(301) /  +.14055382072649964277D-01 /
      DATA P(302) /  +.95773083523463639678D+00 /
      DATA P(303) /  +.14080351962553661325D-01 /
      DATA P(304) /  +.98590611358921753738D+00 /
      DATA P(305) /  +.14092845069160408355D-01 /
C
C     *****    STATEMENT FUNCTIONS     *********************************
C
C TDECR   IS USED TO TRANSFORM AN ABSCISSA FROM THE CURRENT COORDINATE
C         SYSTEM TO ONE IN WHICH NSUB IS DECREMENTED BY A FACTOR OF 2.
C     TDECR(ZL1)=TA+(ZL1-TA)*((ZL1-TA)/TB)
      TDECR(ZL1)=TA*(1.0d0+TA/TB)+ZL1*((ZL1-TA)/TB-TA/TB)
C TINCR   IS USED TO TRANSFORM AN ABSCISSA FROM THE CURRENT COORDINATE
C         SYSTEM TO ONE IN WHICH NSUB IS INCREMENTED BY A FACTOR OF 2.
      TINCR(ZL1)=TA+SIGN(SQRT(ABS(TB*(ZL1-TA))),TB)
C
C     *****    PROCEDURES     ******************************************
C
      IF (WHERE.NE.0) GO TO 2900
C
C     CHECK FOR TRIVIAL CASE
C
      PART=1
      RELEPS=0.0d0
10    ERRT(1)=0.0d0
c--D Next line special: S => D, X => Q
      RESULT(1)=0.0D0
      XCDOBT(1)=.FALSE.
      EDUE2A=0.0d0
      EDUE2B=0.0d0
      IF (AINIT.EQ.BINIT) GO TO 2320
C
C     SETUP FOR COMPLETE INTERVAL.
C     SET INITIAL VALUES.
C
      IFLAG(1)=0
      XJ = 0.D0
      XJP = 0.D0
      FAT(1) = 0.D0
      FAT(2) = 0.D0
      ERRC = 0.0d0
      ERRT(2)=0.0d0
      START(1)=AINIT
      END(1)=BINIT
      EPSR=MAX(EPSO,ABS(RELEPS))
      EPSMAX=XEPS*EPSR
      NSUB=0
      TA=START(1)
      TLEN=ABS(BINIT-AINIT)
      STEP(1)=1.0d0
      WORRY(1)=END(1)
      INIT=.TRUE.
      DID1=.FALSE.
      HAVDIF=.FALSE.
      EPS=EPSR
      IF (TALOC.NE.0) THEN
         I=ABS(TALOC)
         TA=WORK(I)
         NSUB=2
         IF (TALOC.LT.0) NSUB=4
         IF (SIGN(1.0d0,TA-START(1))*(END(1)-TA).GT.0.0d0) THEN
            IF (TA.NE.START(1)) THEN
C              TA IS INSIDE THE INTERVAL.  SUBDIVIDE IMMEDIATELY.
               PART=2
               STEP(2)=-1.0d0
               END(2)=START(1)
               WORRY(2)=END(2)
c--D Next line special: S => D, X => Q
               RESULT(2)=0.0d0
               XCDOBT(2)=.FALSE.
               START(1)=TA
               START(2)=TA
               TB=END(2)-TA
               STEP(1)=DINTSM(TA)
               NSUBSV=NSUB
               TASAVE=TA
               EPS=EPSR*ABS(END(2)-START(2))/TLEN
            END IF
         ELSE
C           TA IS OUTSIDE THE INTERVAL.
            IF (ABS(TA-START(1)).GE.ABS(TA-END(1))) THEN
C              TA IS NEARER BINIT.  TURN THE INTERVAL AROUND.
               X=END(1)
               END(1)=START(1)
               START(1)=X
               STEP(1)=-STEP(1)
            END IF
            IF (TA.NE.START(1)) THEN
               TB=END(1)-TA
               START(1)=TINCR(START(1))
               IF (NSUB.EQ.4) START(1)=TINCR(START(1))
            END IF
         END IF
      END IF
40    RNDC=FER
      EPNOIZ=ENINF
      FATAS=.FALSE.
      DISCF=0
      DISCHK=0
      FSAVED=.FALSE.
      KMAX=KMAXF
      KAIMT=MIN(MAX(5,KMAXF),MAXK-2)
      ALOCAL=START(PART)
      TB=END(PART)-TA
C     GO COMPUTE MINIMUM STEPSIZE
50    DELMIN=DINTSM(ALOCAL)
60    BLOCAL=END(PART)
      TEND=BLOCAL
      IEND=.TRUE.
70    FATBS=.FALSE.
      KMIN=0
      GO TO 80
75    EPS=EPSR*DELTA/TLEN
      INIT=.FALSE.
C
C     INITIALIZE VALUES NECESSARY TO INTEGRATE ONE PANEL.
C
80    LOCAL(3)=ALOCAL
      LOCAL(4)=BLOCAL
      IF (NSUB.NE.0) THEN
         LOCAL(3)=TDECR(LOCAL(3))
         LOCAL(4)=TDECR(LOCAL(4))
         IF (NSUB.NE.2) THEN
            LOCAL(3)=TDECR(LOCAL(3))
            LOCAL(4)=TDECR(LOCAL(4))
         END IF
      END IF
      DELTA=BLOCAL-ALOCAL
      DIFF=0.5d0*DELTA
      ABSDIF=ABS(DIFF)
      PX=DIFF
      DELTA=ABS(DELTA)
      SUM=ALOCAL+0.5d0*(BLOCAL-ALOCAL)
      EPSS=MAX(EPS,EPSMAX)
      DID1=START(PART).NE.ALOCAL
      ROUNDF=.FALSE.
      RE=0.2d0
      REPROD=1.0d0
      IP=1
      IOLD=0
      K=2
      IF (IPRINT.GT.1) CALL DINTO (1,WORK)
C
C     APPLY 1-POINT GAUSS FORMULA (MIDPOINT RULE).
C
      ABSCIS=SUM
      WHERE=1
      GO TO 2850
90    FUNCT(1)=FNCVAL
      PF1=FUNCT(1)
      PF2=FUNCT(1)
      PACUM=FUNCT(1)*DELTA
      PAACUM=ABS(PACUM)
      ERRI=PACUM
      ERR=ERRI
      PEPSMN=ABS(PACUM*RNDC)
      IF (FEA.NE.0) PEPSMN=PEPSMN+ABS(DELTA*T*WORK(FEA))
C
C     GO ON TO NEXT FORMULA.
C
      ACUM=FUNCT(1)*P(1)
      AACUM=ABS(ACUM)
      EPSMIN=0.0d0
      GO TO 130
95    DISCHK=MAX(DISCHK,1)
100   K=K+1
      PACUM=ACUM
      PEPSMN=EPSMIN
      EPSMIN=ENZER
c--D Next line special: S => D, X => Q
      ACUM=0.0D0
      AACUM=0.0D0
C
C     COMPUTE CONTRIBUTION FROM PREVIOUSLY CALCULATED FUNCTION VALUES.
C
      KK=KORECT(K-2)
110      J=MOD(KK,400)
         KK=KK/400
         J1=MOD(J,20)
         J2=J/20
         DO 120 J=J1,J2
            IP=IP+1
            ACUM=ACUM+P(IP)*FUNCT(J)
120      CONTINUE
         IF (KK.NE.0) GO TO 110
C
C     COMPUTE CONTRIBUTION FROM NEW FUNCTION VALUES.
C
130   INEW=IOLD+1
      IOLD=IOLD+INEW
      J2=MOD(FSTORE(K-1),100)
      J1=FSTORE(K-1)/100
      J=INEW
140   X=P(IP+1)*DIFF
      X1=BLOCAL-X
      X2=ALOCAL+X
      ABSCIS=X1
      ERRF=0.0d0
      WHERE=2
      IF (J.EQ.INEW) GO TO 2840
C     THE FIRST FUNCTION VALUE REQUESTED FOR EACH FORMULA IS THE
C     ONE NEAREST THE EDGE OF THE PANEL.
      GO TO 2850
150   F1=FNCVAL
      IF (FEA.NE.0) ERRF=ABS(T*WORK(FEA))
      ABSCIS=X2
      WHERE=3
      IF (J.EQ.INEW) GO TO 2840
      GO TO 2850
160   F2=FNCVAL
      IF (FEA.NE.0) ERRF=ERRF+ABS(T*WORK(FEA))
C     UPDATE TOLERANCE IF A FUNCTION VALUE WAS REQUESTED AT THE EDGE OF
C     THE INITIAL INTERVAL.
      IF (J.EQ.INEW) THEN
         NEWEPS=.FALSE.
         IF (ALOCAL.EQ.AINIT) THEN
            EDUE2A=ABS(F1*ERRINA)
            NEWEPS=.TRUE.
         ELSE IF (ALOCAL.EQ.BINIT) THEN
            EDUE2B=ABS(F1*ERRINB)
            NEWEPS=.TRUE.
         END IF
         IF (BLOCAL.EQ.AINIT) THEN
            EDUE2A=ABS(F2*ERRINA)
            NEWEPS=.TRUE.
         ELSE IF (BLOCAL.EQ.BINIT) THEN
            EDUE2B=ABS(F2*ERRINB)
            NEWEPS=.TRUE.
         END IF
         IF (NEWEPS) THEN
            EPS=EDUE2A+EDUE2B
            EPS=MAX(EPSR-EPS,0.1d0*(ERRT(1)+ERRT(2)+EPS))*DELTA/TLEN
         END IF
      END IF
C     COMPUTE EPSMIN
      X1=(ABS(X1*(F1-PF1))+ABS(X2*(F2-PF2)))*EMEPS
      X1=X1/ABS(X-PX)
      PX=X
      PF1=F1
      PF2=F2
      IP=IP+2
      S=ABS(F2)
      IF (NSUB.NE.0) S=MAX(S,ABS(ANSWER*(TA/(X2-TA))))
      EPSMIN=EPSMIN+P(IP)*(X1+(ABS(F1)+S)*RNDC+ERRF)
      AACUM=AACUM+P(IP)*(ABS(F1)+ABS(F2))
      F2=F2+F1
      ACUM=ACUM+P(IP)*F2
      IF (J1.LE.J2) THEN
         FUNCT(J1)=F2
         J1=J1+1
      END IF
      IF (J.LE.7) FUNCT(17+J)=F1
      J=J+1
      IF (J.LE.IOLD) GO TO 140
c--D Next line special: S => D, X => Q
      ACUM=ABSDIF*ACUM+0.5D0*PACUM
      AACUM=MAX(EDELM1,ABSDIF*AACUM+0.5D0*PAACUM)
      EPSMIN=MAX(EDELM1,ABSDIF*EPSMIN+0.5d0*PEPSMN)
      EPSR=MAX(EPSR,EPSMIN)
C
C     CHECK FOR CONVERGENCE.
C
      EP=MAX(ABS(ERRI),0.125d0*EPSMIN)
      ERRI=ACUM-PACUM
      ERR=ABS(ERRI)
      EXTRA=0.0d0
      REP=RE
      IF (INIT) THEN
         IF (RELTOL.NE.0.0d0) THEN
            IF (RELEPS.LE.0.0d0) THEN
               RELEPS=ABS(ACUM)-ERR
               RELEPS=-ABS(WORK(RELTOL+2))*MAX(RELEPS,0.0d0)
               EPSR=MAX(EPS,ABS(RELEPS))
               EPSMAX=XEPS*EPSR
               EPSS=MAX(EPSMAX,EPSR)
            END IF
         END IF
      END IF
      EPS=MAX(EPSS,(EPSMIN/AACUM)**RELOBT*AACUM)
      RE=ERR/EP
      TPS=RE+RE
      IF (NSUB.NE.0) THEN
        TPS=0.25d0-1000.0d0*ABS((TA-LOCAL(3))/(LOCAL(4)-LOCAL(3)))**NSUB
        TPS=2.0d0*MAX(RE,RECORR(NSUB-1)*TPS)
      END IF
      REPROD=MIN(RE*REPROD,1.0d0)
      FAIL=.TRUE.
      SEARCH=1
      WHERE=2
      HAVDIF=HAVDIF.AND.K.GT.4
      J=MIN(0,K-KAIMT)
      IF (J.LT.0) TPS=MAX(TPS,2.0d0*MIN(5.0d0*TPS,MIN(REP,0.125d0)))
      TP=EPSMIN*(FUDGE(J+5)/RNDC)
      IF (ERR.GT.TP) THEN
C       IF ERROR IS LARGE RELATIVE TO THE OBTAINABLE PRECISION, INCREASE
C       THE ERROR ESTIMATE.  (THIS TEST IS CONSERVATIVE WHEN K IS SMALL)
         TP1=0.25d0/MAX(1.0d0,25.0d0*REP*REP)
190         ERR=4.0d0*ERR
            TP=32.0d0*TP
            TPS=TPS*(1.0d0+TPS)/(TP1+TPS)
            IF (ERR.GT.TP) GO TO 190
C        END OF INCREASING ERROR ESTIMATE WHEN RELATIVE ERROR IS LARGE.
         IF ((ERR*RNDC).GE.EPSMIN) THEN
            IF (ERR.GT.EPS) IF (DISCHK) 220,220,230
C           GOT ACCURACY REQUESTED, INTEGRAND ASSUMED TO BE ALL NOISE.
C           DO NOT BLINDLY ACCEPT RESULT IF ON INITIAL INTERVAL.
            IF (K-2) 220,220,230
         END IF
      END IF
      ERR=ERR*TPS
C     DO NOT ACCEPT RESULT WITHOUT LOOKING AT DIFFERENCE LINE IF
C     WORKING ON THE INITIAL INTERVAL OR THE LEFT BOUNDARY WAS A
C     STOP POINT OR AN END-POINT SINGULARITY, OR IF RE*REP .GT. RELOOK.
      IF (DISCHK.GT.0) GO TO 230
      IF (ERR.GT.EPS) GO TO 220
      IF (K.NE.KAIMT) GO TO 230
      IF (DISCHK.LT.0) GO TO 230
      IF (RE*REP.GT.RELOOK) GO TO 230
      IF (INIT) GO TO 230
      IF (RE.GE.1.0d0) GO TO 220
      WHERE=3
      GO TO 230
220   WHERE=1
230   IF (IPRINT.GT.2) CALL DINTO (2,WORK)
      EP=ERRC
      ERRC=ERR
      ERRI=ABS(ERRI)
      GO TO (240,250,2250), WHERE
240   IF (K.LT.MIN(KMAX,KAIMT)) GO TO 100
C
C     THE ERROR IS NOT SMALL ENOUGH, OR IT SEEMS SMALL ENOUGH BUT WE
C     HAVE A REASON TO BE SUSPICIOUS.
C
250   IF (HAVDIF) GO TO 500
C     COMPUTE ABSCISSAE AND UNSCRAMBLE FUNCTION VALUES INTO
C     USABLE ORDER.
      FT(9)=FUNCT(1)
      XT(9)=SUM
      J=1
      L=8
      KK=MIN(K,4)
      DO 270 J1=2,KK
         INC=L
         L=L/2
         DO 260 I=L,7,INC
            J=J+1
            J2=IH(I)
            BETA=DIFF*P(J2)
            XT(I+1)=ALOCAL+BETA
            XT(17-I)=BLOCAL-BETA
            FT(I+1)=FUNCT(J)-FUNCT(16+J)
            FT(17-I)=FUNCT(16+J)
260      CONTINUE
270   CONTINUE
      ENDPTS=1
      IF (K.NE.2) GO TO 300
      ENDPTS=3
C     COMPUTE F AT (OR NEAR) THE ENDS OF THE PANEL.
      WHERE=4
      WHERE2=1
      GO TO 2810
280   FT(1)=FATA
      XT(1)=X
      WHERE2=2
      GO TO 2810
290   FT(17)=FATB
      XT(17)=X
300   WHERE=0
C     WHERE=0 TELLS DINTDL TO FORM AND ANALYZE DIFFERENCE LINES.
      GO TO 370
C     A NEW FUNCTION VALUE HAS BEEN COMPUTED.
310   WHERE=3
C     NSUB HAS BEEN CHANGED OR THE DIFFERENCE LINES HAVE BEEN OTHERWISE
C     MUTILATED.  UPDATE THE DIFFERENCE LINES.
350   IF (KDIM.EQ.1) ERRF=0.0d0
      IF (FEA.NE.0) ERRF=ABS(T*WORK(FEA))
C     WHERE TELLS DINTDU WHETHER TO ADD A FUNCTION VALUE IN THE
C     MIDDLE, ON THE ALOCAL END OR THE BLOCAL END, OR TO RECOMPUTE THE
C     DIFFERENCE LINES.
      CALL DINTDU
360   WHERE=1
C     WHERE .NE. 0 TELLS DINTDL TO ANALYZE DIFFERENCE LINES THAT
C     HAVE ALREADY BEEN FORMED.
370   CALL DINTDL (WORK)
500   KONVRG=IC1+LENDT-IC2
C
C     EXAMINE THE ISTOP ARRAY.
C
510   JUMPS=1
      IF (ISTOP(2,2)+18-ISTOP(2,1).EQ.0) GO TO 530
      JUMPS=4
      IF (ISTOP(2,2)*(18-ISTOP(2,1)).NE.0) GO TO 530
      JUMPS=3
      IF ((LENDT-ISTOP(2,1))*(1-ISTOP(2,2)).EQ.0) JUMPS=2
      GO TO 530
C     IF WE WERE NOT DOING A JUMP SEARCH AND FOUND JUMPS, START ONE.
520   ESOLD=-1.0d0
      COUNT=0.0d0
      XJUMP=XJUMPS
      SEARCH=5
C     SET A CAUTION POINT.
      WORRY(PART)=XT(LENDT)
      IF (PART.EQ.2.AND.STEP(1)*STEP(2).LT.0.0d0) WORRY(1)=XT(1)
530   CONTINUE
C     DO BLOCK
         J1=ISTOP(2,1)
         J2=ISTOP(2,2)
         ERR=ERRC
         GO TO (540,2150,770,860,800,535), SEARCH
c        One hopes SEARCH==6 and JUMPS==4 can't happen
535      GO TO (2160,620,970,520), JUMPS
540      COUNT=0.0d0
         IF (JUMPS.EQ.1) GO TO 1920
         IF (K.LE.2) GO TO 95
         IF (DELTA.LT.DELMIN+DELMIN) GO TO 2170
         IF (DISCHK.EQ.0) DISCHK=-1
         IF (JUMPS.EQ.3) GO TO 1920
         IF (JUMPS.NE.2) GO TO 790
C   88888888888
C     END POINT JUMP.
C
         IF ((TA-ALOCAL)*SIGN(1.0d0,BLOCAL-TA).GT.0.0d0)
     *   IF (ENDPTS-2) 2750,2760,560
560      CONTINUE
         IF (NSUB.NE.0) THEN
C     NSUB CAN ONLY BE NON-ZERO WHEN .NOT.DID1 IF THE JUMP WAS EXPECTED.
            IF (.NOT.DID1) XJ=ABS(XJ)
         END IF
         IF (XJ.GE.0.0d0) THEN
            IF (XJP.GT.0.0d0) THEN
               IF (J2.EQ.0) GO TO 630
               IF (RE.GT.1.0d0) GO TO 670
               ERR=ERR+EXTRA
               IF (ERR.LE.EPS) GO TO 750
c??            Do we need to ask (re .lt. rep) and compute two ways?
c##            S=RE*(RE/REP)*ERR
c##            IF (S.GT.EPS) GO TO 670
               IF (RE*RE*ERR.GT.REP*EPS) GO TO 670
C         IF KAIMT IS EVER PERMITTED .GT. MAXK-2, MAXK-K MUST BE TESTED.
               IF (K-KAIMT-1) 2740,2740,670
            END IF
         END IF
         IF (ENDPTS-2) 2750,2760,600
600      CONTINUE
         IF (SEARCH.EQ.1) THEN
C           IF THE ERROR IN THE FUNCTION AT THE BOUNDARY OF THE PANEL IS
C           LARGER THAN THE DIFFERENCE BETWEEN THE FUNCTION AND THE
C           VALUE PREDICTED BY EXTRAPOLATING THE DIFFERENCE LINE,
C           REPLACE THE FUNCTION BY THE EXTRAPOLATED VALUE, AND DECREASE
C           THE ERROR BY THE CORRECTION.  THIS DOESN'T AFFECT THE ERROR
C           IN THE INTEGRAL BECAUSE FUNCTION VALUES AT THE ENDS DON'T
C           CONTRIBUTE TO THE INTEGRAL.
            ENDPHI(2)=LENDT
            DO 610 WHERE2 = 1, 2
               IF (FATS(WHERE2)) THEN
           IF (ERRAT(WHERE2).GT.ABS(PHIS(ENDPHI(3-WHERE2),WHERE2))) THEN
               FAT(WHERE2)=FAT(WHERE2)-PHIS(ENDPHI(3-WHERE2),WHERE2)
               FT(ENDPHI(WHERE2))=FAT(WHERE2)
          ERRAT(WHERE2)=ERRAT(WHERE2)-ABS(PHIS(ENDPHI(3-WHERE2),WHERE2))
                     GO TO 310
                  END IF
               END IF
610      CONTINUE
         END IF
         DISCHK=MAX(DISCHK,1)
         IF (RE.LE.1.0d0) THEN
            IF (J2.EQ.0) GO TO 630
            IF (ERR.LE.EPS) GO TO 760
         END IF
         IF (K.LT.4) GO TO 100
C     END BLOCK
620   IF (J2.NE.0) GO TO 670
C
C     IF WE SEE A JUMP NEAR BLOCAL AND WE ARE AT THE END OF A PART,
C     TURN THE INTERVAL AROUND.  IF WE ARE NOT AT THE END OF A
C     PART, TAKE HALF THE CURRENT INTERVAL AND START OVER.
C
      IF (SEARCH.NE.1) GO TO 1440
630   CONTINUE
C     DO BLOCK
         IF (IEND) THEN
            IF (NSUB.EQ.0) GO TO 640
            IF (ABS(ALOCAL-TA).GE.0.125d0*ABS(TB)) THEN
               IF (.NOT.DID1.OR.XJP.LE.0.1d0) GO TO 640
               ISTOP(2,1)=18
C              IF THE JUMP IS WEAK, PRETEND ITS NOT THERE.
               GO TO 510
            END IF
         END IF
         DELTA=0.5d0*DELTA
         GO TO 2610
C     END BLOCK
640   CONTINUE
      if (NSUB.NE.0) then
         CALL DINTNS (NSRB)
         GO TO 640
      END IF
      START(PART)=BLOCAL
      BLOCAL=ALOCAL
      END(PART)=ALOCAL
      ALOCAL=START(PART)
      TA=ALOCAL
      TB=END(PART)-TA
      IF (DISCHK.LE.0) DISCHK=-2
      CALL DINTNS (NSIA)
C     TAKE HALF OF THE INTERVAL TO AVOID GETTING IN A LOOP.
      DELTA=0.5d0*DELTA
      TEND=ALOCAL+SIGN(DELTA,TB)
      STEP(PART)=-STEP(PART)
      GO TO 2610
C
C     IF WE SEE AN EXPECTED AND BELIEVABLE JUMP NEAR ALOCAL, WE
C     SHOULD INCREASE NSUB.  OTHERWISE, GO INTO THE INTERVAL SEARCH.
C
670   CONTINUE
      IF (XJ.GT.0.0d0.AND.NSUB.LT.NSUBMX.AND.(INIT.OR.LENDT-IC2.GE.2))
     *THEN
         IF (RE.GE.1.0d0.AND.DID1.AND.DISCHK.NE.-2) THEN
            IF (K.LT.KAIMT) GO TO 100
         END IF
         IF ((TA-ALOCAL)*(BLOCAL-TA).GT.0.0d0) THEN
            TA=ALOCAL
            TB=END(PART)-TA
         END IF
         KAIMT=MIN(MAX(5,KMAXF),MAXK-2)
         KMAX=KMAXF
         IF (SEARCH.NE.1) THEN
            BLOCAL=XT(LENDT)
            IEND=.FALSE.
         END IF
         CALL DINTNS (NSIB)
         IF (.NOT.IEND.AND.(PART.NE.2.OR.DID1)) THEN
C           SET A CAUTION POINT.
       IF (ABS(BLOCAL-START(PART)).LT.ABS(TEND-START(PART))) TEND=BLOCAL
         END IF
         IF (SEARCH.EQ.1) THEN
            IF (IC2.LT.LENDT-2.OR.K.GT.3) THEN
               DELMIN=DINTSM(ALOCAL)
               EPS=EPSS
               GO TO 70
            END IF
            BLOCAL=ALOCAL+0.5d0*(BLOCAL-ALOCAL)
            DELTA=0.5d0*DELTA
         END IF
         IEND=.FALSE.
         DELMIN=DINTSM(ALOCAL)
         GO TO 75
      END IF
C     WE CAN NOT INCREASE NSUB.  START A SEARCH FOR THE CORRECT
C     INTERVAL LENGTH.
      IF (ENDPTS-2) 2750,2760,2155
C
C     ADD POINTS AT THE END NEAR THE JUMP UNTIL WE ARE SURE THE
C     JUMP IS NOT INSIDE THE INTERVAL.  THIS PROCESS CONTINUES UNTIL
C     THE JUMP VANISHES, MOVES INSIDE THE INTERVAL, OR THE ERROR
C     ESTIMATED ACROSS THE JUMP IS SUFFICIENTLY SMALL.  IF ALOCAL = TA
C     (WHICH MEANS WE ARE FAIRLY CERTAIN THE SINGULARITY IS NOT IN THE
C     INTERVAL), NSUB IS NOT ZERO, AND THE JUMP IS EXPECTED, ACCEPT THE
C     ANSWER WITHOUT DOING A SEARCH.
C
750   IF (ENDPTS-2) 2750,2760,760
760   SEARCH=-3
770   GO TO (1190,780,980,870), JUMPS
780   IF (NSUB.EQ.0) GO TO 980
      IF (XJ.LE.0.0d0) GO TO 980
      IF (ALOCAL-TA) 980,2240,980
C
C     JUMPS WERE NOTICED IN THE ORIGINAL DIFFERENCE LINE ON THIS PANEL.
C
790   IF (K+JUMPS.LT.6) GO TO 100
      IF (MAX(J1,J2).LE.NSUB+2) GO TO 2130
      EPSS=EPS
      IF (NSUB.NE.0) GO TO 810
      SEARCH=-4
      XJUMP=XJUMPS
      GO TO 940
800   IF (NSUB.EQ.0) GO TO 850
C     IF WE HAVE NOT DONE ANYTHING, REDUCE DELTA INSTEAD OF SEARCHING.
C     THERE WAS A REASON FOR SETTING NSUB NON-ZERO.
810   IF (J2.LT.2) J2=J1
      J=MIN(J1,J2)
      IF (ABS(2.0d0*(XT(1)-TA)).LE.ABS(XT(J)-TA)) THEN
         WORRY(PART)=XT(J)
         DELTA=0.75d0*ABS(XT(J)-ALOCAL)
         DISCHK=2
         GO TO 2610
      END IF
C     REDUCE NSUB TO ZERO BEFORE CONTINUING THE SEARCH.
      X1=DBLE(NSUB)
      J=1
      DO 820 I=1,LENDT
         T=(XT(I)-TA)/TB
         IF (NSUB.NE.2) T=T**(NSUB-1)
         FT(I)=FT(I)/(X1*T)
         XT(I)=TA+T*(XT(I)-TA)
         IF (I.NE.1) THEN
            IF (XT(I).NE.XT(I-1)) THEN
               J=J+1
               XT(J)=XT(I)
               FT(J)=FT(I)
            END IF
         END IF
820   CONTINUE
      LENDT=J
      SEARCH=6
C     FINISH DECREASING NSUB
840   IF (NSUB.NE.0) THEN
         CALL DINTNS (NSRB)
         GO TO 840
      END IF
      GO TO 310
C     ARE JUMPS STILL WITHIN RANGE OF PREVIOUS JUMPS?
850   SEARCH=4
860   IF (JUMPS.LT.4) GO TO 870
      IF (J1.LT.MIN(J1OLD,J2OLD)) GO TO 870
      IF (J2.LE.MAX(J1OLD,J2OLD)) GO TO 940
870   IF (NFEVAL.LE.NFJUMP+3) GO TO 2690
      GO TO 310
C
C     THE DIFFERENCE LINE HAS BEEN BALANCED.  CONTINUE THE JUMP SEARCH
C     OR COMPUTE THE INITIAL INTERVAL.
C
910   CONTINUE
      IF (SEARCH.NE.4) THEN
         IF (JUMPS-3) 920,930,915
915      IF (J2-(NSUB+2)) 930,520,520
920      I=6-NSUB/2
         J=I
         IF (IC1.GT.J) J=MIN(I+1,6)
         KK=4
         KMAX=MAX(KMAX,5)
         DELTA=ABS(XT(J+1)-ALOCAL)
         IF (NSUB.EQ.0)
     1      DELTA=MIN(DELTA,ABS(XT(6)-ALOCAL)+2.0d0*ABS(XT(6)-XT(5)))
         DELTA=Z(1)*DELTA
         TPS=EMEPS*EPSMIN/RNDC
C        TPS DOESN'T DEPEND ON TOLERANCE; INIT. INTERVAL DOESN'T EITHER.
         IF (NSUB.NE.0) TPS=MIN(EPSMIN,TPS)
         ALPHA=(DELTA*ERRCF(2)*(ABS(PHI(7))+ABS(PHI(8))))/TPS
         IF (ALPHA.LT.1.0d0) KK=5
         GO TO 2480
930      CONTINUE
      END IF
      IF (NFEVAL.GT.NFJUMP+1) GO TO 310
      IF (JUMPS.EQ.1) GO TO 1190
C
C     JUMPS HAVE NOT ESCAPED.
C
940   IF (J1.GE.J2) GO TO 970
945   CONTINUE
      IF (J2-J1.LE.1) THEN
         I=J1
         J1=J2
         J2=I
      ELSE
C        JUMPS DO NOT OVERLAP, CHOOSE ONE ON WHICH TO CONVERGE.
         IF (PART.NE.1.AND.STEP(1)*STEP(2).LE.0.0d0) THEN
            J1=18
         ELSE
            J2=0
         END IF
      END IF
C     JUMPS OVERLAP NOW.
970   J1OLD=J1
      IF (J1OLD.EQ.LENDT) J1OLD=LENDT-1
980   J2OLD=J2
      IF (J2OLD.EQ.1) J2OLD=2
C     WE ONLY GET HERE IF WE HAVE AT LEAST 1 JUMP.
      IF (J1.EQ.18) J1=J2+1
      IF (J2.EQ.0) J2=J1-1
      IF (MIN(ISTOP(2,1),LENDT+1-ISTOP(2,2)).LT.(LENDT+2)/3) GO TO 1350
      IF (ABS(SEARCH).NE.4) GO TO 990
c     IF (JUMPS.NE.4) GO TO 1350
c     The above test prevented getting out of the search early on
c     rapidly decaying integrands.
      IF (ISTOP(2,1).LT.ISTOP(2,2)) GO TO 1350
C
C     ESTIMATE ERROR ACROSS JUMP.
C
990   WHERE=0
1000  CONTINUE
      IF (ABS(SEARCH).NE.4.AND.ABS(SEARCH).NE.3) THEN
         IF (JUMPS.GE.3) THEN
c           IF (WHERE.EQ.0) GO TO 1350
c           The above test prevented estimating error if JUMPS == 3
            INC=LENDT
            INC2=1
         END IF
      ELSE
         if (jumps .eq. 3) then
            inc = min(lendt,j1old+1)
            inc2 = max(1,j2old-1)
         else
            INC=MIN(LENDT,J1+1)
            INC2=MAX(1,J2-1)
         end if
      END IF
      TPS=0.0d0
      F1=-1.0d0
      DO 1010 I=INC2,INC-1
         F2=MAX(ABS(XT(I+1)-XT(I)),ABS(EEPSM8*XT(I)))
         F2=ABS((FT(I+1)-FT(I))*F2)
C        F1 IS THE ERROR COMMITTED ON THE TRAPEZOID WITH THE LARGEST ERR
         F1=MAX(F1,F2)
         TPS=TPS+F2
1010  CONTINUE
      TPS=0.5d0*TPS
      IF (MAX(XJ,XJP).GT.0.05d0) TPS=4.0d0*TPS
      X=ABS(0.5d0*(XT(INC)-XT(INC2)))
      IF (INC2.GE.INC-2) THEN
         IF (INC2.GT.1) TPS=TPS+ABS(X*(FT(INC2)-FT(INC2-1)))
         IF (INC.LT.LENDT) TPS=TPS+ABS(X*(FT(INC+1)-FT(INC)))
      END IF
C     COMPUTE THE ALLOWABLE ERROR.
C     DO BLOCK
         IF (WHERE.EQ.0) THEN
C           THE ABSCISSAE DID NOT COALESCE DURING INSERTION OF A POINT
C           IN THE DIFFERENCE LINES
            S=MAX(0.25d0*EPS,EPSMAX)
C           TRY TO GET AT LEAST 3 DIGITS IF THE JUMPS LOOK REAL:
            if (jumps .ne. 3) S=MIN(S, 1.0d-3*(EPSMIN/RNDC))
            IF (SEARCH.NE.6) THEN
               IF (ABS(SEARCH)-3) 1060,1080,1100
1060           CONTINUE
            END IF
C           ABS(SEARCH)=2 OR 6 HERE
            TPS=10.0d0*TPS
            IF (IPRINT.GT.3) CALL DINTO (6,WORK)
            IF (TPS.LE.0.01d0*S) IF (ERR-S) 2240, 2240, 1150
            if (jumps .ge. 3) GO TO 1350
            IF (NSUB.NE.0) IF (ABS(XT(4)-ALOCAL)-DELMIN) 1150,1150,1440
            IF (ABS(XT(9)-ALOCAL).GE.8.0d0*DELMIN) GO TO 1440
C           SET STEPSIZE TO DELMIN.
            DELTA=0.0d0
            GO TO 2620
C           ABS(SEARCH)=3 HERE
1080        CONTINUE
            IF (TPS.LE.EPS) THEN
               IF (XJ.GT.0.0d0) THEN
                  IF (XJP.GT.0.0d0) THEN
                     IF (RE.LE.1.0d0) THEN
                        IF (NSUB.NE.0) THEN
C                          AUGMENT EPSMIN DUE TO NOISE INTRODUCED BY
C                          T**2N TRANSFORMATION.
                           X1=ALOCAL-TA
                           IF (X1.EQ.0.0d0) X1=XT(2)-TA
                           X1=ABS(TB/X1)
                           IF (NSUB.NE.2) X1=X1**(NSUB-1)
                           EPSMIN=EPSMIN+ABS(RNDC*FT(1)*TA*X1)
                        END IF
                        GO TO 2240
                     END IF
                  END IF
               END IF
               IF (ISTOP(2,2).NE.0.AND.TPS.LE.S) GO TO 1150
            END IF
1100        CONTINUE
            IF (SEARCH.LE.0) THEN
               IF (SEARCH.EQ.-3) GO TO 1120
C              THE JUMP SEARCH IS JUST STARTING.  IF THE SUM OF THE
C              ERROR ESTIMATED BY THE QUADRATURE STEP, AND THE ERROR
C              ACROSS THE JUMP IS LESS THAN THE TOLERANCE FOR THE PANEL,
C              ACCEPT THE ANSWER.
               ERR=TPS+ERRC
               IF (ERR.GT.EPS) GO TO 1120
               DISCF=1
               S=EPS
               GO TO 1180
            END IF
         END IF
         IF (ESOLD.GT.0.0d0) COUNT=COUNT+MIN((F1-ESOLD)/ESOLD,1.0d0)
C     END BLOCK
1120  CONTINUE
      ESOLD=F1
      IF (IPRINT.GT.3) CALL DINTO (6,WORK)
      SEARCH=ABS(SEARCH)
      IF (WHERE.EQ.0) THEN
C        THE ABSCISSAE DID NOT COALESCE DURING INSERTION OF A POINT IN
C        THE DIFFERENCE LINES.
         IF (SEARCH.EQ.3) GO TO 1350
         IF (TPS-S) 2670,2670,1350
      END IF
      IF (SEARCH.EQ.3) IF (COUNT) 2237,2237,2190
      IF (COUNT.LE.0.0d0) GO TO 1150
      LOCAL(3)=XT(INC)
      GO TO 2190
C
C     THE ESTIMATED ERROR ACROSS THE JUMP IS SMALL ENOUGH TO ACCEPT.
C     ESTIMATE THE INTEGRAL USING THE COMPOSITE TRAPEZOID RULE.
C     THIS IS SATISFACTORY SINCE THE PRESENCE OF THE JUMP NEGATES
C     ANY POSSIBILITY OF HIGH-ORDER CONVERGENCE.
C
1150  EPS=S
1160  ERR=TPS
      ERRI=ERR
c--D Next line special: S => D, X => Q
      ACUM=0.0D0
      X1=0.0d0
      X2=0.0d0
      J=INC-1
C     DO BLOCK
1170     S=ABS(XT(J+1)-XT(J))
C        COMPUTE EPSMIN ALSO
         F1=ABS(FT(J))
         F2=ABS(FT(J+1))
         X1=X1+ABS((XT(J+1)+XT(J))*(F2-F1))
         X2=X2+S*(F2+F1)
         ACUM=ACUM+(FT(J+1)+FT(J))*S
         J=J-1
         IF (J.GE.INC2) GO TO 1170
C     END BLOCK
      EPSMIN=0.5d0*(X1*EMEPS+X2*RNDC)
c--D Next line special: S => D, X => Q
      ACUM=0.5D0*ACUM
      IEND=.FALSE.
      HAVDIF=.FALSE.
C     SET K SMALL FOR STEPSIZE SELECTION PURPOSES.
C     K=2
1180  FAIL=.FALSE.
      I=INC
      J=INC2
      IF (SEARCH-4) 1240,1230,1240
C
C     IF BOTH JUMPS VANISH, AND THE ORIGINAL ERROR WAS SUFFICIENTLY
C     SMALL, AND ENOUGH OF THE INTERVAL IS BEING EXAMINED, ACCEPT
C     THE ANSWER.
C
1190  IF (ERRC.LE.EPS.AND.10.0d0*ABS(XT(LENDT)-XT(1)).GT.DELTA)
     1   GO TO 2237
C
C     COMPUTE A NEW PANEL TO INTEGRATE.
C
      DISCHK=2
      IF (PART.NE.1) THEN
         IF (STEP(1)*STEP(2).LE.0.0d0) THEN
            IF (J1.EQ.18) J1=0
            J=MAX(2,MAX(J2,J1-2))
            I=LENDT
            IF (JUMPS.NE.1) I=MIN(LENDT,J+3)
            GO TO 1240
         END IF
      END IF
      IF (J2.EQ.0) J2=18
      I=MIN(LENDT-1,MIN(J1,J2+2))
      J=1
      IF (JUMPS.NE.1) J=MAX(1,I-3)
      GO TO 1250
C
C     GET THE BOUNDARIES OF THE SUBDIVISION.
C
1230  IF (PART.EQ.1.OR.STEP(1)*STEP(2).GT.0.0d0) GO TO 1250
1240  L=I
      I=J
      J=L
      X=XJ
      XJ=XJP
      XJP=X
1250  X1=XT(I)
      X2=XT(J)
      FATA=FT(I)
      FATB=FT(J)
      FATAS=.TRUE.
      FATBS=.TRUE.
      WHERE=5
      IF (.NOT.FAIL) THEN
         IF (NSUB.EQ.0) THEN
            DISCX(1)=X1
            DISCX(2)=X2
         ELSE
            DISCX(1)=TA+(X1-TA)*((X1-TA)/TB)**(NSUB-1)
            DISCX(2)=TA+(X2-TA)*((X2-TA)/TB)**(NSUB-1)
         END IF
         IF (SEARCH.EQ.-4) GO TO 2240
         DELTA=ABS(X1-X2)
         IF (IPRINT.GT.1) CALL DINTO (7,WORK)
         IF (SEARCH.NE.4) THEN
            BLOCAL=X2
            ABSDIF=0.5d0*DELTA
            IF (MIN(ABS(XJ),ABS(XJP)).LE.1.0D-3) DISCF=1
            IF (SEARCH.EQ.2.AND.NSUB.NE.0) CALL DINTNS (NSRB)
            GO TO 2240
         END IF
      END IF
C     DECREASE NSUB TO ZERO.
1260  IF (NSUB.NE.0) THEN
         CALL DINTNS (NSRB)
         go to 1260
      END IF
      IF (FAIL) THEN
C
C        TRY NOT TO CHOOSE A PANEL TO BE INTEGRATED THAT WILL LEAVE
C        A RIDICULOUSLY SMALL UN-INTEGRATED PIECE.
C
C        DO BLOCK
            IF (ABS(X2-ALOCAL).GE.1.5d0*DELMIN) THEN
               IF (ABS(X2-X1).LE.8.0d0*ABS(X1-ALOCAL)) THEN
                  SMIN=DINTSM(X1)
                  IF (ABS(X1-BLOCAL).GE.1.5d0*SMIN) THEN
                     IF (ABS(X2-X1).LE.8.0d0*ABS(X1-BLOCAL)) GO TO 1280
                  END IF
                  IF (PART.EQ.2) GO TO 1270
                  IF (.NOT.IEND) GO TO 1270
                  X1=X2
               END IF
            END IF
            DELTA=ABS(X1-ALOCAL)
            FATAS=.FALSE.
            GO TO 2610
C        END BLOCK
1270     CONTINUE
         X1=BLOCAL
         FATAS=.FALSE.
      END IF
1280  CONTINUE
      FATBS=FATAS
C
C     SPLIT THE PROBLEM INTO TWO PARTS AND INTEGRATE AWAY FROM
C     THE MIDDLE IN BOTH DIRECTIONS.
C
C     DO BLOCK
         IF (PART.NE.1) THEN
            I=2
c##         EPSR=EPSR+ERRT(2)
            IF (STEP(1)*STEP(2).LE.0.0d0) GO TO 1290
         END IF
         END(2)=ALOCAL
         I=1
C     END BLOCK
1290  CONTINUE
      TLEN=ABS(END(2)-END(1))
      FSAVE=FATA
      FSAVED=FATAS
      S=START(1)
      START(1)=X1
      IF (FAIL) THEN
C        Jump disappeared.  X1..X2 not yet integrated.
         START(2)=X1
! Why was the following ever done ???
!         TA=X1+0.5d0*(X2-X1)
         TA=X1
         TASAVE=TA
      ELSE
C        Jump didn't disappear.  X1..X2 was integrated using trapezoid
         START(2)=X2
         TA=START(1)
         TASAVE=START(2)
      END IF
      TB=END(2)-TA
      F1=WORRY(1)
      WORRY(1)=LOCAL(3-I)
      IF (ABS(F1-S).GE.ABS(X1-S)) THEN
         IF (ABS(F1-S).LT.ABS(WORRY(1)-S)) WORRY(1)=F1
      END IF
      WORRY(2)=LOCAL(I)
      DELTA=ABS(X1-X2)
      STEP(2)=SIGN(1.0d0,SIGN(1.0d0,TB)*(BINIT-AINIT))
      ALOCAL=X1
      BLOCAL=X2
      PART=2
      DID1=.FALSE.
      INIT=.FALSE.
      NSUBSV=0
      IF (.NOT.FAIL) THEN
         IF (XJ.GT.0.0d0) NSUBSV=2
         IF (XJP.GT.0.0d0) CALL DINTNS (NSIB)
      END IF
      ABSDIF=0.5d0*DELTA
      ERRT(2)=0.0d0
      RESULT(2)=0.0d0
      XCDOBT(2)=.FALSE.
      DISCHK=MAX(DISCHK,1)
      IF (.NOT.FAIL) THEN
         SMIN=DINTSM(BLOCAL)
         IEND=ABS(END(2)-BLOCAL).LE.SMIN
C        SOMEDAY, WE SHOULD USE A TRAPEZOIDAL ESTIMATE ACROSS ANY
C        TINY PIECE LEFT UN-INTEGRATED DUE TO THE ABOVE STATEMENT.
         IF (JUMPS .EQ. 3) then
            DISCF = -2
         else
            DISCF=2
         end if
         S=WORRY(1)
         IF (XJ.LE.0.0d0) S=END(1)
C        SET THE RESTARTING STEPSIZE FOR PART 1.
         SMIN=DINTSM(START(1))
         STEP(1)=SIGN(MAX(ABS(S-START(1)),SMIN),STEP(1))
         GO TO 2250
      END IF
C     SET THE RESTARTING STEPSIZE FOR PART 1 TO THE STEPSIZE FOR PART 2.
      SMIN=DINTSM(X1)
      STEP(1)=SIGN(MAX(DELTA,SMIN),STEP(1))
      KMIN=0
      GO TO 2615
C
C     THE ESTIMATED ERROR ACROSS THE INTERVAL IS TOO LARGE.  REDUCE
C     THE INTERVAL LENGTH AND ADD A FUNCTION VALUE IN THE NEIGHBORHOOD
C     OF THE JUMP.  THEN RE-FORM AND RE-EXAMINE THE DIFFERENCES.
C
1350  I=J1
      J=J2
C     DO FOREVER
1360  CONTINUE
C        DO BLOCK
            IF (I.NE.LENDT) THEN
               IF (J.NE.1) THEN
                  IF (JUMPS.GE.3) THEN
                     X1=ABS(XT(I)-XT(I+1))
                     X2=ABS(XT(J)-XT(J-1))
                     IF (0.3d0*X2.GE.MIN(X1,ABS(XT(J+1)-XT(J)))) THEN
                        I=J-1
                        GO TO 1400
                     END IF
                     IF (0.3d0*X1.GE.MIN(X2,ABS(XT(I)-XT(I-1)))) THEN
                        J=I
                        I=I+1
                        GO TO 1400
                     END IF
                  END IF
               END IF
            END IF
            IF (I-J-1) 1390,1450,1410
1390        I=J+1
            IF (ABS(XT(J)-XT(J+1)).LT.ABS(XT(J)-XT(J-1))) I=J-1
C        END BLOCK
1400     CONTINUE
         ABSCIS=XT(J)+0.25d0*(XT(I)-XT(J))
         GO TO 1490
1410     CONTINUE
         IF (ABS(XJ).LE.ABS(XJP)) THEN
            XJ=5.0d0*XJ
            I=I-1
         ELSE
            XJP=5.0d0*XJP
            J=J+1
         END IF
      GO TO 1360
C     END FOREVER
1440  J=1
1450  I=J+1
      ABSCIS=XT(I)+0.5d0*(XT(J)-XT(I))
1490  L=MAX(I,J)
      IF (SEARCH.LE.0) THEN
         SEARCH=-SEARCH
C        IF (SEARCH.EQ.4) THEN
C          DETERMINE WHETHER NODES HAVE COALESCED.  WE ASSUME THAT NODES
C          WILL NOT COALESCE WHEN SEARCH = 2.
C        END IF
      END IF
1500  IF (ABSCIS.EQ.XT(I).OR.ABSCIS.EQ.XT(J)) THEN
         X1=XT(I)+0.5d0*(XT(J)-XT(I))
         IF (X1.EQ.ABSCIS) GO TO 2680
         ABSCIS=X1
         go to 1500
      END IF
1520  CONTINUE
      IF (LENDT.GT.7) THEN
C        DECIDE HOW MUCH OF THE DIFFERENCE LINE TO THROW AWAY.
C        THROW INDICATES WHICH END TO TRIM - 1=ALOCAL, 2=BLOCAL, 3=BOTH.
         THROW=0
         IF (JUMPS.GE.3) IF (ABS(SEARCH)-4) 1555,1530,1555
         IF (SEARCH.EQ.2.OR.SEARCH.EQ.6) IF (LENDT-9) 1720,1680,1630
1530     CONTINUE
         IF (ISTOP(2,2).GT.ISTOP(2,1)+1) THEN
C           IF JUMPS DO NOT OVERLAP, THROW AWAY A LOT.
            IF (PART.NE.1) THEN
               IF (STEP(1)*STEP(2).LE.0.0d0) THEN
                  J1=ISTOP(2,2)-4
                  J2=0
                  IF (J1) 1610,1610,1640
               END IF
            END IF
            J2=LENDT-ISTOP(2,1)-3
            J1=0
            IF (J2) 1610,1610,1640
         END IF
         IF (JUMPS.EQ.1) GO TO 1610
1555     J=J2
         J2=MAX(MIN(LENDT-L,(IC2-J1-1)/2),0)
         J1=MAX(MIN(L-1,(J-IC1-1)/2),0)
         IF (SEARCH.NE.4) THEN
            IF (LENDT.LT.17) J1=0
         END IF
C        DO BLOCK
            IF (J2-1) 1580,1570,1640
1570        THROW=2
1580        IF (J1-1) 1600,1590,1640
1590        THROW=THROW+1
1600        IF (THROW.NE.0) IF (THROW-2) 1670,1680,1670
C        END BLOCK
1610     CONTINUE
         IF (LENDT.LT.15) GO TO 1740
         IF (JUMPS.NE.1) THEN
            IF (ISTOP(2,1).EQ.18) GO TO 1680
            IF (ISTOP(2,2).EQ.0) GO TO 1670
         END IF
         IF (LENDT+1-J1OLD-J2OLD) 1670,1660,1680
1630     J1=0
         J2=LENDT-9
C        THROW OUT J1 NODES FROM THE HEAD AND J2 NODES FROM THE TAIL
C        OF THE DIFFERENCE LINES.
1640     LENDT=LENDT-J1-J2
C        FORCE RE-FORMATION OF DIFFERENCES FROM XT AND FT.
         NFJUMP=NFEVAL-10
         IF (J1.EQ.0) GO TO 1740
         IF (SEARCH.NE.4) GO TO 520
         IF (J1OLD.NE.18) J1OLD=J1OLD-J1
         J2OLD=J2OLD-J1
         L=L-J1
         DO 1050 I=1,LENDT
            XT(I)=XT(I+J1)
            FT(I)=FT(I+J1)
1050     CONTINUE
         GO TO 1740
1660     THROW=3
1670     IF (SEARCH.NE.4) GO TO 520
C        THROW OUT X(1) AND F(1), RE-FORM THE DIFFERENCE LINES.
         KK=1
         INC=1
         J=1
         J1=1
         J2=LENDT
         FATAS=.FALSE.
         L=L-1
         IF (J1OLD.NE.0) J1OLD=J1OLD-1
         IF (J2OLD.NE.0) J2OLD=J2OLD-1
         GO TO 1690
C        THROW OUT X(LENDT) AND F(LENDT), RE-FORM THE DIFFERENCE LINES.
1680     KK=2
         INC=-1
         J=17+LENDT
         J1=LENDT
         J2=1
         FATBS=.FALSE.
1690     BETA=1.0d0
         I=J1
         PHI(J)=PHI(J)-PHI(J+INC)
         J=J+INC
         I=I+INC
C        DO BLOCK
1700        BETA=BETA*((XT(J1+INC)-XT(I+INC))/(XT(J1)-XT(I)))
            PHI(J)=BETA*(PHI(J)-PHI(J+INC))
            J=J+INC
            I=I+INC
            IF (I.NE.J2) GO TO 1700
C        END BLOCK
         LENDT=LENDT-1
         DO 1710 I=1,LENDT
            PHIT(I)=PHIT(I+1)
            IF (INC.GE.0) THEN
               XT(I)=XT(I+1)
               FT(I)=FT(I+1)
            END IF
1710     CONTINUE
         THROW=THROW-KK
         IF (THROW.GT.0) GO TO 1680
C        THE NEXT STATEMENT OCCASSIONALLY CAUSES A DIFFERENCE LINE TO
C        BE SHORTENED, BUT NO NEW POINT ADDED.
         IF (KK+L.EQ.2) GO TO 360
1720     CONTINUE
      END IF
1740  CONTINUE
      IF (L.GT.LENDT) GO TO 310
C     ADD A POINT IN THE INTERIOR OF THE INTERVAL.
      X=ABSCIS
      WHERE=6
      GO TO 2850
C
C     NO JUMPS, OR ONLY ONE INTERIOR JUMP, IN THE ORIGINAL
C     DIFFERENCE LINE FOR THIS PANEL.
C
1920  CONTINUE
C     DO BLOCK
         IF (RE.LE.1.0d0) THEN
            IF (K.LE.5.AND.KAIMT.GT.5.AND.KONVRG.LT.4)
     1         ERR=MAX(ERR,ERRI*(ERRI/EPSMIN))
            IF (ERR.GT.EPS) THEN
               IF (K.LT.4) GO TO 100
               IF (K.GT.5) GO TO 2060
               IF (KONVRG-6) 2060,2060,2130
            END IF
C
C           ERROR CRITERIA HAVE APPARENTLY BEEN SATISFIED.  LOOK AT THE
C           DIFFERENCE LINES TO AVOID BEING FOOLED BY ANSWERS FROM
C           SUCCESSIVE FORMULAE BEING NEARLY EQUAL BUT BOTH WRONG.
C
            IF(KONVRG.LE.LENDT.AND.K.LT.KMIN.AND.PHISUM.LE.5.0d0*PHTSUM)
     1         GO TO 100
            IF (.NOT.INIT.AND.NSUB.EQ.0) THEN
               IF (ERR.GT.EPSMIN) THEN
                  I=KAIMT-2
                  IF (KONVRG.GT.LENDT+LENDT-5) I=KAIMT-3
                  IF (K-I) 100,1950,1960
1950              IF (ENDPTS-2) 2750,2760,1960
1960              CONTINUE
               END IF
            END IF
            IF (K.LE.4) THEN
               EXTRA=ERRCF(K-1)
               IF (K.NE.2) THEN
                  IF (KONVRG.LT.K+K-1+ENDPTS/2) EXTRA=10.0d0*EXTRA
               END IF
               EXTRA=EXTRA*DELTA*(ABS(PHI(LENDT-1))+ABS(PHIT(2)))
               J=(K-2)*(KONVRG-7)
               IF (J.GT.20) THEN
                  EXTRA=0.0d0
               ELSE IF (J.GT.0) THEN
                  EXTRA=EXTRA*(0.1d0**J)
               END IF
               IF (JUMPS.GE.3) GO TO 2140
               IF (EXTRA.GE.32.0d0*EPSMIN) THEN
                  S=ERR+EXTRA
                  IF (S.GT.EPS) IF (K-KAIMT) 2730,2050,2050
                  ERR=S
               END IF
            END IF
            IF (JUMPS.EQ.3) GO TO 750
C           DO BLOCK
               IF (DISCHK.LE.0) THEN
                  IF (DISCHK.GE.-1) THEN
                     IF (.NOT.INIT) GO TO 2235
                     IF (K.GE.4) GO TO 2235
                     KMAX=MAXK
                     GO TO 2030
                  END IF
               END IF
C              COMPUTE F AT (OR NEAR) THE ENDS OF THE PANEL, UNLESS
C              THERE IS A JUMP AT THE OTHER END OF THE PANEL.
               KMAX=KAIMT+1
C           END BLOCK
2030        CONTINUE
            IF (ENDPTS-2) 2750,2760,2235
         END IF
C
C        CHECK FOR NOISE IN THE INTEGRAND.
C
2050     IF (K.LT.4) GO TO 100
2060     CONTINUE
         IF (RE.GE.0.01d0) THEN
C           DO BLOCK
               IF (RE.LT.1.0d0) THEN
                  IF (NSUB.NE.0.AND.
     1            ABS((LOCAL(3)-TA)/(LOCAL(4)-LOCAL(3)))
     2            .LE.1.0d0) GO TO 2130
               ELSE IF (ERRI.LT.100.0d0*EPSMIN) THEN
                  IF (4.0d0*ERR.GT.EPSMIN) GO TO 2100
                  IF (.NOT.INIT.AND.DISCHK.LE.0) GO TO 2250
                  IF (ENDPTS-2) 2750,2760,2250
               END IF
               IF (REPROD.GE.1.0D-4) GO TO 2130
               IF (K.LT.KAIMT.AND.ERR.GT.100.0d0*EPSMIN) GO TO 2130
               IF (ERRI*REPROD.GT.0.01d0*EPNOIZ) GO TO 2130
               IF (RE*ERR.LT.EPS) GO TO 2130
C           END BLOCK
2100        CONTINUE
            IF (ENDPTS-2) 2750,2760,2120
C           NOISE IN THE INTEGRAND - RESET THE ROUND-OFF CONSTANT.
2120        RNDC=4.0d0*RNDC*(ERR/EPSMIN)
            IF (IPRINT.GT.0.and.100.0d0*err.gt.eps) CALL DINTO (8,WORK)
            EPS=ERR
            ROUNDF=.TRUE.
            DELTA=DELTA+DELTA
            GO TO 2250
         END IF
C     END BLOCK
C     NO NOISE DETECTED.
2130  IF (DELTA.LT.DELMIN+DELMIN) GO TO 2250
      KMAX=MAXK
      IF (DISCHK.GT.0.OR.DISCHK.LT.-1) KMAX=KAIMT+1
C     IF WE HAVE A LOT OF CONVERGENCE, WE SHOULD MAKE IT.
C     IF WE DONT MAKE IT, SOMETHING FISHY IS HAPPENING.
C     ENDPTS IS INCLUDED BECAUSE CONVERGENCE IS EASIER WHEN THE END
C     POINTS HAVE BEEN ADDED TO THE DIFFERENCE LINES.
      IF (KONVRG.GE.ENDPTS+7.AND.RE.GT.1.0D-3) KMAX=5
      IF (10.0d0*RE*RE*ERR.LT.REP*EPS) KMAX=MIN(MAXK,MAX(KMAX,K+1))
      IF (JUMPS.LT.3) THEN
         IF (K.LT.KMAX) THEN
            IF (DISCHK.GE.-1) THEN
      IF (ABS(WORRY(PART)-START(PART))-DELMIN.GT.ABS(BLOCAL-START(PART))
     1.OR.DID1) GO TO 2730
               IF (K.LT.KAIMT) THEN
                  S=RE*ERR
                  IF (K.LT.KAIMT-1) S=S*RE/REP
                  IF (S.LE.EPS) GO TO 2730
               END IF
            END IF
         END IF
      END IF
C
C     USE THE SEARCH MECHANISM TO FIND THE PROPER INTERVAL LENGTH.
C
2140  IF (ENDPTS.LT.2) GO TO 2750
      ENDPTS=3
      SEARCH=2
2150  GO TO (2160,620,970,2165), JUMPS
2155  SEARCH=2
2160  INC=4
      INC2=1
      IF (JUMPS.NE.1) GO TO 990
      IF (IC1+IC1+NSUB-12) 990,2690,2690
2165  IF (J1.LT.J2) GO TO 945
      XJUMP=XJUMPS
      GO TO 870
C
C     ABSCISSAE HAVE COALESCED.
C
2170  IF (K.LT.4) GO TO 100
2180  DELTA=100.0d0*DELTA
      IF (RE.LT.0.875d0) GO TO 2200
      IF (ERRI.LT.MAX(EPSMIN,PEPSMN)) GO TO 2200
      IF (ERR.LT.MAX(EPSMIN,PEPSMN)) GO TO 2200
      IF (REPROD.LE.1.0D-4) GO TO 2200
      IF (HAVDIF.AND.ABS(JUMPS-3).NE.1) GO TO 2200
2190  CONTINUE
      CALL DINTO (9,WORK)
      IFLAG(1)=5+KDIM
      WORK(1)=EMINF
      ANSWER=ALOCAL
      GO TO 2340
2200  CONTINUE
      IF (WHERE.LT.0) THEN
C        WHERE .LT. 0 MEANS ABSCISSAE COALESCED DURING THE QUADRATURE
C        STEP.  ACUM IS GARBAGE DURING THE QUADRATURE STEP SO USE PACUM.
         ACUM=PACUM
         AACUM = abs(ACUM)
         EPSMIN=PEPSMN
      ELSE
         TPS=ACUM-PACUM
         ACUM=PACUM
         TPS=(RE*TPS)/(1.0d0-RE)
         IF (ABS(RE-REP).GE.0.125d0*RE) THEN
C           DO NOT TRUST THE ERROR ESTIMATE.
            ERR=ERR+ABS(TPS)
         ELSE
C           TRUST THE ERROR ESTIMATE, APPLY AITKEN ACCELERATION.
            ACUM=ACUM-TPS
         END IF
      END IF
      IF (IPRINT.GT.1) CALL DINTO (10,WORK)
      GO TO 2250
C
C     INTEGRATION WAS SUCCESSFUL ON LAST INTERVAL.  CHECK IF THE
C     COMPLETE INTERVAL HAS BEEN INTEGRATED.
C
2235  IF ((ERR*RNDC).GT.EPSMIN) DELTA=1000.0d0*DELTA
      GO TO 2240
2237  ERR=ERRC
2240  DISCHK=MAX(DISCHK-1,0)
2250  SEARCH=1
      I=PART
      IF (PART.NE.1) THEN
         IF (STEP(1)*STEP(2).GT.0.0d0) I=1
      END IF
      ERRT(I)=ERRT(I)+MAX(ERR,EPSMIN)
c     At one time, we had max(err,epsmin,edmeps*result(i)), where edmeps
c     is 0.5*d1mach(4), to account for round-off introduced during
c     accumulation.  This effect should be negligible.
      XCDOBT(I)=XCDOBT(I).OR.ERR.GT.(EPSMIN/AACUM)**RELOBT*AACUM
      TPS=MAX(EPSO,ABS(RELEPS))
      epsr=errt(1)+errt(2)+edue2a+edue2b
      EPSR=MAX(TPS-epsr,MAX(0.1d0*epsr,0.0d0))
      FAIL=.FALSE.
      KMIN=0
      RESULT(I)=RESULT(I)+ACUM
      IF (IPRINT.GT.1) CALL DINTO (11,WORK)
      IF (DISCF.gt.0) THEN
         IF (PART.EQ.1) THEN
            DISCF=0
         ELSE
            IF (.NOT.IEND) GO TO 2350
         END IF
         IF (IPRINT.GT.0) CALL DINTO (13,WORK)
      END IF
      IF (.NOT.IEND) GO TO 2350
      IF (PART.NE.1) GO TO 2612
2320  ANSWER=RESULT(1)
      IF (BINIT.LT.AINIT) ANSWER=-ANSWER
      WORK(1)=ERRT(1)+EDUE2A+EDUE2B
C     DO BLOCK
         IF (WORK(1).GT.EPSO) THEN
            IFLAG(1)=-2
            IF (XCDOBT(1)) IFLAG(1)=-3
            IF (RELTOL.EQ.0) GO TO 2340
            TPS=ABS(RELEPS)
            RELEPS=ABS(ANSWER)-WORK(1)
            RELEPS=MAX(ABS(WORK(RELTOL+2))*RELEPS,ENZER)
            IF (WORK(1).GT.RELEPS) THEN
               IF (WORK(1).GT.TPS) GO TO 2340
               IF (IFLAG(1).GE.-2) GO TO 2340
               GO TO 10
            END IF
         END IF
         IFLAG(1)=-1
C     END BLOCK
2340  WHERE=0
2345  IF (NFINDX.NE.0.AND.KDIM.EQ.1) IFLAG(NFINDX)=NFEVAL
      RETURN
C
C     THE INTEGRATION IS COMPLETE THROUGH THE PANEL ENDING AT BLOCAL.
C
2350  FATA=FATB
      FATAS=FATBS
      TLEN=TLEN-ABSDIF-ABSDIF
      IF (abs(DISCF).GE.2) THEN
         DISCF=sign(1,discf)
         ALOCAL=BLOCAL
         EPS=EPSR*(ABS(END(PART)-ALOCAL)/TLEN)
         IF (PART-1) 40,40,50
      END IF
      KMAX=MAXK
      KK=K
      IF (DID1) THEN
         IF (K.LT.KAIMT) THEN
            IF (ERR.GT.ERSQE6*EPSMIN) KK=KK+1
         END IF
      ELSE IF (K.LT.6) THEN
         KK=MIN(6,K+NSUB/2)
      END IF
      IF (NSUB.NE.0) THEN
C        DO BLOCK
            IF (ABS(BLOCAL-START(PART)).LT.ABS(TEND-START(PART))-DELMIN)
     1         THEN
               IF (.NOT.HAVDIF) GO TO 2390
               IF (PHISUM.LT.DBLE(10*(K-KAIMT))*PHTSUM) GO TO 2380
               IF (PHISUM.GE.16.0d0*PHTSUM) GO TO 2390
               IF (PHISUM.GE.PHTSUM*0.5d0**(KAIMT-K+1)) GO TO 2390
               GO TO 2380
            END IF
            TEND=END(PART)
            KMAX=KMAXF
C        END BLOCK
2380     CONTINUE
         HAVDIF=.FALSE.
         CALL DINTNS (NSRB)
      END IF
2390  CONTINUE
      ALOCAL=BLOCAL
      ERRC=MAX(ERR,ERRI)
      EPNOIZ=MIN(ERRC,4.0d0*EPNOIZ)
      IF (.NOT.ROUNDF) THEN
         IF (REP.GT.0.3d0) KMIN=K-1
         S=10.0d0*ERRI
         IF (S.LE.EPSMIN.AND.RNDC.GT.FER) RNDC=MAX(FER,S*RNDC/EPSMIN)
      END IF
      S=MAX(REP,1.0D-5)
      TP1=1.0d0
C     DO BLOCK
2420     TP1=TP1*1.125d0
         S=S*6.0d0
         IF (S.LT.0.4d0) GO TO 2420
C     END BLOCK
C     DO BLOCK
         IF (K.GE.MAX(KAIMT-1,3).AND.K.LE.AIMFOR.AND.ERRI.LE.EPSMIN.AND.
     1      REP.LE.0.1d0) THEN
            I=AIMFOR-K+4
            KK=K-1
         ELSE
            I=AIMFOR-KK+3
C           1 .LE.  I  .LE. 7
            IF (I.GE.3) THEN
               IF (I.EQ.3.AND.ERRI.GE.EPSMIN.AND.REP.GT.0.3d0)GO TO 2450
               IF (ERR.EQ.0.0d0) GO TO 2450
               IF (ERR.GT.EPSMIN) GO TO 2450
               IF (ERRI.GT.10.0d0*EPSMIN) GO TO 2450
C              ADJUST STEP SIZE WHEN APPROACHING NOISE OR ROUND-OFF
C              LIMIT.
               DELTA=DELTA*TP1
               GO TO 2450
            END IF
         END IF
         ERRC=MAX(ERRC,ABS(EP))
         ERR=ERRC
         RE=REP
         IF (REP.GT.1.0d0) DELTA=0.875d0*DELTA
C     END BLOCK
2450  CONTINUE
      TPS=EPS
      IF (K.GE.KK) TPS=MIN(ESQEPS*EPSMIN/RNDC,EPS)
      ALPHA=0.0d0
      IF (TPS.NE.0.0d0) ALPHA=MAX(ERRC*MIN(1.0d0,SQRT(RE)),EPSMIN)/TPS
2480  DELTA=DELTA*XSTEP(I)
      TP=ABS(ALOCAL-START(PART))
      KAIMT=AIMFOR
C     DO BLOCK
      IF (SEARCH.NE.1.OR.TP.LT.ABS(WORRY(PART)-START(PART))-DELMIN) THEN
C
C           CHOOSE THE NEXT STEPSIZE BASED UPON THE ORDER OF THE CURRENT
C           FORMULA (SOME FUNCTION OF K), AND THE RATIO OF THE ERROR
C           COMMITTED TO THE ERROR REQUESTED ON THE CURRENT PANEL.
C           THE ERROR BEHAVES AS  H**ORDER(K-1) WHERE H IS THE STEPSIZE.
C           THEN SET GAMMA = 1.125**ORDER(K-1), AND CHANGE THE ERROR
C           RATIO BY GAMMA AND THE STEPSIZE BY 1.125 UNTIL THE ERROR
C           RATIO IS GREATER THAN 2.0.
C
            IF (ALPHA.NE.1.0d0) THEN
               IF (ALPHA.LT.1.0d0) THEN
                  IF (ALPHA.EQ.0.0d0) GO TO 2500
                  BETA=1.125d0
                  ALPHA=ALPHA*GAMMA(KK-1)
               ELSE
                  ALPHA=1.0d0/ALPHA
                  BETA=1.0d0/1.125d0
                  DELTA=DELTA*BETA
               END IF
2490           CONTINUE
                  ALPHA=ALPHA*GAMMA(KK-1)
                  IF (ALPHA.GT.2.0d0) GO TO 2500
                  DELTA=DELTA*BETA
               GO TO 2490
            END IF
2500        CONTINUE
            IF (SEARCH.NE.1) THEN
               DISCHK=MAX(DISCHK,0)
               IF (IPRINT.GT.3) CALL DINTO (12,WORK)
               IF (SEARCH.EQ.6) GO TO 2560
               KMAX=MIN(MAXK-1,KMAX)
               IF (DELTA.LE.Z(2)*ABSDIF
     1         .OR.(K.LT.6.AND.DELTA.LT.Z(3)*ABSDIF)) GO TO 2610
               IF (DELTA.GT.Z(4)*ABSDIF) KMAX=MIN(KMAX,5)
               KAIMT=MIN(KMAX,KAIMT)
               IF (K.GE.KMAX) GO TO 2560
               DELTA=ABSDIF+ABSDIF
               ERR=ERRC+EXTRA
               IF (K.GE.4) DISCHK=-1
               IF (ERR.GT.EPS) GO TO 100
               IF (RE.GT.1.0d0) GO TO 100
               IF (KONVRG-6) 100,100,2240
            END IF
        IF (KK.LT.K) DELTA=MAX(DELTA,2.0d0*TP1*XSTEP(AIMFOR-K+3)*ABSDIF)
            IF (DELTA.GT.ABSDIF+ABSDIF) KMIN=K
            IF (TP+DELTA.LT.ABS(END(PART)-START(PART))) THEN
         IF (MIN(ABS(TB),TP+DELTA).GE.ABS(WORRY(PART)-START(PART))) THEN
                  DISCHK=-1
                  KMAX=KMAXF
                  IF (TP+1.25d0*DELTA.GE.ABS(TB)) GO TO 2540
               END IF
            END IF
C           DO BLOCK
               IF (TP.LT.ABS(TEND-START(PART))-DELMIN) THEN
                  IF (K.GE.KAIMT) GO TO 2560
                  IF (.NOT.HAVDIF) GO TO 2560
                  IF (PHISUM.LE.10.0d0*PHTSUM) GO TO 2560
                  IF (NSUB.GE.NSUBMX) GO TO 2560
                  BLOCAL=ALOCAL+SIGN(DELTA,TB)
                  CALL DINTNS (NSIB)
                  GO TO 2620
               END IF
               TEND=ALOCAL+SIGN(DELTA,TB)
               KMAX=KMAXF
C           END BLOCK
            IF (DID1.OR.PART.EQ.1.OR.FAIL) GO TO 2610
C           AFTER THE FIRST SUCCESS IN PART 2, TREAT THE NEXT INTERVAL
C           LIKE A RESTART.  ALSO, INCREASE THE STARTING STEPSIZE FOR
C           RESTARTING PART 1 IF POSSIBLE.
            KAIMT=3
            KMAX=KMAXF
            X1=ALOCAL+SIGN(DELTA,TB)
            X2=ALOCAL
            IF (NSUB.NE.0) THEN
               X1=TDECR(X1)
               X2=TDECR(X2)
               IF (NSUB.NE.2) THEN
                  X1=TDECR(X1)
                  X2=TDECR(X2)
               END IF
            END IF
            STEP(1)=SIGN(MAX(ABS(STEP(1)),ABS(X1-X2)),STEP(1))
            GO TO 2610
         END IF
2540     continue
         DISCHK=-1
         KMAX=KMAXF
         DELTA=MAX(DELTA,ABS(END(PART)-ALOCAL))
         IF (DELTA.GT.3.0d0*ABSDIF) DISCHK=-2
         WORRY(PART)=END(PART)
         GO TO 2610
C     END BLOCK
2560  CONTINUE
      DELTA=MIN(ABSDIF,DELTA)
      WORRY(PART)=BLOCAL
      KMAX=KMAXF
      KAIMT=MIN(MAX(5,KMAXF),AIMFOR)
2610  FATBS=.FALSE.
      GO TO 2620
2612  IEND=.FALSE.
      PART=1
      NSUB=0
      DISCF=0
      DISCHK=MAX(DISCHK,1)
      RESULT(1)=RESULT(1)+RESULT(2)
      ERRT(1)=ERRT(1)+ERRT(2)
      ERRT(2)=0.0d0
      XCDOBT(1)=XCDOBT(1).OR.XCDOBT(2)
      ALOCAL=START(1)
      TA=TASAVE
      TB=END(1)-TA
      TLEN=ABS(END(1)-START(1))
      DELMIN=DINTSM(ALOCAL)
2614  IF (NSUB.LT.NSUBSV) THEN
         CALL DINTNS (NSIA)
         GO TO 2614
      END IF
      IF (DID1) DELTA=0.0d0
      DELTA=MAX(ABS(STEP(1)),MAX(10.0d0*DELMIN,DELTA))
      STEP(1)=SIGN(1.0d0,SIGN(1.0d0,TB)*(BINIT-AINIT))
      FATBS=.FALSE.
      FATAS=FSAVED
      FATA=FSAVE
      FSAVED=.FALSE.
      DID1=.FALSE.
2615  KAIMT=3
      KMAX=KMAXF
2620  X2=ABS(END(PART)-ALOCAL)
C     X2 IS THE LENGTH OF THE REST OF THE PART.
      INIT=.FALSE.
      EPS=EPSR*X2/TLEN
      DELMIN=DINTSM(ALOCAL)
      IF (DELTA.LT.X2) THEN
         IEND=.FALSE.
         IF (DELTA.LT.DELMIN) THEN
            DELTA=DELMIN
            FATBS=.FALSE.
         END IF
C        DO NOT LEAVE A PIECE TOO SMALL.
         IF (X2-DELTA.LE.DELMIN) GO TO 60
         BLOCAL=ALOCAL+SIGN(DELTA,TB)
         GO TO 75
      END IF
C     REDUCE KAIMT IF DELTA IS VERY MUCH GREATER THAN THE REMAINING
C     PART OF THE MAJOR SUBDIVISION.
C     DO FOREVER
2640  CONTINUE
         X2=3.0d0*X2
         IF (X2.GT.DELTA) GO TO 60
         KAIMT=KAIMT-1
         IF (KAIMT.LT.3) GO TO 60
      GO TO 2640
C     END FOREVER
C
C     *****     INTERNAL SUBROUTINES     *******************************
C
C     COMPUTE MINIMUM RATIO OF STEPSIZE CHANGES.  IF THE STEPSIZE
C     CHANGES MORE RAPIDLY THAN 0.3, ADD A POINT IN THE LARGEST
C     INTERVAL.
C
2670  WHERE=3
      KK=MAX(INC2+1,3)
      J=MAX(INC,3)
      EPS=S
      GO TO 2710
2680  WHERE=2
C     ABSCISSAE HAVE COALESCED.
      J=LENDT
      GO TO 2700
2690  WHERE=1
      J=LENDT
      IF (JUMPS.NE.1) GO TO 2700
      IF (SEARCH.NE.2.AND.SEARCH.NE.6)
     1 IF (KONVRG+KONVRG-LENDT) 2700,2700,1190
      KK=4
      J=6
      GO TO 2710
2700  KK=3
2710  TP=0.3d0
      TP1=ABS(XT(KK-1)-XT(KK-2))
      DO 2720 I=KK,J
         S=ABS(XT(I)-XT(I-1))
         IF (TP*S.GE.TP1) THEN
            TP=TP1/S
            L=I
         ELSE IF (TP*TP1.GT.S) THEN
            TP=S/TP1
            L=I-1
         END IF
         TP1=S
2720  CONTINUE
      IF (TP.EQ.0.3d0) GO TO (910,1000,1160), WHERE
      ABSCIS=XT(L-1)+0.5d0*(XT(L)-XT(L-1))
      IF (WHERE.NE.2) THEN
         IF (ABSCIS.EQ.XT(L-1)) GO TO 1000
         IF (ABSCIS.EQ.XT(L)) GO TO 1000
      END IF
      IF (JUMPS.EQ.1) GO TO 1520
      IF (J1.EQ.18) J1=J2+1
      IF (J2.EQ.0) J2=J1-1
      GO TO 1520
C
C     REQUEST A FUNCTION VALUE AT THE EDGE OF A PANEL.
C
2730  IF (DISCHK.LE.0) GO TO 100
2740  IF (ENDPTS-2) 2750,2760,100
2750  IF (J1.EQ.LENDT) GO TO 2760
      WHERE=5
      WHERE2=1
      ENDPTS=2
      IF (J2.EQ.1) ENDPTS=3
      GO TO 2810
2760  WHERE=5
      WHERE2=2
      ENDPTS=3
      IF (J2.EQ.1) GO TO 530
2810  ABSCIS=LOCAL(WHERE2)
      X=ABSCIS
      IF (FATS(WHERE2)) GO TO 2830
      GO TO 2850
2820  FAT(WHERE2)=FNCVAL
      FATS(WHERE2)=.TRUE.
      ERRAT(WHERE2)=FER*ABS(FNCVAL)
      IF (FEA.NE.0) ERRAT(WHERE2)=ERRAT(WHERE2)+ABS(T*WORK(FEA))
      IF (LOCAL(WHERE2).EQ.AINIT) THEN
         EDUE2A=FNCVAL*ERRINA
         EPS=MAX(EPSR-(EDUE2A+EDUE2B),EDUE2A+EDUE2B)*DELTA/TLEN
      ELSE IF (LOCAL(WHERE2).EQ.BINIT) THEN
         EDUE2B=FNCVAL*ERRINB
         EPS=MAX(EPSR-(EDUE2A+EDUE2B),EDUE2A+EDUE2B)*DELTA/TLEN
      END IF
      EPS=MAX(EPS,(EPSMIN/AACUM)**RELOBT*AACUM)
2830  J=WHERE2+2*(WHERE-4)
      GO TO (280,290,350,350), J
C
C     COMPUTE F USING FORWARD OR REVERSE COMMUNICATION.
C
2840  IF (ABS(ABSCIS-LOCAL(4-WHERE)).GT.ESMALL) GO TO 2850
      WHERE=-WHERE
      GO TO 2180
2850  IF (NFMAX.LE.0) GO TO 2860
      IF (NFEVAL.LT.NFMAX) GO TO 2860
      WHERE=7

      IFLAG(1)=5
      GO TO 2345
2860  NFEVAL=NFEVAL+1
      T=1.0d0
      IF (NSUB.NE.0) THEN
c                 Put ABSCIS in user coordinates
         X1=TA/TB
         T=(ABSCIS-TA)/TB
c?       ABSCIS=TA*(1.0d0+X1)+ABSCIS*(T-X1)
         ABSCIS=TA+ABSCIS*T+X1*(TA-ABSCIS)
         IF (NSUB.NE.2) THEN
            F2=(ABSCIS-TA)/TB
c?          ABSCIS=TA*(1.0d0+X1)+ABSCIS*(F2-X1)
            ABSCIS=TA+ABSCIS*F2+X1*(TA-ABSCIS)
            T=T*F2
            T=T+T
         END IF
         T=T+T
      END IF
      IF (ABSCIS.EQ.ALOCAL) THEN
c##      F2=DINTSM(ABSCIS)/EDELM2
c##      ABSCIS=ABSCIS+SIGN(F2,TB)
         abscis=abscis+sign(max(enzer,emeps*abs(abscis)),tb)
         X=ABSCIS
      ELSE IF (ABSCIS.EQ.BLOCAL) THEN
c##      F2=DINTSM(ABSCIS)/EDELM2
c##      ABSCIS=ABSCIS-SIGN(F2,TB)
         abscis=abscis-sign(max(enzer,emeps*abs(abscis)),tb)
         X=ABSCIS
      END IF
      WORK(KDIM)=ABSCIS
      IF (REVERS.NE.0) RETURN
      CALL DINTF (ANSWER,WORK,IFLAG(1))
2900  FNCVAL=T*ANSWER
      GO TO (90,150,160,2820,2820,350,2860), WHERE
C
      END
