SUBROUTINE DCUTET(FUNSUB,NUMFUN,VER,NUMTET,MINPTS,MAXPTS,EPSABS,
     +                  EPSREL,LENVER,NW,RESTAR,RESULT,ABSERR,NEVAL,
     +                  IFAIL,WORK,IWORK)
C***BEGIN PROLOGUE DCUTET
C***DATE WRITTEN   900220   (YYMMDD)
C***REVISION DATE  900307   (YYMMDD)
C***CATEGORY NO. H2B2A1/H2B2A2
C***KEYWORDS QUADRATURE, THREE DIMENSIONAL, ADAPTIVE, CUBATURE
C***AUTHORS
C            Jarle Berntsen, The Computing Centre,
C            University of Bergen, Thormohlens gt. 55,
C            N-5008 Bergen, Norway
C            Phone..  47-5-544055
C            Email..  jarle@eik.ii.uib.no
C
C            Ronald Cools, Dept. of Computer Science,
C            Katholieke Universiteit Leuven, Celestijnenlaan 200A,
C            B-3030 Heverlee, Belgium
C            Phone..  32-16-201015 (3562)
C            Email..  ronald@cs.kuleuven.ac.be
C
C            Terje O. Espelid, Department of Informatics,
C            University of Bergen, Thormohlens gt. 55,
C            N-5008 Bergen, Norway
C            Phone..  47-5-544180
C            Email..  terje@eik.ii.uib.no
C
C***PURPOSE  To compute the three-dimensional integral over a region
C            consisting of NUMTET tetrahedrons.
C***DESCRIPTION
C
C            The routine calculates an approximation to a given
C            vector of definite integrals
C
C            I  I  I (F ,F ,...,F )      DX(3)DX(2)DX(1),
C                      1  2      NUMFUN
C
C            where the region of integration is a collection of
C            NUMTET tetrahedrons and
C            where F = F (X(1),X(2),X(3)), J = 1,2,...,NUMFUN.
C                   J   J
C
C
C            A globally adaptive strategy is applied in order to
C            compute approximations RESULT(K)
C            hopefully satisfying, for each component of I, the
C            following claim for accuracy:
C            ABS(I(K)-RESULT(K)).LE.MAX(EPSABS,EPSREL*ABS(I(K)))
C
C            DCUTET is a driver for the integration routine
C            DADTET.
C
C            DADTET repeatedly
C            subdivides the tetrahedrons with greatest estimated  errors
C            and estimates the integrals and the
C            errors over the new subtetrahedrons
C            until the error request
C            is met or MAXPTS function evaluations have been used.
C
C            A 43 point integration rule
C            with all evaluation points inside the tetrahedron
C            is applied. The rule has polynomial degree 8.
C
C            If the values of the input parameters EPSABS
C            or EPSREL are selected great enough,
C            an integration rule is applied over each tetrahedron and
C            the results are added up to give the approximations
C            RESULT(K). No further subdivision
C            of the tetrahedrons will then be applied.
C
C            When DCUTET computes estimates to a vector of
C            integrals, all components of the vector are given
C            the same treatment. That is, I(F ) and I(F ) for
C                                            J         K
C            J not equal to K, are estimated with the same
C            subdivision of the region of integration.
C            For integrals with enough similarity, we may save
C            time by applying DCUTET to all integrands in one call.
C            For integrals that varies continuously as functions of
C            some parameter, the estimates produced by DCUTET will
C            also vary continuously when the same subdivision is
C            applied to all components. This will generally not be
C            the case when the different components are given
C            separate treatment.
C
C            On the other hand this feature should be used with
C            caution when the different components of the integrals
C            require clearly different subdivisions.
C
C   ON ENTRY
C
C     FUNSUB Externally declared subroutine for computing
C            all components of the integrand at the given
C            evaluation point.
C            It must have parameters (X,NUMFUN,FUNVLS)
C            Input parameters:
C              X(1)   The x-coordinate of the evaluation point.
C              X(2)   The y-coordinate of the evaluation point.
C              X(3)   The z-coordinate of the evaluation point.
C              NUMFUN Integer that defines the number of
C                     components of I.
C            Output parameter:
C              FUNVLS Real array of dimension NUMFUN
C                     that defines NUMFUN components of the integrand.
C
C     NUMFUN Integer.
C            Number of components of the integral.
C     VER    Real array of dimension (3,4,LENVER).
C            VER(1,K,L), VER(2,K,L) and VER(3,K,L) are the x, y and z
C            coordinates respectively of vertex K in tetrahedron L.
C            On entry VER(*,*,L) must contain the vertices of the
C            NUMTET user specified tetrahedrons for L = 1,2,...,NUMTET.
C     NUMTET Integer.
C            The number of tetrahedrons.
C     MINPTS Integer.
C            Minimum number of function evaluations.
C     MAXPTS Integer.
C            Maximum number of function evaluations.
C            The number of function evaluations over each subregion
C            is 43.
C
C            MAXPTS >= 43*NUMTET and MAXPTS >= MINPTS
C
C     EPSABS Real.
C            Requested absolute error.
C     EPSREL Real.
C            Requested relative error.
C     LENVER Integer.
C            Defines the length of the array VER.
C
C            We let
C            MAXSUB denote the maximum allowed number of subregions
C            for the given value of MAXPTS.
C            MAXSUB = 7*(MAXPTS-43*NUMTET)/(8*43) + NUMTET
C            LENVER should be greater or equal to MAXSUB.
C
C     NW     Integer.
C            Defines the length of the working array WORK.
C
C            We let
C            MAXSUB denote the maximum allowed number of subregions
C            for the given values of MAXPTS.
C            MAXSUB = 7*(MAXPTS-43*NUMTET)/(8*43) + NUMTET
C            NW should then be greater or equal to
C            MAXSUB*(2*NUMFUN+1) + MAX(8*6*MDIV,6*NUMTET)*NUMFUN + 1
C            where MDIV is the number of tetrahedrons that are divided
C            in each subdivision step.
C            MDIV is default set internally in DCUTET equal to 1.
C            For efficient execution on parallel computers
C            with NPROC processors MDIV should be set equal to
C            the smallest integer such that MOD(8*MDIV,NPROC) = 0.
C
C     RESTAR Integer.
C            If RESTAR = 0, this is the first attempt to compute
C            the integral.
C            If RESTAR = 1,
C            then we restart a previous attempt.
C            In this case the only parameters for DCUTET that may
C            be changed (with respect to the previous call of DCUTET)
C            are MINPTS, MAXPTS, EPSABS, EPSREL and RESTAR.
C
C   ON RETURN
C
C     RESULT Real array of dimension NUMFUN.
C            Approximations to all components of the integral.
C     ABSERR Real array of dimension NUMFUN.
C            Estimates of absolute errors.
C     NEVAL  Integer.
C            Number of function evaluations used by DCUTET.
C     IFAIL  Integer.
C            IFAIL = 0 for normal exit.
C
C              ABSERR(K) <=  EPSABS or
C              ABSERR(K) <=  ABS(RESULT(K))*EPSREL with MAXPTS or less
C              function evaluations for all values of K,
C              1 <= K <= NUMFUN .
C
C            IFAIL = 1 if MAXPTS was too small for DCUTET
C              to obtain the required accuracy. In this case DCUTET
C              returns values of RESULT with estimated absolute
C              errors ABSERR.
C            IFAIL = 2 if NUMFUN is less than 1.
C            IFAIL = 3 if the volume of one of the initially given
C                      tetrahedrons is zero.
C            IFAIL = 4 if MAXPTS is less than 43*NUMTET.
C            IFAIL = 5 if MAXPTS is less than MINPTS.
C            IFAIL = 6 if EPSABS <= 0 and EPSREL <= 0.
C            IFAIL = 7 if LENVER is less than MAXSUB.
C            IFAIL = 8 if NW is less than MAXSUB*(2*NUMFUN+1) +
C                      NUMFUN*MAX(8*6*MDIV,6*NUMTET) + 1.
C            IFAIL = 9 if unlegal RESTAR.
C     VER    Real array of dimension (3,4,LENVER).
C            On exit VER
C            contains the vertices of the tetrahedrons used to produce
C            the approximations to the integrals.
C     WORK   Real array of dimension NW.
C            Used as working storage.
C            WORK(NW) = NSUB, the number of subregions in the data
C            structure.
C            Let WRKSUB=(NW-1-NUMFUN*7*MAX(8*MDIV,NUMTET))/
C                        (2*NUMFUN+1).
C            WORK(1),...,WORK(NUMFUN*WRKSUB) contain
C              the estimated components of the integrals over the
C              subregions.
C            WORK(NUMFUN*WRKSUB+1),...,WORK(2*NUMFUN*WRKSUB) contain
C              the estimated errors over the subregions.
C            WORK(2*NUMFUN*WRKSUB+1),...,WORK(2*NUMFUN*
C              WRKSUB+WRKSUB) contain the greatest errors
C              in each subregion.
C            WORK((2*NUMFUN+1)*WRKSUB),...,WORK(NW - 1) is used as
C              temporary storage in DADTET.
C     IWORK  Integer array of dimension LENVER + MDIV.
C            Used as working storage.
C
C           SAMPLE PROGRAM
C   The following program will approximate the integral of
C                     exp(x*x+y*y+z*z)
C   over the tetrahedron (0.,0.,0.),(1.,0.,0.),(0.,1.,0.),(0.,0.,1.)
C   and print out the values of the estimated integral, the estimated
C   error, the number of function evaluations and IFAIL.
C     PROGRAM DTEST1
C     DOUBLE PRECISION VER(3,4,50),EPSABS,EPSREL,RESULT(1),ABSERR(1),
C    +                 WORK(1000)
C     INTEGER NUMFUN,NUMTET,MINPTS,MAXPTS,LENVER,NW,RESTAR,NEVAL,
C    +        IFAIL,IWORK(51)
C     EXTERNAL F
C     VER(1,1,1) = 0.D0
C     VER(1,2,1) = 1.D0
C     VER(1,3,1) = 0.D0
C     VER(1,4,1) = 0.D0
C     VER(2,1,1) = 0.D0
C     VER(2,2,1) = 0.D0
C     VER(2,3,1) = 1.D0
C     VER(2,4,1) = 0.D0
C     VER(3,1,1) = 0.D0
C     VER(3,2,1) = 0.D0
C     VER(3,3,1) = 0.D0
C     VER(3,4,1) = 1.D0
C     EPSABS = 0.D0
C     EPSREL = 1.D-5
C     NUMFUN = 1
C     NUMTET = 1
C     MINPTS = 0
C     MAXPTS = 1000
C     LENVER = 50
C     NW = 1000
C     RESTAR = 0
C     CALL DCUTET(F,NUMFUN,VER,NUMTET,MINPTS,MAXPTS,EPSABS,
C    +            EPSREL,LENVER,NW,RESTAR,RESULT,ABSERR,NEVAL,
C    +            IFAIL,WORK,IWORK)
C     WRITE(*,*)'RESULT=',RESULT(1)
C     WRITE(*,*)'ERROR =',ABSERR(1)
C     WRITE(*,*)'NEVAL =',NEVAL
C     WRITE(*,*)'IFAIL =',IFAIL
C     END
C     SUBROUTINE F(X,NUMFUN,FUNVLS)
C     INTEGER NUMFUN
C     DOUBLE PRECISION X(3),FUNVLS(NUMFUN)
C     FUNVLS(1) = EXP(X(1)*X(1)+X(2)*X(2)+X(3)*X(3))
C     RETURN
C     END
C
C   Output produced on a SUN 3/50
C
C       RESULT=   0.22779999698986
C       ERROR =    1.9752829160298D-06
C       NEVAL =  43
C       IFAIL =  0
C
C***LONG DESCRIPTION
C
C   The information for each tetrahedron is contained in the
C   data structures VER,  WORK and IWORK.
C   VER contains the coordinates of the tetrahedrons.
C   When passed on to DADTET, WORK is split into four
C   arrays VALUES, ERRORS, GREATE and WORK2.
C   VALUES contains the estimated values of the integrals.
C   ERRORS contains the estimated errors.
C   GREATE contains the greatest estimated error for each tetrahedron.
C   The data structures for the tetrahedrons are in DADTET organized
C   as a heap and the size of GREATE(I) defines the position of
C   tetrahedron I in the heap. The heap is maintained by the program
C   DTRTET and we use a partially ordered list of pointers, saved in
C   IWORK. When passed on to DADTET, IWORK is split into two
C   arrays LIST and VACANT. LIST is a partially ordered list
C   such that GREATE(LIST(1)) is the maximum error estimate
C   for all subtetrahedrons in our datastructure.
C   VACANT is a work space needed in th updating process of the list.
C
C   The subroutine DADTET is written for efficient execution on shared
C   memory parallel computers. On a computer with NPROC processors we
C   will in each subdivision step divide MDIV tetrahedrons, where MDIV
C   is chosen such that MOD(8*MDIV,NPROC) = 0, in totally 8*MDIV new
C   tetrahedrons.  Each processor will then compute estimates of the
C   integrals and errors over 8*MDIV/NPROC tetrahedrons in each
C   subdivision step.
C   The subroutine for estimating the integral and the error over
C   each subregion, DRLTET, uses WORK2 as a work array.
C   We must make sure that each processor writes its results to
C   separate parts of the memory, and therefore the sizes of WORK and
C   WORK2 are functions of MDIV.
C   In order to achieve parallel processing of tetrahedrons, compiler
C   directives should be placed in front of the DO 20 and the DO 200
C   loops in DADTET on machines like Alliant and CRAY.
C
C***REFERENCES
C   J.Berntsen, R. Cools and T.O.Espelid, An Algorithm for Adaptive
C   Cubature Over a Collection of Tetrahedrons, to be published.
C
C***ROUTINES CALLED DCHTET,DADTET
C***END PROLOGUE DCUTET
C
C   Parameters
C
C   MDIV   Integer.
C          MDIV is the number of tetrahedrons that are divided in
C          each subdivision step in DADTET.
C          MDIV is chosen default to 1.
C          For efficient execution on parallel computers
C          with NPROC processors MDIV should be set equal to
C          the smallest integer such that MOD(8*MDIV,NPROC) = 0.
C
      INTEGER MDIV
      PARAMETER (MDIV=1)

