      SUBROUTINE DINTM (NDIMI,ANSWER,WORK,NWORK,IOPT)
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 DINTM  Krogh  Initialized lots of variables.
C>> 1996-03-31 DINTM  Krogh  Removed unused variable in common.
c>> 1995-11-20 DINTM  Krogh  Converted from SFTRAN to Fortran 77.
C>> 1994-11-14 DINTM  Krogh  Declared all vars.
c>> 1994-10-19 DINTM  Krogh  Changes to use M77CON
c>> 1994-07-07 DINTM  Snyder set up for CHGTYP.
C>> 1993-05-18 DINTM  Krogh -- Changed "END" to "END PROGRAM"
C>> 1994-05-02 DINTM  Snyder corrected some comments
C>> 1991-09-20 DINTM  Krogh converted '(1)' dimensioning to '(*)'.
C>> 1987-11-19 DINTM  Snyder  Initial code.
c
c--D replaces "?": ?INTA,?intc,?intec,?INTF,?INTM,?INTMA,?INTNC,?INTOP
C
C     ******************************************************************
C
C     THIS SUBROUTINE ATTEMPTS TO CALCULATE THE INTEGRAL OF A FUNCTION
C     PROVIDED BY THE USER.  THE FUNCTION IS PROVIDED BY THE USER VIA A
C     SUBROUTINE REFERENCED BY CALL DINTF(F,X,IFLAG) OR BY REVERSE
C     COMMUNICATION.  ALL ABSCISSAE ARE WITHIN THE REGION OF
C     INTEGRATION.
C
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 IN THE REGION OF
C     INTEGRATION, AND NON-INTEGRABLE SINGULARITIES NEAR THIS REGION.
C     (INCLUDING COMPLEX VALUES).  THE EFFICIENCY AND RELIABILITY
C     OF INTEGRATING SUCH FUNCTIONS MAY BE GREATLY IMPROVED BY MANUALLY
C     SUBDIVIDING THE REGION OF INTEGRATION AT THE DISCONTINUITY,
C     SINGULARITY OR CLOSEST POINT TO A 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 NDIMI   IS THE NUMBER OF DIMENSIONS OF INTEGRATION.
C ANSWER  IS THE ESTIMATE OF THE INTEGRAL.  WHEN REVERSE COMMUNICATION
C         IS SPECIFIED IT IS USED TO PASS FUNCTION VALUES FROM
C         THE USER PROGRAM TO DINTA OR DINTMA.
      DOUBLE PRECISION ANSWER
C WORK    CONTAINS THE LIMITS, WORKING SPACE AND MAY BE REFERENCED BY
C         THE OPTION VECTOR (SEE IOPT BELOW).  WHEN THE INTEGRATION
C         IS COMPLETE, WORK(1) CONTAINS THE ESTIMATED ABSOLUTE ERROR.
C         WHEN REVERSE COMMUNICATION IS SPECIFIED, WORK(1) IS USED TO
C         PASS ABSCISSAE FROM DINTA OR DINTMA TO THE USER PROGRAM.
C         WORK(1), ..., WORK(NDIMI) ARE COMPONENTS OF THE VECTOR OF
C         ABSCISSAE.  WORK(NDIMI+1), ..., WORK(2*NDIMI) ARE LOWER
C         LIMITS.  WORK(2*NDIMI+1), ..., WORK(3*NDIMI) ARE UPPER LIMITS.
C         THE ABSCISSAE AND LIMITS ARE STORED INNERMOST FIRST.  DINTF
C         IS CALLED TO REQUEST THE LIMITS OF EVERY DIMENSION, BUT
C         CONSTANT LIMITS (THE LIMITS OF THE OUTER DIMENSION ARE ALWAYS
C         CONSTANT) MAY INSTEAD BE STORED BEFORE INTEGRATION BEGINS.
C         WORK(3*NDIMI+1), ..., WORK(3*NDIMI+KWORK*(NDIMI-1)) ARE
C         WORKING STORAGE, WHERE KWORK DEPENDS ON THE MACHINE AND THE
C         PRECISION OF THE PROGRAM VERSION.  KWORK MUST SPECIFY ENOUGH
C         STORAGE FOR 4 VARIABLES WHICH ARE ALWAYS DOUBLE PRECISION,
C         139 VARIABLES WHICH ARE DOUBLE PRECISION IN THE DOUBLE
C         PRECISION VERSION, 30 VARIABLES WHICH ARE DOUBLE PRECISION
C         IN THE DOUBLE PRECISION PROGRAM IF THE LARGEST REPRESENTABLE
C         NUMBERS IN SINGLE AND DOUBLE PRECISION ARE DIFFERENT BY MORE
C         THAN A FACTOR OF 4, 29 INTEGER VARIABLES AND 11 LOGICAL
C         VARIABLES.  REMEMBER TO CONSIDER THE DATA TYPE OF WORK WHEN
C         SPECIFYING THE VALUE OF KWORK.  EXAMPLES ARE
C-- Begin mask code changes
C           IF PRECISION IS REAL, AND INTEGER AND LOGICAL VARIABLES
C           OCCUPY THE SAME SPACE AS REAL VARIABLES, KWORK=217.
C           IF PRECISION IS REAL, AND INTEGER AND LOGICAL VARIABLES
C           OCCUPY HALF AS MUCH SPACE AS REAL VARIABLES, KWORK=?
C           IF PRECISION IS PARTLY DOUBLE, AND INTEGER AND LOGICAL
C           VARIABLES OCCUPY THE SAME SPACE AS REAL VARIABLES,
C           KWORK=?
C           IF PRECISION IS PARTLY DOUBLE, AND INTEGER AND LOGICAL
C           VARIABLES OCCUPY HALF AS MUCH SPACE AS REAL VARIABLES,
C           KWORK=?
C           IF PRECISION IS DOUBLE, AND INTEGER AND LOGICAL VARIABLES
C           OCCUPY THE SAME SPACE AS REAL VARIABLES, KWORK=193.
C           IF PRECISION IS DOUBLE, AND INTEGER AND LOGICAL VARIABLES
C           OCCUPY HALF AS MUCH SPACE AS REAL VARIABLES, KWORK=?
C           FOR PORTABILITY, ASSUME THE WORST CASE.  THAT IS, IF THE
C           PRECISION IS SINGLE, KWORK=217.  IF THE PRECISION IS DOUBLE,
C           KWORK=193.
C         WORK(3*NDIMI+KWORK*(NDIMI-1)+1), ... ARE AVAILABLE FOR
C         REFERENCE THROUGH THE OPTION VECTOR.
C-- End mask code changes
      DOUBLE PRECISION WORK(*)
C NWORK   IS THE NUMBER OF ELEMENTS ALLOCATED FOR WORK.
C         NWORK MUST BE GREATER THAN 3*NDIMI+KWORK*(NDIMI-1).
C IOPT    IS A VECTOR OF INTEGERS USED TO RETURN THE STATUS AND TO
C         SPECIFY OPTIONS.
      INTEGER IOPT(*)