C
C   Global variables.
C
      EXTERNAL FUNSUB
      common/ierdcu/iero
      INTEGER NUMFUN,NUMTET,MINPTS,MAXPTS,LENVER,NW,RESTAR,NEVAL,IFAIL,
     +        IWORK(LENVER+MDIV)
      DOUBLE PRECISION VER(3,4,LENVER),EPSABS,EPSREL,RESULT(NUMFUN),
     +                 ABSERR(NUMFUN),WORK(NW)
C
C   Local variables.
C
C   MAXSUB Integer.
C          The maximum allowed number of subdivisions
C          for the given values of MAXPTS.
C   MINSUB Integer.
C          The minimum allowed number of subregions for the given
C          values of MINPTS.
C   WRKSUB Integer.
C          The maximum allowed number of subregions as a function
C          of NW, NUMFUN and MDIV. This determines the length
C          of the main work arrays.
C
      INTEGER MAXSUB,MINSUB,NSUB,LENW
      INTEGER WRKSUB,I1,I2,I3,I4,I5
C
C***FIRST EXECUTABLE STATEMENT DCUTET
C
C   Call DCHTET to compute MAXSUB and MINSUB,
C   and to check the input parameters.
C
      iero=0
      CALL DCHTET(NUMFUN,MDIV,VER,NUMTET,MINPTS,MAXPTS,EPSABS,EPSREL,
     +            LENVER,NW,RESTAR,MAXSUB,MINSUB,IFAIL)
C
C   All tests on input data passed?
C
      IF (IFAIL.NE.0) THEN
          RETURN
      END IF
C
C   Split up the work space.
C
      WRKSUB = (NW-1-NUMFUN*7*MAX(8*MDIV,NUMTET))/ (2*NUMFUN+1)
      I1 = 1
      I2 = I1 + WRKSUB*NUMFUN
      I3 = I2 + WRKSUB*NUMFUN
      I4 = I3 + WRKSUB
      I5 = I4 + 6*MAX(8*MDIV,NUMTET)*NUMFUN
C
C   On restart runs the number of subregions from the
C   previous call is assigned to NSUB.
C
      IF (RESTAR.EQ.1) THEN
          NSUB = WORK(NW)
      END IF
C
C   Compute the size of the temporary work space needed in DADTET.
C
      LENW = MAX(8*MDIV,NUMTET)*NUMFUN
      CALL DADTET(NUMFUN,MDIV,VER,NUMTET,MINSUB,MAXSUB,FUNSUB,EPSABS,
     +            EPSREL,LENVER,RESTAR,LENW,RESULT,ABSERR,NEVAL,NSUB,
     +            IFAIL,WORK(I1),WORK(I2),WORK(I3),WORK(I4),WORK(I5),
     +            IWORK(1),IWORK(1+LENVER))
      WORK(NW) = NSUB
      RETURN
C
C***END DCUTET
C
      END

      SUBROUTINE DCHTET(NUMFUN,MDIV,VER,NUMTET,MINPTS,MAXPTS,EPSABS,
     +                  EPSREL,LENVER,NW,RESTAR,MAXSUB,MINSUB,IFAIL)
C***BEGIN PROLOGUE DCHTET
C***REFER TO DCUTET
C***PURPOSE  DCHTET checks the validity of the
C            input parameters to DCUTET.
C***DESCRIPTION
C            DCHTET computes MAXSUB, MINSUB and IFAIL as
C            functions of the input parameters to DCUTET,
C            and checks the validity of the input parameters to DCUTET.
C
C   ON ENTRY
C
C     NUMFUN Integer.
C            Number of components of the integral.
C     MDIV   Integer.
C            MDIV is the number of tetrahedrons that are divided in
C            each subdivision step in DADTET.
C            MDIV is chosen default to 1.
C            For efficient execution on parallel computers
C            with NPROC processors MDIV should be set equal to
C            the smallest integer such that MOD(8*MDIV,NPROC) = 0.
C     VER    Real array of dimension (3,4,LENVER).
C            VER(1,K,L), VER(2,K,L) and VER(3,K,L) are the x, y and z
C            coordinates respectively of vertex K in tetrahedron L.
C            On entry VER(*,*,L) must contain the vertices of the
C            NUMTET user specified tetrahedrons for L = 1,2,...,NUMTET.
C     NUMTET Integer.
C            The number of tetrahedrons.
C     MINPTS Integer.
C            Minimum number of function evaluations.
C     MAXPTS Integer.
C            Maximum number of function evaluations.
C            The number of function evaluations over each subregion
C            is NUM.
C            MAXPTS >= NUM*NUMTET and MAXPTS >= MINPTS
C
C     EPSABS Real.
C            Requested absolute error.
C     EPSREL Real.
C            Requested relative error.
C     LENVER Integer.
C            Defines the length of the array VER.
C
C            We let
C            MAXSUB denote the maximum allowed number of subregions
C            for the given values of MAXPTS.
C            MAXSUB = 7*(MAXPTS-NUM*NUMTET)/(8*NUM) + NUMTET
C            LENVER should then be greater or equal to MAXSUB.
C
C     NW     Integer.
C            Defines the length of the working array WORK.
C
C            We let
C            MAXSUB denote the maximum allowed number of subregions
C            for the given values of MAXPTS.
C            MAXSUB = 7*(MAXPTS-NUM*NUMTET)/(8*NUM) + NUMTET
C            NW should then be greater or equal to
C            MAXSUB*(2*NUMFUN+1) + 6*MAX(8*MDIV,NUMTET)*NUMFUN + 1
C
C     RESTAR Integer.
C            If RESTAR = 0, this is the first attempt to compute
C            the integral.
C            If RESTAR = 1,
C            then we restart a previous attempt.
C            In this case the only parameters for DCUTET that may
C            be changed (with respect to the previous call of DCUTET)
C            are MINPTS, MAXPTS, EPSABS, EPSREL and RESTAR.
C
C
C   ON RETURN
C
C     MAXSUB Integer.
C            The maximum allowed number of subregions for the
C            given values of MAXPTS.
C     MINSUB Integer.
C            The minimum allowed number of subregions for the given
C            values of MINPTS.
C     IFAIL  Integer.
C            IFAIL = 0 for normal exit.
C            IFAIL = 2 if NUMFUN is less than 1.
C            IFAIL = 3 if the volume of one of the initially given
C                      tetrahedrons is zero.
C            IFAIL = 4 if MAXPTS is less than NUM*NUMTET.
C            IFAIL = 5 if MAXPTS is less than MINPTS.
C            IFAIL = 6 if EPSABS <= 0 and EPSREL <= 0.
C            IFAIL = 7 if LENVER is less than MAXSUB.
C            IFAIL = 8 if NW is less than MAXSUB*(2*NUMFUN+1) +
C                      NUMFUN*6*MAX(8*MDIV,NUMTET) + 1.
C            IFAIL = 9 if illegal value of RESTAR.
C
C***ROUTINES CALLED-NONE
C***END PROLOGUE DCHTET
C
C   Global variables.
C
      INTEGER NUMFUN,MDIV,NUMTET,MINPTS,MAXPTS,NW,MAXSUB,MINSUB,RESTAR,
     +        LENVER,IFAIL
      DOUBLE PRECISION VER(3,4,NUMTET),EPSABS,EPSREL

C
C   Local variables. 
C   NUM is set equal to the number of evaluation points used
C   by the cubature rule over each tetrahedron.
C
      INTEGER LIMIT,J,NUM
      DOUBLE PRECISION VOLUME
C
C***FIRST EXECUTABLE STATEMENT DCHTET
C
      IFAIL = 0
      NUM = 43
C
C   Compute MAXSUB and MINSUB.
C
      MAXSUB = 7* (MAXPTS-NUM*NUMTET)/ (8*NUM) + NUMTET
      MINSUB = 7* (MINPTS-NUM*NUMTET)/ (8*NUM) + NUMTET
      IF (MOD(7* (MINPTS-NUM*NUMTET),8*NUM).GT.0) THEN
          MINSUB = MINSUB + 1
      END IF
      MINSUB = MAX(NUMTET,MINSUB)
C
C   Check on positive NUMFUN.
C
      IF (NUMFUN.LT.1) THEN
          IFAIL = 2
          RETURN
      END IF
C
C   Check on legal volume of each subregion of integration.
C
      DO 10 J = 1,NUMTET
           VOLUME = (VER(1,2,J)-VER(1,1,J))*
     +     ((VER(2,3,J)-VER(2,1,J))*(VER(3,4,J)-VER(3,1,J)) -
     +      (VER(2,4,J)-VER(2,1,J))*(VER(3,3,J)-VER(3,1,J))) -
     +              (VER(2,2,J)-VER(2,1,J))*
     +     ((VER(1,3,J)-VER(1,1,J))*(VER(3,4,J)-VER(3,1,J)) - 
     +      (VER(1,4,J)-VER(1,1,J))*(VER(3,3,J)-VER(3,1,J))) +
     +              (VER(3,2,J)-VER(3,1,J))* 
     +     ((VER(1,3,J)-VER(1,1,J))*(VER(2,4,J)-VER(2,1,J)) -  
     +      (VER(1,4,J)-VER(1,1,J))*(VER(2,3,J)-VER(2,1,J)))
           VOLUME = VOLUME/6
          IF (VOLUME.EQ.0) THEN
              IFAIL = 3
              RETURN
          END IF
   10 CONTINUE
C
C   Check on MAXPTS >= NUM*NUMTET.
C
      IF (MAXPTS.LT.NUM*NUMTET) THEN
          IFAIL = 4
          RETURN
      END IF
C
C   Check on MAXPTS >= MINPTS.
C
      IF (MAXPTS.LT.MINPTS) THEN
          IFAIL = 5
          RETURN
      END IF
C
C   Check on legal accuracy requests.
C
      IF (EPSABS.LE.0 .AND. EPSREL.LE.0) THEN
          IFAIL = 6
          RETURN
      END IF
C
C   Check on big enough LENVER.
C
      IF (LENVER.LT.MAXSUB) THEN
          IFAIL = 7
          RETURN
      END IF
C
C   Check on big enough NW.
C
      LIMIT = MAXSUB* (2*NUMFUN+1) + 7*MAX(8*MDIV,NUMTET)*NUMFUN + 1
      IF (NW.LT.LIMIT) THEN
          IFAIL = 8
          RETURN
      END IF
C
C    Check on legal RESTAR.
C
      IF (RESTAR.NE.0 .AND. RESTAR.NE.1) THEN
          IFAIL = 9
      END IF
      RETURN