C         IOPT(1) IS USED TO RETURN A STATUS INDICATOR TO THE USER.
C         VALUES OF THIS INDICATOR ARE
C         -NDIM   - NORMAL TERMINATION, EITHER THE ABSOLUTE OR THE
C                   RELATIVE ERROR CRITERIA ARE SATISFIED.
C         -NDIM-1 - NORMAL TERMINATION, NEITHER THE ABSOLUTE NOR THE
C                   RELATIVE ERROR CRITERIA ARE SATISFIED, BUT THE
C                   ERROR RELATIVE TO THE OBTAINABLE PRECISION
C                   PRECISION CRITERION IS SATISFIED.
C         -NDIM-2 - NORMAL TERMINATION, BUT NONE OF THE ERROR CRITERIA
C                   ARE SATISFIED.
C         -NDIM-3 - NOT ENOUGH SPACE IN WORK, NWORK IS TOO SMALL.
C         -NDIM-4 - BAD IOPT VALUE.
C         -NDIM-5 - TOO MANY FUNCTION VALUES NEEDED.
C         -NDIM-KDIM-5 - APPARENT NON-INTEGRABLE SINGULARITY IN
C                   DIMENSION KDIM.  ANSWER CONTAINS THE APPROXIMATE
C                   ABSCISSA OF THE SINGULARITY IN THE KDIM-TH
C                   DIMENSION, WORK(KDIM+1) THROUGH WORK(NDIM)
C                   CONTAIN THE ABSCISSAE OF EXTERIOR INTEGRALS.
C         -2*NDIM-KDIM-5 - INCORRECT INNER INTEGRAL DIMENSIONALITY
C                   SPECIFIED DURING EVALUATION OF INTEGRAL AT
C                   DIMENSION KDIM, THAT IS, IABS(IFLAG(1)) = KDIM.
C         ENTRIES IN IOPT STARTING WITH IOPT(2) ARE DESCRIBED BELOW.
C         IOPT(I)   ENTRY IN IOPT(I) MEANS
C           0       NO MORE OPTIONS, BEGIN INTEGRATION.
C           1       IOPT(I+1) CONTAINS THE UNIT NUMBER FOR OUTPUT.
C           2       IOPT(I+1) IS AN NDIMI DIGIT INTEGER, WHERE
C                   EACH DIGIT IS THE DIAGNOSTIC PRINT LEVEL FOR
C                   ONE DIMENSION.  THE LOW ORDER DIGIT IS THE PRINT
C                   LEVEL FOR THE INNER DIMENSION.
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           3       WORK(IOPT(I+1)) CONTAINS THE ABSOLUTE TOLERANCE,
C                   WORK(IOPT(I+1)+1) CONTAINS TOLERANCE RELATIVE TO
C                   THE LOCALLY OBTAINABLE PRECISION, AND
C                   WORK(IOPT(I+1)+2) CONTAINS TOLERANCE RELATIVE TO
C                   THE VALUE OF THE INTEGRAL.  THE TOLERANCE RELATIVE
C                   TO THE LOCALLY OBTAINABLE PRECISION IS SPECIFIED AS
C                   THE FRACTION OF LOCALLY OBTAINABLE DIGITS THAT ARE
C                   PERMITTED TO BE WRONG, AND IS INTERNALLY BOUNDED
C                   BETWEEN 0.0 AND 1.0.  IF ANY OF THE ERROR CONTROL
C                   CRITERIA ARE SATISFIED, THE ANSWER IS ACCEPTED.  IF
C                   THIS OPTION IS NOT USED, THE ABSOLUTE AND RELATIVE
C                   TOLERANCES ARE SET TO ZERO, AND THE TOLERANCE
C                   RELATIVE TO THE LOCALLY OBTAINABLE PRECISION IS SET
C                   TO 0.25.
C           4       WORK(IOPT(I+1)) CONTAINS ABSOLUTE ERROR COMMITTED
C                   COMPUTING F.  WORK(IOPT(I+1)) MAY BE CHANGED DURING
C                   THE INTEGRATION.
C           5       WORK(IOPT(I+1)) CONTAINS RELATIVE ERROR EXPECTED TO
C                   BE COMMITTED COMPUTING F.  CHANGES TO THIS VALUE
C                   DURING INTEGRATION WILL NOT BE DETECTED.
C           6       USE REVERSE COMMUNICATION -
C                   CALL DINTM (NDIMI,ANSWER,WORK,NWORK,IOPT)
C                 1 CALL DINTMA (ANSWER,WORK,IOPT(1))
C                   IF (IOPT(1)) 3,2,4
C                 2 COMPUTE THE INNERMOST INTEGRAND -
C                        ANSWER=F(WORK(K)), K=1,NDIMI
C                   THEN EITHER
C                        GO TO 1
C                   OR (FOR BETTER EFFICIENCY)
C                        CALL DINTA (ANSWER,WORK,IOPT(1))
C                        IF (IOPT(1)) 1,2,5
C                   VALUES OF IOPT(1) PRODUCED BY DINTA AND DINTMA
C                   HAVE DIFFERENT MEANINGS.
C                 3 IF (IOPT(1)+NDIMI.LE.0) GO TO 6
C                   THE INTEGRAL OVER THE (-IOPT(1))TH DIMENSION IS
C                   CONTAINED IN ANSWER.  IF A TRANSFORMATION OF
C                   THE INTEGRAL IS NOT NECESSARY BEFORE IT IS USED
C                   AS AN INTEGRAND FOR THE (1-IOPT(1))TH DIMENSION
C                        GO TO 1
C                   ELSE COMPUTE ANSWER=G(ANSWER,WORK(K)),
C                   K=1-IOPT(1),NDIMI AND G IS THE TRANSFORMATION.
C                   COMPUTE ALSO WORK(1)=PARTIAL DERIVATIVE OF G
C                   WITH RESPECT TO ANSWER THEN
C                        GO TO 1
C                 4 COMPUTE WORK(NDIMI+IOPT(1))=A(WORK(K)),
C                   K=IOPT(1)+1,NDIMI, WORK(2*NDIMI+IOPT(1))=B(WORK(K))
C                   K=IOPT(1)+1,NDIM, WHERE A IS THE LOWER LIMIT OF
C                   INTEGRATION FOR THE (IOPT(1))TH DIMENSION, AND B IS
C                   THE UPPER LIMIT.  IF A TRANSFORMATION AS DISCUSSED
C                   ABOVE IS TO BE APPLIED WHEN THE INTEGRATION OF THIS
C                   DIMENSION IS COMPLETE, COMPUTE WORK(1)=AN ESTIMATE
C                   OF THE UPPER BOUND OF THE PARTIAL DERIVATIVE OF G
C                   WITH RESPECT TO THE INTEGRAL.  IF G IS LINEAR IN
C                   THE INTEGRAL, THIS DERIVATIVE CAN BE COMPUTED
C                   WITHOUT KNOWING THE INTEGRAL, RATHER THAN
C                   ESTIMATED.  IF THERE IS A SINGULARITY IN THE
C                   INTEGRAND FOR THIS DIMENSION (PERHAPS INTRODUCED
C                   BY THE BOUNDARY OR THE TRANSFORMATION G), SET
C                   IOPT(1) TO THE LOCATION IN WORK OF THE ABSCISSA IN
C                   THIS DIMENSION OF THE SINGULARITY, AND WITH THE
C                   SIGN SET AS DESCRIBED FOR OPTION 11 BELOW.  THE
C                   MAGNITUDE OF IOPT(1) MUST BE GREATER THAN NDIM FOR
C                   THIS REQUEST TO BE DETECTED.  IN PARTICULAR, IF
C                   THERE IS A SINGULARITY AT ONE OF THE LIMITS, SET
C                   IOPT(1) TO IOPT(1)+NDIM OR -(IOPT(1)+2*NDIM).  THEN
C                        GO TO 1.
C                 5 IOPT(1)=-(IOPT(1)+NDIM)
C                 6 CONTINUE
C                   AT THIS POINT, IERR CONTAINS THE STATUS DESCRIBED
C                   BELOW, FOR OPTION 8.
C           7       SET MINIMUM INDEX OF QUADRATURE FORMULA TO USE
C                   BEFORE SUBDIVISION TO IOPT(I+1).
C           8       NO EFFECT.  IOPT(I+1) MAY BE USED TO PASS
C                   INFORMATION TO DINTF.  INCREMENT I BY 2.
C           9       IOPT(I+1) IS THE MAXIMUM NUMBER OF FUNCTION VALUES
C                   TO USE TO INTEGRATE THE ENTIRE PROBLEM.  IF
C                   IOPT(I+1) .LE. 0, THE NUMBER OF FUNCTION VALUES
C                   IS NOT CONTROLLED.
C          10       IOPT(I+1) IS USED TO RETURN THE NUMBER OF FUNCTION
C                   VALUES USED TO INTEGRATE THE ENTIRE PROBLEM.
C          11       WORK(IABS(IOPT(I+1))) IS THE LOCATION OF A
C                   SINGULARITY OR DISCONTINUITY IN THE OUTERMOST
C                   INTEGRAL.  IF THE LOCATION IS INSIDE THE INTERVAL,
C                   THE INTERVAL IS IMMEDIATELY SUBDIVIDED.  IF
C                   IOPT(I+1) .GT. 0, A T**2 SUBSTITUTION BASED AT
C                   WORK(...) WILL BE USED.  IF IOPT(I+1) .LT. 0,
C                   A T**4 SUBSTITUTION WILL BE USED.  IF IOPT(I+1)
C                   .EQ. 0, NO SUBSTITUTION WILL BE USED.  TO NOTIFY
C                   THE PROGRAM OF SINGULARITIES OR DISCONTINUITIES IN
C                   INTERIOR INTEGRANDS, SET IFLAG NEGATIVE WHEN
C                   ASKED TO COMPUTE LIMITS.  SEE THE DESCRIPTION
C                   OF REVERSE COMMUNICATION ABOVE, OR THE DESCRIPTION
C                   OF DINTF BELOW.
C          12       WORK(IOPT(I+1)) CONTAINS THE ABSOLUTE ERROR IN THE
C                   LOWER LIMIT OF THE CURRENT DIMENSION.  THE ERROR IN
C                   THE UPPER LIMIT IS IN WORK(IOPT(I+1)+1).
C          13       IOPT(I+1) IS THE LOCATION IN IOPT IN WHICH THE USER
C                   IS TO STORE NON-STANDARD CHANGES TO THE DIMENSION OF
C                   INNER INTEGRALS.  THE DEFAULT VALUE IS 1, ALLOWING
C                   SUCH NONSTANDARD CHANGES TO BE STORED IN THE CELL OF
C                   IFLAG (=IOPT) USUALLY USED FOR COMMUNICATION WITH
C                   DINTF, IN THE CASE WHEN IT IS NOT SIMULTANEOUSLY
C                   NECESSARY TO INDICATE A SINGULARITY.  SEE THE
C                   DESCRIPTION OF THE INTERFACE TO DINTF BELOW.
C
C     ALL OPTIONS ARE SET TO THEIR NOMINAL VALUES BEFORE THE OPTION
C     VECTOR IS PROCESSED.
C
C    *****     DINTF     ***********************************************
C
C     DINTF IS REFERENCED VIA CALL DINTF (ANSWER,WORK,IFLAG).
C     WORK AND ANSWER ARE AS DESCRIBED ABOVE.  VALUES OF IFLAG ARE
C     IFLAG .LT. 0
C                    THE INTEGRAL OVER THE (-IFLAG)TH DIMENSION IS
C                    CONTAINED IN ANSWER.  IF A TRANSFORMATION OF
C                    THE INTEGRAL IS NOT NECESSARY BEFORE IT IS USED
C                    AS AN INTEGRAND FOR THE (1-IFLAG)TH DIMENSION
C                         RETURN
C                    ELSE IF THE INTEGRAND DEPENDS ON ONLY ONE INNER
C                    INTEGRAL COMPUTE ANSWER=G(ANSWER,WORK(K)),
C                    ALSO MULTIPLY WORK(1) BY THE PARTIAL DERIVATIVE
C                    OF G WITH RESPECT TO ANSWER THEN
C                         RETURN
C                    ELSE IF THE INTEGRAND DEPENDS ON MORE INTEGRALS
C                    SAVE ANSWER AND WORK(1), AND SET IFLAG(IXKDIM) TO
C                    THE DIMENSION OF INTEGRAL NEXT TO BE COMPUTED
C                    (IXKDIM IS THE VALUE SPECIFIED BY OPTION 12, OR 1
C                    IF OPTION 12 IS NOT SELECTED), THEN
C                         RETURN
C                    ELSE (ALL INNER INTEGRALS UPON WHICH THE INNER
C                    INTEGRAND DEPENDS HAVE BEEN EVALUATED) CALCULATE
C                    ANSWER=G(ANSWER,WORK(K),SAVED INTEGRALS)
C                    K=1-IFLAG, NDIMI (G IS THE TRANSFORMATION),
C                    MULTIPLY WORK(1) BY THE PARTIAL DERIVATIVE OF G
C                    WITH RESPECT TO THE LAST INTEGRAL, AND ADD THE
C                    PRODUCT OF EACH OF THE SAVED VALUES OF WORK(1) AND
C                    PARTIAL DERIVATIVE OF G WITH RESPECT TO THE
C                    CORRESPONDING SAVED INTEGRAL ONTO WORK(1).  THAT
C                    IS, WORK(1) IS THE TOTAL ERROR IN G, CALCULATED AS
C                    THE INNER PRODUCT OF THE ERRORS IN THE INDIVIDUAL
C                    ARGUMENTS, AND THE PARTIAL DERIVATIVES OF G WITH
C                    RESPECT TO THOSE ARGUMENTS.
C     IFLAG .EQ. 0
C                    COMPUTE THE INNERMOST INTEGRAND -
C                         ANSWER=F(WORK(K)), K=1,NDIMI
C                    AND
C                         RETURN
C     IFLAG .GT. 0
C                    COMPUTE WORK(NDIMI+IFLAG)=A(WORK(K)),
C                    K=IFLAG+1,NDIMI, WORK(2*NDIMI+IFLAG)=B(WORK(K)),
C                    K=IFLAG+1,NDIM, WHERE A IS THE LOWER LIMIT OF
C                    INTEGRATION FOR THE (IFLAG)TH DIMENSION, AND B IS
C                    THE UPPER LIMIT.  IF A TRANSFORMATION AS DISCUSSED
C                    ABOVE IS TO BE APPLIED WHEN THE INTEGRATION OF THIS
C                    DIMENSION IS COMPLETE, COMPUTE WORK(1)=AN ESTIMATE
C                    OF THE UPPER BOUND OF THE PARTIAL DERIVATIVE
C                    OF G WITH RESPECT TO THE INTEGRAL.  IF G IS
C                    LINEAR IN THE INTEGRAL, THIS DERIVATIVE MAY BE
C                    COMPUTED WITHOUT KNOWING THE INTEGRAL, RATHER
C                    THAN ESTIMATED.  IF THERE IS A SINGULARITY IN
C                    THE INTEGRAND FOR THIS DIMENSION (PERHAPS
C                    INTRODUCED BY THE BOUNDARY OR THE TRANSFORMATION),
C                    SET IFLAG TO THE LOCATION IN WORK OF THE ABSCISSA
C                    IN THIS DIMENSION OF THE SINGULARITY, AND SET THE
C                    SIGN AS FOR OPTION 11 ABOVE.  NOTE THAT THE
C                    MAGNITUDE OF IFLAG MUST BE GREATER THAN NDIM FOR
C                    THIS REQUEST TO BE DETECTED.  IN PARTICULAR, IF
C                    THERE IS A SINGULARITY AT ONE OF THE LIMITS, SET
C                    THE MAGNITUDE OF IFLAG TO IFLAG+NDIM OR
C                    IFLAG+2*NDIM, WITH THE SIGN AS DESCRIBED FOR
C                    OPTION 11.  THEN
C                         RETURN.
C                    IF THE FIRST INTEGRAL TO BE COMPUTED DOES NOT HAVE
C                    DIMENSION IFLAG(1)-1, SET IFLAG(IXKDIM) TO THE
C                    DIMENSION OF INTEGRAL TO BE NEXT COMPUTED.  SEE
C                    OPTION 12 AND THE DISCUSSION OF IFLAG .LT. 0 ABOVE
C                    FOR INFORMATION ON IXKDIM.
C
C     *****    EXTERNAL REFERENCES     *********************************
C
C DINTMA  TO DO THE INTEGRATION.
C DINTOP  TO SET OPTIONS.
C
C     *****    LOCAL VARIABLES     *************************************
C
      INTEGER NDIMI, NWORK
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 For initialization of common blocks
      integer KREAL, KINT, KLOG, ITEMP
      parameter (KINT=29, KREAL=169, KLOG=11)
      logical LMOVE(KLOG)
      integer IMOVE(KINT)
      double precision RMOVE(KREAL)
      equivalence (AACUM, RMOVE)
      equivalence (DISCF, IMOVE)
      equivalence (DID1, LMOVE)
C
C
C     *****    PROCEDURES     ******************************************
C

c Initialize common blocks to avoid references to undefined variables.

      do 5 ITEMP = 1, KLOG
        LMOVE(ITEMP) = .false.
 5    continue
      do 10 ITEMP = 1, KREAL
        RMOVE(ITEMP) = 0.D0
 10   continue
      do 15 ITEMP = 1, KINT
        IMOVE(ITEMP) = 0
 15   continue

      WHEREM=0
      NFEVAL=0
      NDIM=NDIMI
      KDIM=NWORK
      X = 0.D0


C     KDIM IS TEMPORARILY USED IN DINTMA TO CHECK THE DIMENSION OF
C     WORK.
      CALL DINTOP (IOPT,WORK)
      IF (IOPT(1).EQ.0) THEN
C
C     CALL DINTMA TO DO THE INTEGRATION.
C
         IF (REVERM.EQ.0) CALL DINTMA (ANSWER,WORK,IOPT)
      ELSE
         IOPT(1)=-NDIM-4
      END IF
      RETURN
      END