C
C***END DCHTET
C
      END

      SUBROUTINE DADTET(NUMFUN,MDIV,VER,NUMTET,MINSUB,MAXSUB,FUNSUB,
     +                  EPSABS,EPSREL,LENVER,RESTAR,LENW,RESULT,ABSERR,
     +                  NEVAL,NSUB,IFAIL,VALUES,ERRORS,GREATE,WORK2,
     +                  WORK3,LIST,VACANT)
C***BEGIN PROLOGUE DADTET
C***REFER TO DCUTET
C***DESCRIPTION Computation of integrals over regions consisting of
C            a collection of tetrahedrons.
C
C            DADTET repeatedly
C            subdivides the tetrahedrons with greatest estimated  errors
C            and estimates the integrals and the errors over the new
C            new sub-tetrahedrons until the error request is met
C            or MAXPTS function evaluations have been used.
C
C            A 43 point integration rule
C            with all evaluation points inside the tetrahedron
C            is applied. The rule has polynomial degree 8.
C
C   ON ENTRY
C
C     NUMFUN Integer.
C            Number of components of the integral.
C     MDIV   Integer.
C            MDIV is the number of tetrahedrons that are divided in
C            each subdivision step in DADTET.
C            MDIV is chosen default to 1.
C            For efficient execution on parallel computers
C            with NPROC processors MDIV should be set equal to
C            the smallest integer such that MOD(8*MDIV,NPROC) = 0.
C     VER    Real array of dimension (3,4,LENVER).
C            VER(1,K,L), VER(2,K,L) and VER(3,K,L) are the x, y and z
C            coordinates respectively of vertex K in tetrahedron L.
C            On entry VER(*,*,L) must contain the vertices of the
C            NUMTET user specified tetrahedrons for L = 1,2,...,NUMTET.
C     NUMTET Integer.
C            The number of tetrahedrons.
C     MINSUB Integer.
C            The minimum allowed number of subregions.
C     MAXSUB Integer.
C            The maximum allowed number of subregions.
C     FUNSUB Externally declared subroutine for computing
C            all components of the integrand at the given
C            evaluation point.
C            It must have parameters (X,NUMFUN,FUNVLS)
C            Input parameters:
C              X(1)   The x-coordinate of the evaluation point.
C              X(2)   The y-coordinate of the evaluation point.
C              X(3)   The z-coordinate of the evaluation point.
C              NUMFUN Integer that defines the number of
C                     components of I.
C            Output parameter:
C              FUNVLS Real array of dimension NUMFUN
C                     that defines NUMFUN components of the integrand.
C
C     EPSABS Real.
C            Requested absolute error.
C     EPSREL Real.
C            Requested relative error.
C
C     LENVER Integer.
C            Defines the length of the array VER.
C
C            We let
C            MAXSUB denote the maximum allowed number of subregions
C            for the given values of MAXPTS.
C            MAXSUB = 7*(MAXPTS-43*NUMTET)/(8*43) + NUMTET
C            LENVER should then be greater or equal to MAXSUB.
C
C     RESTAR Integer.
C            If RESTAR = 0, this is the first attempt to compute
C            the integral.
C            If RESTAR = 1,
C            then we restart a previous attempt.
C            In this case the only parameters for DCUTET that may
C            be changed (with respect to the previous call of DCUTET)
C            are MINPTS, MAXPTS, EPSABS, EPSREL and RESTAR.
C     LENW   Integer.
C            Length of the workspace WORK2.
C
C   ON RETURN
C
C     RESULT Real array of dimension NUMFUN.
C            Approximations to all components of the integral.
C     ABSERR Real array of dimension NUMFUN.
C            Estimates of absolute errors.
C     NEVAL  Integer.
C            Number of function evaluations used by DCUTET.
C     NSUB   Integer.
C            The number of tetrahedrons in the data structure.
C     IFAIL  Integer.
C            IFAIL = 0 for normal exit.
C
C              ABSERR(K) <=  EPSABS or
C              ABSERR(K) <=  ABS(RESULT(K))*EPSREL with MAXPTS or less
C              function evaluations for all values of K,
C              1 <= K <= NUMFUN .
C
C            IFAIL = 1 if MAXSUB was too small for DADTET
C              to obtain the required accuracy. In this case DADTET
C              returns values of RESULT with estimated absolute
C              errors ABSERR.
C     VALUES Real array of dimension (NUMFUN,MAXSUB).
C            The estimated components of the integrals over the
C            subregions.
C     ERRORS Real array of dimension (NUMFUN,MAXSUB).
C            The estimated errors over the subregions.
C     GREATE Real array of dimension MAXSUB.
C            The greatest errors in each subregion.
C     WORK2  Real array of dimension 6*LENW.
C            Work array used in DRLTET.
C     WORK3  Real array of dimension LENW.
C            Work array used in DRLTET.
C     LIST   Integer array used in DTRTET of dimension LENVER.
C            Is a partially sorted list, where LIST(1) is the top
C            element in a heap of subregions.
C     VACANT Integer array of dimension MDIV
C            Used as intermediate storage of pointers.
C***ROUTINES CALLED DTRTET,DRLTET
C***END PROLOGUE DADTET
C
C   Global variables.
C
      EXTERNAL FUNSUB
      common/ierdcu/iero

      INTEGER NUMFUN,MDIV,NUMTET,MINSUB,MAXSUB,LENW,RESTAR,LENVER,NEVAL,
     +        NSUB,IFAIL,LIST(LENVER),VACANT(MDIV)
      DOUBLE PRECISION VER(3,4,LENVER),EPSABS,EPSREL,RESULT(NUMFUN),
     +                 ABSERR(NUMFUN),VALUES(NUMFUN,MAXSUB),
     +                 ERRORS(NUMFUN,MAXSUB),GREATE(MAXSUB),WORK3(LENW),
     +                 WORK2(6*LENW)

C
C   Local variables.
C
C   SBRGNS is the number of stored subregions.
C   NDIV   The number of subregions to be divided in each main step.
C   POINTR Pointer to the position in the data structure where
C          the new subregions are to be stored.
C   G      is the homogeneous coordinates of the integration rule.
C   W      is the weights of the integration rule and the null rules.
C   TOP    is a pointer to the top element in the heap of subregions.
C   COUNT  is the number of times the termination criterium is
C          satisfied.
      INTEGER I,J,K,L,TOP,SIZE
      INTEGER SBRGNS,I1,I2,I3,I4,I5,I6,I7,I8
      INTEGER L1,K1
      INTEGER COUNT,GROUP,LEFT,NDIV,POINTR,INDEX,NUM
      DOUBLE PRECISION VEROLD(3,4),VERNEW(3,6),SUMVAL(2),SUMERR(2)
C
C***FIRST EXECUTABLE STATEMENT DADTET
C
      COUNT = 0
      NUM = 43
      IF (RESTAR.EQ.1) THEN
          SBRGNS = NSUB
          GO TO 110
      END IF
C
C   Initiate RESULT and ABSERR.
C
      DO 10 J = 1,NUMFUN
          RESULT(J) = 0
          ABSERR(J) = 0
   10 CONTINUE
C
C   Apply DRLTET over the NUMTET tetrahedrons.
C   This loop may be run in parallel.
C
      DO 20 I = 1,NUMTET
          L1 = 1 + (I-1)*6*NUMFUN
          K1 = 1 + (I-1)*NUMFUN
          CALL DRLTET(VER(1,1,I),NUMFUN,FUNSUB,WORK2(L1),VALUES(1,I),
     +                ERRORS(1,I),GREATE(I),WORK3(K1))
          if(iero.ne.0) return
   20 CONTINUE
      NEVAL = NUM*NUMTET
      SBRGNS = NUMTET
C
C   Add the computed values to RESULT and ABSERR.
C
      DO 40 I = 1,NUMTET
          DO 30 J = 1,NUMFUN
              RESULT(J) = RESULT(J) + VALUES(J,I)
              ABSERR(J) = ABSERR(J) + ERRORS(J,I)
   30     CONTINUE
   40 CONTINUE
C
C   Store results in heap.
C
      DO 50 I = 1,NUMTET
          INDEX = I
          CALL DTRTET(2,INDEX,GREATE,LIST,(INDEX))
   50 CONTINUE

C
C   We check for termination.
C
      IF (SBRGNS.LT.MINSUB) THEN
          GO TO 110
      END IF
      DO 60 J = 1,NUMFUN
          IF (ABSERR(J).GT.EPSREL*ABS(RESULT(J)) .AND.
     +        ABSERR(J).GT.EPSABS) THEN
              GO TO 110
          END IF
   60 CONTINUE
      IFAIL = 0
      COUNT = COUNT + 1
      GO TO 499
C
C***End initiation.
C
C***Begin loop while the error is too great,
C   and SBRGNS+7 is less than MAXSUB.
C
  110 IF (SBRGNS+7.LE.MAXSUB) THEN
C
C   If we are allowed to divide further,
C   prepare to apply basic rule over  the tetrahedrons produced
C   by dividing the
C   NDIV sub-tetrahedrons with greatest errors.
C   If MAXSUB and SBRGNS are great enough, NDIV = MDIV.
C
          NDIV = MAXSUB - SBRGNS
          NDIV = MIN(NDIV,MDIV,SBRGNS)
C
C   Divide each of the NDIV sub-tetrahedrons in eight new
C   sub-tetrahedrons, and compute integral and error over each.
C   When we pick a tetrahedron to divide in eight one of the
C   new sub-tetrahedron is stored in the original tetrahedrons'
C   position in the datastructure. Thus we get 7*NDIV new
C   elements in the datastructure after the subdivision. The size
C   of the datastructure before the sub-division
C   is stored in the variable SIZE, while SBRGNS is the size of
C   heap at any time.
          SIZE = SBRGNS
          DO 150 I = 1,NDIV
              POINTR = SIZE + 7* (NDIV+1-I)
C
C   Adjust RESULT and ABSERR. TOP is a pointer to the top of the heap.
C
              TOP = LIST(1)
              VACANT(I) = TOP
              DO 115 J = 1,NUMFUN
                  RESULT(J) = RESULT(J) - VALUES(J,TOP)
                  ABSERR(J) = ABSERR(J) - ERRORS(J,TOP)
  115         CONTINUE
C
C   Save the vertices.
C
              DO 130 L = 1,3
                  DO 135 K = 1,4
                      VEROLD(L,K) = VER(L,K,TOP)
  135             CONTINUE
  130         CONTINUE
C
C   Adjust the heap.
C
              CALL DTRTET(1,SBRGNS,GREATE,LIST,K)
C
C   Compute the eight new tetrahedrons. Store 7 in consecutive
C   positions and 1 in the vacant position.
C
              I1 = TOP
              I2 = POINTR - 6
              I3 = POINTR - 5
              I4 = POINTR - 4
              I5 = POINTR - 3
              I6 = POINTR - 2
              I7 = POINTR - 1
              I8 = POINTR
C
C   Compute the 6 division points.
C
              DO 140 L = 1,3
                  VERNEW(L,1) = (VEROLD(L,1)+VEROLD(L,2))/2
                  VERNEW(L,2) = (VEROLD(L,2)+VEROLD(L,3))/2
                  VERNEW(L,3) = (VEROLD(L,1)+VEROLD(L,3))/2
                  VERNEW(L,4) = (VEROLD(L,1)+VEROLD(L,4))/2
                  VERNEW(L,5) = (VEROLD(L,2)+VEROLD(L,4))/2
                  VERNEW(L,6) = (VEROLD(L,3)+VEROLD(L,4))/2
  140         CONTINUE
C
              DO 145 L = 1,3
                  VER(L,1,I1) = VEROLD(L,1)
                  VER(L,2,I1) = VERNEW(L,1)
                  VER(L,3,I1) = VERNEW(L,3)
                  VER(L,4,I1) = VERNEW(L,4)
                  VER(L,1,I2) = VEROLD(L,2)
                  VER(L,2,I2) = VERNEW(L,1)
                  VER(L,3,I2) = VERNEW(L,2)
                  VER(L,4,I2) = VERNEW(L,5)
                  VER(L,1,I3) = VEROLD(L,3)
                  VER(L,2,I3) = VERNEW(L,2)
                  VER(L,3,I3) = VERNEW(L,3)
                  VER(L,4,I3) = VERNEW(L,6)
                  VER(L,1,I4) = VEROLD(L,4)
                  VER(L,2,I4) = VERNEW(L,4)
                  VER(L,3,I4) = VERNEW(L,5)
                  VER(L,4,I4) = VERNEW(L,6)
                  VER(L,1,I5) = VERNEW(L,1)
                  VER(L,2,I5) = VERNEW(L,3)
                  VER(L,3,I5) = VERNEW(L,4)
                  VER(L,4,I5) = VERNEW(L,5)
                  VER(L,1,I6) = VERNEW(L,1)
                  VER(L,2,I6) = VERNEW(L,2)
                  VER(L,3,I6) = VERNEW(L,3)
                  VER(L,4,I6) = VERNEW(L,5)
                  VER(L,1,I7) = VERNEW(L,2)
                  VER(L,2,I7) = VERNEW(L,3)
                  VER(L,3,I7) = VERNEW(L,5)
                  VER(L,4,I7) = VERNEW(L,6)
                  VER(L,1,I8) = VERNEW(L,3)
                  VER(L,2,I8) = VERNEW(L,4)
                  VER(L,3,I8) = VERNEW(L,5)
                  VER(L,4,I8) = VERNEW(L,6)
  145         CONTINUE
  150     CONTINUE
C
C   Apply basic rule.
C   This loop may be run in parallel.
C
          DO 200 I = 1,8*NDIV
              IF (I.LE.NDIV) THEN
                  INDEX = VACANT(I)
              ELSE
                  INDEX = SBRGNS + I
              END IF
              L1 = 1 + (I-1)*6*NUMFUN
              K1 = 1 + (I-1)*NUMFUN
              CALL DRLTET(VER(1,1,INDEX),NUMFUN,FUNSUB,WORK2(L1),
     +                    VALUES(1,INDEX),ERRORS(1,INDEX),GREATE(INDEX),
     +                    WORK3(K1))
              if(iero.ne.0) return
  200     CONTINUE
          NEVAL = NEVAL + 8*NDIV*NUM
C
C   Add new contributions to RESULT and ABSERR.
C
          DO 220 I = 1,8*NDIV
              IF (I.LE.NDIV) THEN
                  INDEX = VACANT(I)
              ELSE
                  INDEX = SBRGNS + I
              END IF
              DO 210 J = 1,NUMFUN
                  RESULT(J) = RESULT(J) + VALUES(J,INDEX)
                  ABSERR(J) = ABSERR(J) + ERRORS(J,INDEX)
  210         CONTINUE
  220     CONTINUE
C
C   Store results in heap.
C
          DO 250 I = 1,8*NDIV
              IF (I.LE.NDIV) THEN
                  INDEX = VACANT(I)
              ELSE
                  INDEX = SBRGNS + I
              END IF
              J = SBRGNS + I
              CALL DTRTET(2,J,GREATE,LIST,INDEX)
  250     CONTINUE
          SBRGNS = SBRGNS + 8*NDIV
C
C   Check for termination.
C
          IF (SBRGNS.LT.MINSUB) THEN
              GO TO 110
          END IF
          DO 255 J = 1,NUMFUN
              IF (ABSERR(J).GT.EPSREL*ABS(RESULT(J)) .AND.
     +            ABSERR(J).GT.EPSABS) THEN
                  GO TO 110
              END IF
  255     CONTINUE
          IFAIL = 0
          COUNT = COUNT + 1
C
C   Else we did not succeed with the
C   given value of MAXSUB.
C
      ELSE
          IFAIL = 1
      END IF
C
C   Compute more accurate values of RESULT and ABSERR. This is done
C   to reduce the effect of roundoff on final results. Large
C   intermediate sums in the computation may course large, unnecessary
C   roundoff errors. Thus recomputing the sums of errors and estimates
C   and in addition grouping the sums in three groups should remove this
C   problem.
C
  499 CONTINUE
      DO 550 J = 1,NUMFUN
          RESULT(J) = 0
          ABSERR(J) = 0
          GROUP = SBRGNS** (1.d0/3.d0)
          IF (GROUP**3.NE.SBRGNS) THEN
              GROUP = GROUP + 1
          END IF
          LEFT = SBRGNS
          DO 520 I2 = 0,SBRGNS - 1,GROUP**2
              SUMVAL(2) = 0
              SUMERR(2) = 0
              DO 510 I1 = 0,MIN(LEFT,GROUP**2) - 1,GROUP
                  SUMVAL(1) = 0
                  SUMERR(1) = 0
                  DO 500 I = 1 + I1 + I2,MIN(LEFT,GROUP) + I1 + I2
                      SUMVAL(1) = SUMVAL(1) + VALUES(J,I)
                      SUMERR(1) = SUMERR(1) + ERRORS(J,I)
  500             CONTINUE
                  LEFT = LEFT - MIN(LEFT,GROUP)
                  SUMVAL(2) = SUMVAL(2) + SUMVAL(1)
                  SUMERR(2) = SUMERR(2) + SUMERR(1)
  510         CONTINUE
              RESULT(J) = RESULT(J) + SUMVAL(2)
              ABSERR(J) = ABSERR(J) + SUMERR(2)
  520     CONTINUE
  550 CONTINUE
C
C   Check again for termination
C
      IF (IFAIL.EQ.0) THEN
          IF (COUNT.LE.1) THEN
              DO 560 J = 1,NUMFUN
                  IF (ABSERR(J).GT.EPSREL*ABS(RESULT(J)) .AND.
     +                ABSERR(J).GT.EPSABS) THEN
                      GO TO 110
                  END IF
  560         CONTINUE
          END IF
      END IF
C
      NSUB = SBRGNS
      RETURN
C
C***END DADTET
C
      END

      SUBROUTINE DTRTET(DVFLAG,SBRGNS,GREATE,LIST,NEW)
C***BEGIN PROLOGUE DTRTET
C***REFER TO DCUTET
C***PURPOSE  DTRTET maintains a heap of subregions.
C***DESCRIPTION DTRTET maintains a heap of subregions.
C            The subregions are stored in a partially sorted
C            binary tree, ordered according to the size of the
C            greatest error estimates of each subregion(GREATE).
C            The subregion with greatest error estimate is in the
C            first position of the heap.
C
C   PARAMETERS
C
C     DVFLAG Integer.
C            If DVFLAG = 1, we remove the subregion with
C            greatest error from the heap.
C            If DVFLAG = 2, we insert a new subregion in the heap.
C     SBRGNS Integer.
C            Number of subregions in the heap.
C     GREATE Real array of dimension SBRGNS.
C            Used to store the greatest estimated errors in
C            all subregions.
C     LIST   Integer array of dimension SBRGNS.
C            Used as a partially ordered list of pointers to the
C            different subregions. This list is a heap where the
C            element on top of the list is the subregion with the
C            greatest error estimate.
C     NEW    Integer.
C            Index to the new region to be inserted in the heap.
C***ROUTINES CALLED-NONE
C***END PROLOGUE DTRTET
C
C   Global variables.
C
      INTEGER DVFLAG,NEW,SBRGNS,LIST(*)
      DOUBLE PRECISION GREATE(*)
C
C   Local variables.
C
C   GREAT  is used as intermediate storage for the greatest error of a
C          subregion.
C   SUBRGN Position of child/parent subregion in the heap.
C   SUBTMP Position of parent/child subregion in the heap.
      INTEGER SUBRGN,SUBTMP
      DOUBLE PRECISION GREAT
C
C***FIRST EXECUTABLE STATEMENT DTRTET
C
C    If DVFLAG = 1, we will reduce the partial ordered list by the
C    element with greatest estimated error. Thus the element in
C    in the heap with index LIST(1) is vacant and can be used later.
C    Reducing the heap by one element implies that the last element
C    should be re-positioned.
C
      IF (DVFLAG.EQ.1) THEN
          GREAT = GREATE(LIST(SBRGNS))
          SBRGNS = SBRGNS - 1
          SUBRGN = 1
   20     SUBTMP = 2*SUBRGN
          IF (SUBTMP.LE.SBRGNS) THEN
              IF (SUBTMP.NE.SBRGNS) THEN
C
C   Find max. of left and right child.
C
                  IF (GREATE(LIST(SUBTMP)).LT.
     +                GREATE(LIST(SUBTMP+1))) THEN
                      SUBTMP = SUBTMP + 1
                  END IF
              END IF
C
C   Compare max.child with parent.
C   If parent is max., then done.
C
              IF (GREAT.LT.GREATE(LIST(SUBTMP))) THEN
C
C   Move the pointer at position subtmp up the heap.
C
                  LIST(SUBRGN) = LIST(SUBTMP)
                  SUBRGN = SUBTMP
                  GO TO 20
              END IF
          END IF
C
C   Update the pointer.
C
          IF (SBRGNS.GT.0) THEN
              LIST(SUBRGN) = LIST(SBRGNS+1)
          END IF
      ELSE IF (DVFLAG.EQ.2) THEN
C
C   If DVFLAG = 2, find the position for the NEW region in the heap.
C
          GREAT = GREATE(NEW)
          SUBRGN = SBRGNS
   40     SUBTMP = SUBRGN/2
          IF (SUBTMP.GE.1) THEN
C
C   Compare max.child with parent.
C   If parent is max, then done.
C
              IF (GREAT.GT.GREATE(LIST(SUBTMP))) THEN
C
C   Move the pointer at position subtmp down the heap.
C
                  LIST(SUBRGN) = LIST(SUBTMP)
                  SUBRGN = SUBTMP
                  GO TO 40
              END IF
          END IF
C
C    Set the pointer to the new region in the heap.
C
          LIST(SUBRGN) = NEW
      END IF
C
C***END DTRTET
C
      RETURN
      END

      SUBROUTINE DRLTET(VER,NUMFUN,FUNSUB,NULL,BASVAL,RGNERR,GREATE,
     +                  SUMVAL)
C***BEGIN PROLOGUE DRLTET
C***REFER TO DCUTET
C***PURPOSE  To compute basic integration rule values and
C            corresponding error estimates.
C***DESCRIPTION DRLTET computes basic integration rule values
C            for a vector of integrands over a tetrahedron.
C            DRLTET also computes estimates for the errors by
C            using several null rule approximations.
C   ON ENTRY
C
C   VER    Real array of dimension (3,4).
C          The coordinates of the vertices of the tetrahedron.
C          vertex i -> ( ver(1,i),ver(2,i),ver(3,i) )
C   NUMFUN Integer.
C          Number of components of the vector integrand.
C   FUNSUB Externally declared subroutine for computing
C            all components of the integrand at the given
C            evaluation point.
C            It must have parameters (X,NUMFUN,FUNVLS)
C            Input parameters:
C              X(1)   The x-coordinate of the evaluation point.
C              X(2)   The y-coordinate of the evaluation point.
C              X(3)   The z-coordinate of the evaluation point.
C              NUMFUN Integer that defines the number of
C                     components of the vector integrand.
C            Output parameter:
C              FUNVLS Real array of dimension NUMFUN
C                     that defines NUMFUN components of the integrand.
C
C   NULL   Real array of dimension (NUMFUN,6).
C          A work array.
C   SUMVAL Real array of dimension (NUMFUN).
C          A work array
C
C   ON RETURN
C
C   BASVAL Real array of dimension NUMFUN.
C          The values for the basic rule for each component
C          of the integrand.
C   RGNERR Real array of dimension NUMFUN.
C          The error estimates for each component of the integrand.
C   GREATE Real.
C          The greatest error component of the integrand.
C
C***REFERENCES Beckers M., A. Haegemans,
C              The construction of cubature formula for the tetrahedron
C              Report TW128, K.U. Leuven (1990).
C***ROUTINES CALLED
C                   DLAMCH,DORTET,FUNSUB
C***END PROLOGUE DRLTET
C
C   Parameters
C
C   ORBITS Integer
C          The number of orbits of the cubature formula and null rules
C   CRIVAL Real
C          The decision to choose the optimistic part of the error
C          estimator is based on CRIVAL
C   FACOPT Real
C          FACOPT is the safety coefficient used in the optimistic part
C          of the error estimator.
C   FACMED Real
C          FACMED is the safety coefficient used in the non-optimistic
C          part of the error estimator. FACMED is related to CRIVAL 
C          and FACOPT.
C
      INTEGER ORBITS
      DOUBLE PRECISION CRIVAL,FACOPT,FACMED
      PARAMETER (ORBITS=7)
      PARAMETER (CRIVAL=0.5,FACMED=5,FACOPT=FACMED/CRIVAL)
C
C   Global variables.
C
      INTEGER NUMFUN
      DOUBLE PRECISION VER(3,4),BASVAL(NUMFUN),RGNERR(NUMFUN),GREATE,
     +                 SUMVAL(NUMFUN),NULL(NUMFUN,6)
      EXTERNAL FUNSUB
      common/ierdcu/iero
C
C   Local variables.
C
C   K      Integer array of dimension (0:3) that contains the structure
C          parameters. K(I) = number of orbits of type I.
C   TYPE1  Real array of dimension (K(1)).
C          Contains the first homogeneous coordinate of the generators
C          of type 1
C   TYPE2  Real array of dimension (K(2)).
C          Contains the first homogeneous coordinate of the generators
C          of type 2
C   TYPE3  Real array of dimension (2,K(2)).
C          Contains the first two homogeneous coordinates of
C          the generators of type 3.
C   WEIGHT Real array of dimension (9,ORBITS).
C          The weights of the cubature formula and the null rules.
C          WEIGHT(1,1) ,..., WEIGHT(1,ORBITS) are the weights of the
C                cubature formula
C          WEIGHT(I,1) ,..., WEIGHT(I,ORBITS) for I > 1, are the weights
C                of the null rules
C
C
      INTEGER TYPE,NR,K(0:3),I,J,P
      DOUBLE PRECISION Z(3),TYPE1(3),TYPE2(1),TYPE3(2,2),NOISE,DLAMCH,
     +                 WEIGHT(7,ORBITS),VOLUME,DEG4,DEG3,DEG1,R2,R1,R,
     +                 TRES
C
C  Cubature formula of degree 8 with 43 points
C
C  Structure parameters
C
      DATA (K(I),I=0,3)/1,3,1,2/
C
C  Information for the generators
C
      DATA (TYPE1(I),I=1,3)/0.379510205167980387748057300876D+00,
     +     0.753689235068359830728182577696D+00,
     +     0.982654148484406008240470085259D+00/
      DATA (TYPE2(I),I=1,1)/0.449467259981105775574375471447D+00/
      DATA ((TYPE3(I,J),I=1,2),J=1,2)/
     +     0.506227344977843677082264893876D+00,
     +     0.356395827885340437169173969841D-01,
     +     0.736298458958971696943019005441D+00,
     +     0.190486041934633455699433285302D+00/
C
C  Weights of the cubature formula
C
      DATA (WEIGHT(1,I),I=1,ORBITS)/-
     +     0.123001131951839495043519102752D+00,
     +     0.855018349372014074906384482699D-01,
     +     0.118021998788034059253768205083D-01,
     +     0.101900465455732427902646736855D-02,
     +     0.274781029468036908044610867719D-01,
     +     0.342269148520915110408153517904D-01,
     +     0.128431148469725555789001180031D-01/
C
C  Weights of the null rule of degree 5
C
      DATA (WEIGHT(2,I),I=1,ORBITS)/0.211921237628032658308230999090D+00
     +     ,-0.660207516445726284649283745987D-01,
     +     0.225058824086711710443385047042D-01,
     +     -0.375962972067425589765730699401D-03,
     +     0.710066020561055159657284834784D-02,
     +     0.156515256061747694921427149028D-02,
     +     -0.814530839643584660306807872526D-02/
C
C  Weights of null rule of degree 4
C
      DATA (WEIGHT(3,I),I=1,ORBITS)/-
     +     0.508105488137100551376844924797D-01,
     +     0.104596681151665328209751420525D-01,
     +     0.927471438532788763594989973184D-01,
     +     0.210489990008917994323967321174D-02,
     +     0.379184172251962722213408547663D-01,
     +     -0.111747242913563605790923001557D-01,
     +     -0.386541758762774673113423570465D-01/
C
C  Weights of first null rule of degree 3
C
      DATA (WEIGHT(4,I),I=1,ORBITS)/-
     +     0.775992773232808462404390159802D-01,
     +     -0.527453289659022924847298408064D-01,
     +     0.145876238555932704488677626554D-01,
     +     0.739374873393616192857532718429D-02,
     +     -0.374618791364332892611678523428D-01,
     +     0.538502846550653076078817013885D-01,
     +     -0.183980865177843057548322735665D-01/
C
C  Weights of second null rule of degree 3
C
      DATA (WEIGHT(5,I),I=1,ORBITS)/0.181767621501470154602720474731D-01
     +     ,0.179938831310058580533178529022D-01,
     +     0.713210362750414891598257378898D-01,
     +     -0.443935688958258805893448212636D-01,
     +     -0.657639036547720234169662790056D-01,
     +     -0.101551807522541414699808460583D-01,
     +     0.265486188970540796821750584204D-01/
C
C  Weights of null rule of degree 2
C
      DATA (WEIGHT(6,I),I=1,ORBITS)/-
     +     0.867629853722843888927184699428D-01,
     +     -0.715881271235661902772072127812D-01,
     +     0.886720767790426261677273459523D-02,
     +     -0.577885573028655167063092577589D-01,
     +     0.430310167581202031805055255554D-01,
     +     -0.606467834856775537069463817445D-02,
     +     0.319492443333738343104163265406D-01/
C
C   Weights of null rule of degree 1
C
      DATA (WEIGHT(7,I),I=1,ORBITS)/0.510374015624925451319499382594D-01
     +     ,0.463998830432033721597269299429D-01,
     +     -0.191086148397852799983451475821D-01,
     +     -0.973768821003670776204287367278D-01,
     +     0.180352562073914141268335496511D-01,
     +     0.277129527093489643801598303110D-01,
     +     -0.176218263109360550515567818653D-01/
C
C  The number of points used by the cubature formula is
C     NUM    = K(0) + 4*K(1) + 6*K(2) + 12*K(3) = 43
C
C***FIRST EXECUTABLE STATEMENT DRLTRI
      TRES = 50*DLAMCH('p')
C
C  Compute the volume of the tetrahedron
C
      VOLUME = (VER(1,2)-VER(1,1))* ((VER(2,3)-VER(2,1))*
     +         (VER(3,4)-VER(3,1))- (VER(2,4)-VER(2,1))*
     +         (VER(3,3)-VER(3,1))) - (VER(2,2)-VER(2,1))*
     +         ((VER(1,3)-VER(1,1))* (VER(3,4)-VER(3,1))-
     +         (VER(1,4)-VER(1,1))* (VER(3,3)-VER(3,1))) +
     +         (VER(3,2)-VER(3,1))* ((VER(1,3)-VER(1,1))*
     +         (VER(2,4)-VER(2,1))- (VER(1,4)-VER(1,1))*
     +         (VER(2,3)-VER(2,1)))
      VOLUME = ABS(VOLUME)/6
C
C  Initialise BASVAL and NULL
C
      DO 10 J = 1,NUMFUN
          BASVAL(J) = 0.d0
          DO 20 I = 1,6
              NULL(J,I) = 0.d0
   20     CONTINUE
   10 CONTINUE
      P = 1
C
C  Compute contributions from orbits with 1, 4, 6 and 12  points
C
      DO 4 TYPE = 0,3
          DO 3 NR = 1,K(TYPE)
              IF (TYPE.EQ.1) THEN
C                 Generator ( z(1) , z(2), z(2) ; z(2) )
                  Z(1) = TYPE1(NR)
                  Z(2) = (1-Z(1))/3
              ELSE IF (TYPE.EQ.2) THEN
C                 Generator ( z(1) , z(1), z(2) ; z(2) )
                  Z(1) = TYPE2(NR)
                  Z(2) = (1-2*Z(1))/2
              ELSE
C                 Generator ( z(1) , z(2), z(3) ; z(3) )
                  Z(1) = TYPE3(1,NR)
                  Z(2) = TYPE3(2,NR)
                  Z(3) = (1-Z(1)-Z(2))/2
              END IF
C                                RGNERR is work array for DORTHE
              CALL DORTET(TYPE,Z,VER,NUMFUN,FUNSUB,SUMVAL,RGNERR)
              if(iero.ne.0) return
              DO 2 J = 1,NUMFUN
                  BASVAL(J) = BASVAL(J) + WEIGHT(1,P)*SUMVAL(J)
                  DO 1 I = 1,6
                      NULL(J,I) = NULL(J,I) + WEIGHT(I+1,P)*SUMVAL(J)
    1             CONTINUE
    2         CONTINUE
              P = P + 1
    3     CONTINUE
    4 CONTINUE
C
C  Compute error estimates
C
      GREATE = 0
      DO 30 J = 1,NUMFUN
          NOISE = ABS(BASVAL(J))*TRES
          DEG4 = SQRT(NULL(J,1)**2+NULL(J,2)**2)
          DEG3 = SQRT(NULL(J,3)**2+NULL(J,4)**2)
          IF (DEG4.LE.NOISE) THEN
              RGNERR(J) = NOISE
          ELSE
              DEG1 = SQRT(NULL(J,5)**2+NULL(J,6)**2)
              IF (DEG3.NE.0) THEN
                  R1 = (DEG4/DEG3)**2
              ELSE
                  R1 = 1
              END IF
              IF (DEG1.NE.0) THEN
                  R2 = DEG3/DEG1
              ELSE
                  R2 = 1
              END IF
              R = MAX(R1,R2)
              IF (R.GE.CRIVAL) THEN
                  RGNERR(J) = FACMED*R*DEG4
              ELSE
                  RGNERR(J) = FACOPT* (R**2)*DEG4
              END IF
              RGNERR(J) = MAX(NOISE,RGNERR(J))
          END IF
          RGNERR(J) = VOLUME*RGNERR(J)
          BASVAL(J) = VOLUME*BASVAL(J)
          GREATE = MAX(GREATE,RGNERR(J))
   30 CONTINUE
C
C***END DRLTET
C
      RETURN
      END

      SUBROUTINE DORTET(TYPE,GENER,VER,NUMFUN,FUNSUB,SUMVAL,WORK)
C***BEGIN PROLOGUE DORTET
C***PURPOSE  To compute the sum of function values over all points
C            of an orbit.
C***DESCRIPTION
C   ON ENTRY
C
C   TYPE   Integer
C          The type of the orbit.
C   GENER  Integer array of dimension (3).
C          The generator for the orbit in homogeneous coordinates.
C   VER    Real array of dimension (3,4).
C          The coordinates of the vertices of the tetrahedron.
C          vertex i -> ( ver(1,i),ver(2,i),ver(3,i) )
C   NUMFUN Integer.
C          Number of components of the vector integrand.
C   FUNSUB Externally declared subroutine for computing
C            all components of the integrand at the given
C            evaluation point.
C            It must have parameters (X,NUMFUN,FUNVLS)
C            Input parameters:
C              X(1)   The x-coordinate of the evaluation point.
C              X(2)   The y-coordinate of the evaluation point.
C              X(3)   The z-coordinate of the evaluation point.
C              NUMFUN Integer that defines the number of
C                     components of the vector integrand.
C            Output parameter:
C              FUNVLS Real array of dimension NUMFUN
C                     that defines NUMFUN components of the integrand.
C   WORK   Real array of dimension (NUMFUN).
C          A work array
C
C   ON RETURN
C
C   SUMVAL Real array of dimension (NUMFUN).
C          The sum of function values over all points
C          of the given orbit.
C
C***END PROLOGUE  DORTET
C
C  Global variables
C
      INTEGER NUMFUN,TYPE
      DOUBLE PRECISION GENER(3),SUMVAL(NUMFUN),VER(3,4),WORK(NUMFUN)
      EXTERNAL FUNSUB
      common/ierdcu/iero
C
C  Local variables
C
      INTEGER I,J,NUMBER
      DOUBLE PRECISION X(3,12),Z1,Z2,Z3
C
C  The array 'gener' contains the necessary information
C   for the generator
C
C***FIRST EXECUTABLE STATEMENT  DORTET
      GO TO (100,110,120,130),TYPE + 1
C
C  Generator with homogeneous coordinates (1/4,1/4,1/4;1/4)
C
  100 NUMBER = 1
      X(1,1) = (VER(1,1)+VER(1,2)+VER(1,3)+VER(1,4))/4
      X(2,1) = (VER(2,1)+VER(2,2)+VER(2,3)+VER(2,4))/4
      X(3,1) = (VER(3,1)+VER(3,2)+VER(3,3)+VER(3,4))/4
      GO TO 200
C
C  Generator with homogeneous coordinates (z1,z2,z2;z2)
C
  110 NUMBER = 4
      Z1 = GENER(1)
      Z2 = GENER(2)
      DO 115 J = 1,3
          X(J,1) = Z1*VER(J,1) + Z2* (VER(J,2)+VER(J,3)+VER(J,4))
          X(J,2) = Z1*VER(J,2) + Z2* (VER(J,1)+VER(J,3)+VER(J,4))
          X(J,3) = Z1*VER(J,3) + Z2* (VER(J,2)+VER(J,1)+VER(J,4))
          X(J,4) = Z1*VER(J,4) + Z2* (VER(J,2)+VER(J,3)+VER(J,1))
  115 CONTINUE
      GO TO 200
C
C  Generator with homogeneous coordinates (z1,z1,z2;z2)
C
  120 NUMBER = 6
      Z1 = GENER(1)
      Z2 = GENER(2)
      DO 125 J = 1,3
          X(J,1) = Z1* (VER(J,1)+VER(J,2)) + Z2* (VER(J,3)+VER(J,4))
          X(J,2) = Z1* (VER(J,1)+VER(J,3)) + Z2* (VER(J,2)+VER(J,4))
          X(J,3) = Z1* (VER(J,1)+VER(J,4)) + Z2* (VER(J,3)+VER(J,2))
          X(J,4) = Z1* (VER(J,2)+VER(J,3)) + Z2* (VER(J,1)+VER(J,4))
          X(J,5) = Z1* (VER(J,2)+VER(J,4)) + Z2* (VER(J,1)+VER(J,3))
          X(J,6) = Z1* (VER(J,3)+VER(J,4)) + Z2* (VER(J,1)+VER(J,2))
  125 CONTINUE
      GO TO 200
C
C  Generator with homogeneous coordinates (z1,z2,z3;z3)
C
  130 NUMBER = 12
      Z1 = GENER(1)
      Z2 = GENER(2)
      Z3 = GENER(3)
      DO 135 J = 1,3
          X(J,1) = Z1*VER(J,1) + Z2*VER(J,2) + Z3* (VER(J,3)+VER(J,4))
          X(J,2) = Z1*VER(J,1) + Z2*VER(J,3) + Z3* (VER(J,2)+VER(J,4))
          X(J,3) = Z1*VER(J,1) + Z2*VER(J,4) + Z3* (VER(J,2)+VER(J,3))
          X(J,4) = Z1*VER(J,2) + Z2*VER(J,1) + Z3* (VER(J,3)+VER(J,4))
          X(J,5) = Z1*VER(J,2) + Z2*VER(J,3) + Z3* (VER(J,1)+VER(J,4))
          X(J,6) = Z1*VER(J,2) + Z2*VER(J,4) + Z3* (VER(J,1)+VER(J,3))
          X(J,7) = Z1*VER(J,3) + Z2*VER(J,1) + Z3* (VER(J,2)+VER(J,4))
          X(J,8) = Z1*VER(J,3) + Z2*VER(J,2) + Z3* (VER(J,1)+VER(J,4))
          X(J,9) = Z1*VER(J,3) + Z2*VER(J,4) + Z3* (VER(J,1)+VER(J,2))
          X(J,10) = Z1*VER(J,4) + Z2*VER(J,1) + Z3* (VER(J,2)+VER(J,3))
          X(J,11) = Z1*VER(J,4) + Z2*VER(J,2) + Z3* (VER(J,1)+VER(J,3))
          X(J,12) = Z1*VER(J,4) + Z2*VER(J,3) + Z3* (VER(J,1)+VER(J,2))
  135 CONTINUE
  200 CONTINUE
      CALL FUNSUB(X(1,1),NUMFUN,SUMVAL)
      if(iero.ne.0) return
      DO 220 J = 2,NUMBER
          CALL FUNSUB(X(1,J),NUMFUN,WORK)
          if(iero.ne.0) return
          DO 210 I = 1,NUMFUN
              SUMVAL(I) = SUMVAL(I) + WORK(I)
  210     CONTINUE
  220 CONTINUE
C
C***END DORTET
C
      RETURN
      END