MODULE adapt_quad
! Module for adaptive quadrature.
! Adapted from TOMS algorithm 691.
! This version uses the ELF90 subset of Fortran 90.
! Conversion by Alan Miller
! amiller @ bigpond.net.au

CONTAINS

SUBROUTINE qpsrt(limit, last, maxerr, ermax, elist, iord, nrmax)
!     ..................................................................

! 1.   QPSRT
!      ORDERING ROUTINE
!         STANDARD FORTRAN SUBROUTINE
!         REAL VERSION

! 2.   PURPOSE
!         THIS ROUTINE MAINTAINS THE DESCENDING ORDERING IN THE LIST OF THE
!         LOCAL ERROR ESTIMATES RESULTING FROM THE INTERVAL SUBDIVISION
!         PROCESS.  AT EACH CALL TWO ERROR ESTIMATES ARE INSERTED USING THE
!         SEQUENTIAL SEARCH METHOD, TOP-DOWN FOR THE LARGEST ERROR ESTIMATE
!         AND BOTTOM-UP FOR THE SMALLEST ERROR ESTIMATE.

! 3.   CALLING SEQUENCE
!         CALL QPSRT(LIMIT, LAST, MAXERR, ERMAX, ELIST, IORD, NRMAX)

!      PARAMETERS (MEANING AT OUTPUT)
!         LIMIT  - INTEGER
!                  MAXIMUM NUMBER OF ERROR ESTIMATES THE LIST CAN CONTAIN

!         LAST   - INTEGER
!                  NUMBER OF ERROR ESTIMATES CURRENTLY IN THE LIST

!         MAXERR - INTEGER
!                  MAXERR POINTS TO THE NRMAX-TH LARGEST ERROR ESTIMATE
!                  CURRENTLY IN THE LIST

!         ERMAX  - REAL
!                  NRMAX-TH LARGEST ERROR ESTIMATE
!                  ERMAX = ELIST(MAXERR)

!         ELIST  - REAL
!                  VECTOR OF DIMENSION LAST CONTAINING THE ERROR ESTIMATES

!         IORD   - INTEGER
!                  VECTOR OF DIMENSION LAST, THE FIRST K ELEMENTS OF
!                  WHICH CONTAIN POINTERS TO THE ERROR ESTIMATES,
!                  SUCH THAT ELIST(IORD(1)), ... , ELIST(IORD(K))
!                  FORM A DECREASING SEQUENCE, WITH K = LAST IF
!                  LAST <= (LIMIT/2+2), AND K = LIMIT+1-LAST OTHERWISE

!         NRMAX  - INTEGER
!                  MAXERR = IORD(NRMAX)

! 4.   NO SUBROUTINES OR FUNCTIONS NEEDED

!     ..................................................................

USE constants_NSWC
IMPLICIT NONE

INTEGER, INTENT(IN)                 :: limit, last
REAL (dp), DIMENSION(:), INTENT(IN) :: elist
INTEGER, INTENT(IN OUT)             :: nrmax
INTEGER, DIMENSION(:), INTENT(OUT)  :: iord
INTEGER, INTENT(OUT)                :: maxerr
REAL (dp), INTENT(OUT)              :: ermax

REAL (dp) :: errmax, errmin
INTEGER   :: i, ibeg, ido, isucc, j, jbnd, jupbn, k

!           CHECK WHETHER THE LIST CONTAINS MORE THAN TWO ERROR ESTIMATES.

!***FIRST EXECUTABLE STATEMENT  QPSRT
IF(last > 2) GO TO 10
iord(1) = 1
iord(2) = 2
GO TO 90

!           THIS PART OF THE ROUTINE IS ONLY EXECUTED IF,
!           DUE TO A DIFFICULT INTEGRAND, SUBDIVISION INCREASED
!           THE ERROR ESTIMATE.   IN THE NORMAL CASE THE INSERT PROCEDURE
!           SHOULD START AFTER THE NRMAX-TH LARGEST ERROR ESTIMATE.

10 errmax = elist(maxerr)
IF(nrmax == 1) GO TO 30
ido = nrmax-1
DO i = 1, ido
  isucc = iord(nrmax-1)
! ***JUMP OUT OF DO-LOOP
  IF(errmax <= elist(isucc)) EXIT
  iord(nrmax) = isucc
  nrmax = nrmax-1
END DO

!           COMPUTE THE NUMBER OF ELEMENTS IN THE LIST TO BE MAINTAINED
!           IN DESCENDING ORDER.  THIS NUMBER DEPENDS ON THE NUMBER OF
!           SUBDIVISIONS STILL ALLOWED.

30 jupbn = last
IF(last > (limit/2+2)) jupbn = limit+3-last
errmin = elist(last)

!           INSERT ERRMAX BY TRAVERSING THE LIST TOP-DOWN,
!           STARTING COMPARISON FROM THE ELEMENT ELIST(IORD(NRMAX+1)).

jbnd = jupbn-1
ibeg = nrmax+1
DO i=ibeg, jbnd
  isucc = iord(i)
! ***JUMP OUT OF DO-LOOP
  IF(errmax >= elist(isucc)) GO TO 60
  iord(i-1) = isucc
END DO
iord(jbnd) = maxerr
iord(jupbn) = last
GO TO 90

!           INSERT ERRMIN BY TRAVERSING THE LIST BOTTOM-UP.

60 iord(i-1) = maxerr
k = jbnd
DO j=i, jbnd
  isucc = iord(k)
! ***JUMP OUT OF DO-LOOP
  IF(errmin < elist(isucc)) GO TO 80
  iord(k+1) = isucc
  k = k-1
END DO
iord(i) = last
GO TO 90
80 iord(k+1) = last

!           SET MAXERR AND ERMAX.

90 maxerr = iord(nrmax)
ermax = elist(maxerr)
RETURN
END SUBROUTINE qpsrt


SUBROUTINE qelg (n, epstab, result, abserr, res3la, nres, epmach, oflow)
!-----------------------------------------------------------------------

! 1.   PURPOSE
!         THE ROUTINE DETERMINES THE LIMIT OF A GIVEN SEQUENCE OF
!         APPROXIMATIONS, BY MEANS OF THE EPSILON ALGORITHM OF P. WYNN.
!         AN ESTIMATE OF THE ABSOLUTE ERROR IS ALSO GIVEN.
!         THE CONDENSED EPSILON TABLE IS COMPUTED. ONLY THOSE ELEMENTS
!         NEEDED FOR THE COMPUTATION OF THE NEXT DIAGONAL ARE PRESERVED.

! 2.   PARAMETERS
!         N      - INTEGER
!                  EPSTAB(N) CONTAINS THE NEW ELEMENT IN THE
!                  FIRST COLUMN OF THE EPSILON TABLE.

!         EPSTAB - REAL
!                  VECTOR OF DIMENSION 52 CONTAINING THE ELEMENTS OF THE
!                  TWO LOWER DIAGONALS OF THE TRIANGULAR EPSILON TABLE.
!                  THE ELEMENTS ARE NUMBERED STARTING AT THE RIGHT-HAND
!                  CORNER OF THE TRIANGLE.

!         RESULT - REAL
!                  RESULTING APPROXIMATION TO THE INTEGRAL

!         ABSERR - REAL
!                  ESTIMATE OF THE ABSOLUTE ERROR COMPUTED FROM
!                  RESULT AND THE 3 PREVIOUS RESULTS

!         RES3LA - REAL
!                  VECTOR OF DIMENSION 3 CONTAINING THE LAST 3 RESULTS

!         NRES   - INTEGER
!                  NUMBER OF CALLS TO THE ROUTINE
!                  (SHOULD BE ZERO AT FIRST CALL)

!         EPMACH - REAL
!                  THE RELATIVE PRECISION OF THE FLOATING ARITHMETIC
!                  BEING USED.

!         OFLOW  - REAL
!                  THE LARGEST POSITIVE MAGNITUDE.

! 3.   NO SUBROUTINES OR FUNCTIONS USED

!-----------------------------------------------------------------------
USE constants_NSWC
IMPLICIT NONE

INTEGER, INTENT(IN OUT)                 :: n, nres
REAL (dp), INTENT(IN)                   :: epmach, oflow
REAL (dp), INTENT(OUT)                  :: abserr, result
REAL (dp), DIMENSION(:), INTENT(IN OUT) :: epstab, res3la
!---------------------

!      LIST OF MAJOR VARIABLES
!      -----------------------

!      E0     - THE 4 ELEMENTS ON WHICH THE
!      E1       COMPUTATION OF A NEW ELEMENT IN
!      E2       THE EPSILON TABLE IS BASED
!      E3                 E0
!                   E3    E1    NEW
!                         E2
!      NEWELM - NUMBER OF ELEMENTS TO BE COMPUTED IN THE NEW DIAGONAL
!      ERROR  - ERROR = ABS(E1-E0)+ABS(E2-E1)+ABS(NEW-E2)
!      RESULT - THE ELEMENT IN THE NEW DIAGONAL WITH LEAST VALUE OF ERROR

!      LIMEXP IS THE MAXIMUM NUMBER OF ELEMENTS THE EPSILON
!      TABLE CAN CONTAIN. IF THIS NUMBER IS REACHED, THE UPPER
!      DIAGONAL OF THE EPSILON TABLE IS DELETED.

REAL (dp) :: delta1, delta2, delta3, epsinf, error, err1, err2, err3, e0, &
             e1, e1abs, e2, e3, res, ss, tol1, tol2, tol3
INTEGER   :: i, ib, ib2, ie, indx, k1, k2, k3, limexp, newelm, num

nres = nres + 1
abserr = oflow
result = epstab(n)
IF (n < 3) GO TO 100
limexp = 50
epstab(n + 2) = epstab(n)
newelm = (n - 1)/2
epstab(n) = oflow
num = n
k1 = n
DO i = 1, newelm
  k2 = k1 - 1
  k3 = k1 - 2
  res = epstab(k1 + 2)
  e0 = epstab(k3)
  e1 = epstab(k2)
  e2 = res
  e1abs = ABS(e1)
  delta2 = e2 - e1
  err2 = ABS(delta2)
  tol2 = MAX(ABS(e2),e1abs)*epmach
  delta3 = e1 - e0
  err3 = ABS(delta3)
  tol3 = MAX(e1abs,ABS(e0))*epmach
  IF (err2 > tol2 .OR. err3 > tol3) GO TO 10
  
!      IF E0, E1 AND E2 ARE EQUAL TO WITHIN MACHINE
!      ACCURACY, CONVERGENCE IS ASSUMED.
!      RESULT = E2
!      ABSERR = ABS(E1-E0) + ABS(E2-E1)
  
  result = res
  abserr = err2 + err3
! ***JUMP OUT OF DO-LOOP
  GO TO 100
  10 e3 = epstab(k1)
  epstab(k1) = e1
  delta1 = e1 - e3
  err1 = ABS(delta1)
  tol1 = MAX(e1abs,ABS(e3))*epmach
  
!      IF TWO ELEMENTS ARE VERY CLOSE TO EACH OTHER, OMIT
!      A PART OF THE TABLE BY ADJUSTING THE VALUE OF N
  
  IF (err1 <= tol1 .OR. err2 <= tol2 .OR. err3 <= tol3) GO TO 20
  ss = 1.0D0/delta1 + 1.0D0/delta2 - 1.0D0/delta3
  epsinf = ABS(ss*e1)
  
!      TEST TO DETECT IRREGULAR BEHAVIOUR IN THE TABLE, AND EVENTUALLY
!      OMIT A PART OF THE TABLE ADJUSTING THE VALUE OF N.
  
  IF (epsinf > 0.1D-03) GO TO 30
  20 n = i + i - 1
! ***JUMP OUT OF DO-LOOP
  GO TO 50
  
!      COMPUTE A NEW ELEMENT AND EVENTUALLY ADJUST THE VALUE OF RESULT.
  
  30 res = e1 + 1.0D0/ss
  epstab(k1) = res
  k1 = k1 - 2
  error = err2 + ABS(res - e2) + err3
  IF (error > abserr) CYCLE
  abserr = error
  result = res
END DO

!      SHIFT THE TABLE.

50 IF (n == limexp) n = 2*(limexp/2) - 1
ib = 1
IF ((num/2)*2 == num) ib = 2
ie = newelm + 1
DO i = 1, ie
  ib2 = ib + 2
  epstab(ib) = epstab(ib2)
  ib = ib2
END DO
IF (num == n) GO TO 80
indx = num - n + 1
DO i = 1, n
  epstab(i) = epstab(indx)
  indx = indx + 1
END DO
80 IF (nres >= 4) GO TO 90
res3la(nres) = result
abserr = oflow
GO TO 100

!      COMPUTE ERROR ESTIMATE

90 abserr = ABS(result - res3la(3)) + ABS(result - res3la(2)) +  &
            ABS(result - res3la(1))
res3la(1) = res3la(2)
res3la(2) = res3la(3)
res3la(3) = result
100 abserr = MAX(abserr, 5.0D0*epmach*ABS(result))
RETURN
END SUBROUTINE qelg


SUBROUTINE qxgs (f, a, b, epsabs, epsrel, result, abserr, ier, limit, last)

!     THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
!     DEFINITE INTEGRAL  I = INTEGRAL OF F OVER (A,B),
!     HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
!     ABS(I-RESULT) <= MAX(EPSABS, EPSREL*ABS(I)).

! PARAMETERS
!  ON ENTRY
!     F      - REAL
!              FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
!              FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
!              DECLARED E X T E R N A L IN THE DRIVER PROGRAM.

!     A      - REAL
!              LOWER LIMIT OF INTEGRATION

!     B      - REAL
!              UPPER LIMIT OF INTEGRATION

!     EPSABS - REAL
!              ABSOLUTE ACCURACY REQUESTED

!     EPSREL - REAL
!              RELATIVE ACCURACY REQUESTED

!  ON RETURN
!     RESULT - REAL
!              APPROXIMATION TO THE INTEGRAL

!     ABSERR - REAL
!              ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
!              WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)

!     IER    - INTEGER
!              IER = 0 NORMAL AND RELIABLE TERMINATION OF THE ROUTINE.
!                      IT IS ASSUMED THAT THE REQUESTED ACCURACY HAS
!                      BEEN ACHIEVED.
!              IER > 0 ABNORMAL TERMINATION OF THE ROUTINE
!                      THE ESTIMATES FOR INTEGRAL AND ERROR ARE
!                      LESS RELIABLE. IT IS ASSUMED THAT THE
!                      REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
!     ERROR MESSAGES
!              IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED HAS BEEN
!                      ACHIEVED.  ONE CAN ALLOW MORE SUB-DIVISIONS BY
!                      INCREASING THE VALUE OF LIMIT (AND TAKING THE ACCORDING
!                      DIMENSION ADJUSTMENTS INTO ACCOUNT.  HOWEVER, IF THIS
!                      YIELDS NO IMPROVEMENT IT IS ADVISED TO ANALYZE THE
!                      INTEGRAND IN ORDER TO DETERMINE THE INTEGRATION
!                      DIFFICULTIES.  IF THE POSITION OF A LOCAL DIFFICULTY
!                      CAN BE DETERMINED (E.G. SINGULARITY, DISCONTINUITY
!                      WITHIN THE INTERVAL) ONE WILL PROBABLY GAIN FROM
!                      SPLITTING UP THE INTERVAL AT THIS POINT AND CALLING THE
!                      INTEGRATOR ON THE SUBRANGES.  IF POSSIBLE, AN
!                      APPROPRIATE SPECIAL-PURPOSE INTEGRATOR SHOULD BE USED,
!                      WHICH IS DESIGNED FOR HANDLING THE TYPE OF DIFFICULTY
!                      INVOLVED.
!                  = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS DETECTED,
!                      WHICH PREVENTS THE REQUESTED TOLERANCE FROM BEING
!                      ACHIEVED.
!                      THE ERROR MAY BE UNDER-ESTIMATED.
!                  = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR
!                      OCCURS AT SOME POINTS OF THE INTEGRATION INTERVAL.
!                  = 4 THE ALGORITHM DOES NOT CONVERGE.
!                      ROUNDOFF ERROR IS DETECTED IN THE EXTRAPOLATION TABLE.
!                      IT IS PRESUMED THAT THE REQUESTED TOLERANCE CANNOT BE
!                      ACHIEVED, AND THAT THE RETURNED RESULT IS THE BEST
!                      WHICH CAN BE OBTAINED.
!                  = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR SLOWLY CONVERGENT.
!                      IT MUST BE NOTED THAT DIVERGENCE CAN OCCUR WITH ANY
!                      OTHER VALUE OF IER.
!                  = 6 THE INPUT IS INVALID BECAUSE EPSABS OR EPSREL IS
!                      NEGATIVE, LIMIT < 1, LENW < 46*LIMIT, OR LENIW < 3*LIMIT.
!                      RESULT, ABSERR, LAST ARE SET TO ZERO.
!                      EXCEPT WHEN LIMIT OR LENW OR LENIW IS INVALID, IWORK(1),
!                      WORK(LIMIT*2+1) AND WORK(LIMIT*3+1) ARE SET TO ZERO,
!                      WORK(1) IS SET TO A, AND WORK(LIMIT+1) TO B.

!  DIMENSIONING PARAMETERS
!     LIMIT - INTEGER
!             LIMIT DETERMINES THE MAXIMUM NUMBER OF SUBINTERVALS IN THE
!             PARTITION OF THE GIVEN INTEGRATION INTERVAL (A,B), LIMIT >= 1.
!             IF LIMIT < 1, THE ROUTINE WILL END WITH IER = 6.

!     LAST  - INTEGER
!             ON RETURN, LAST EQUALS THE NUMBER OF SUBINTERVALS PRODUCED
!             IN THE SUBDIVISION PROCESS, WHICH DETERMINES THE NUMBER OF
!             SIGNIFICANT ELEMENTS ACTUALLY IN THE WORK ARRAYS.

USE constants_NSWC
IMPLICIT NONE

REAL (dp), INTENT(IN)  :: a, b, epsabs, epsrel
REAL (dp), INTENT(OUT) :: result, abserr
INTEGER, INTENT(IN)    :: limit
INTEGER, INTENT(OUT)   :: ier, last

INTERFACE
  FUNCTION f(x) RESULT(fx)
    USE constants_NSWC
    IMPLICIT NONE
    REAL (dp), INTENT(IN) :: x
    REAL (dp)             :: fx
  END FUNCTION f
END INTERFACE

!         CHECK VALIDITY OF LIMIT

ier = 6
last = 0
result = 0.0D0
abserr = 0.0D0
IF (limit < 1) RETURN

!         PREPARE CALL FOR QXGSE.

CALL qxgse(f, a, b, epsabs, epsrel, limit, result, abserr, ier, last)

RETURN
END SUBROUTINE qxgs


SUBROUTINE qxgse(f, a, b, epsabs, epsrel, limit, result, abserr, ier, last)

!       THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A
!       DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B),
!       HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
!       ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).

!   PARAMETERS
!    ON ENTRY
!       F      - REAL
!                FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
!                FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
!                DECLARED E X T E R N A L IN THE DRIVER PROGRAM.

!       A      - REAL
!                LOWER LIMIT OF INTEGRATION

!       B      - REAL
!                UPPER LIMIT OF INTEGRATION

!       EPSABS - REAL
!                ABSOLUTE ACCURACY REQUESTED

!       EPSREL - REAL
!                RELATIVE ACCURACY REQUESTED

!       LIMIT  - INTEGER
!                GIVES AN UPPERBOUND ON THE NUMBER OF SUBINTERVALS
!                IN THE PARTITION OF (A,B)

!    ON RETURN
!       RESULT - REAL
!                APPROXIMATION TO THE INTEGRAL

!       ABSERR - REAL
!                ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
!                WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)

!       IER    - INTEGER
!                IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
!                        ROUTINE. IT IS ASSUMED THAT THE REQUESTED
!                        ACCURACY HAS BEEN ACHIEVED.
!                IER > 0 ABNORMAL TERMINATION OF THE ROUTINE
!                        THE ESTIMATES FOR INTEGRAL AND ERROR ARE
!                        LESS RELIABLE. IT IS ASSUMED THAT THE
!                        REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
!       ERROR MESSAGES
!                    = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
!                        HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB-
!                        DIVISIONS BY INCREASING THE VALUE OF LIMIT
!                        (AND TAKING THE ACCORDING DIMENSION
!                        ADJUSTMENTS INTO ACCOUNT). HOWEVER, IF
!                        THIS YIELDS NO IMPROVEMENT IT IS ADVISED
!                        TO ANALYZE THE INTEGRAND IN ORDER TO
!                        DETERMINE THE INTEGRATION DIFFICULTIES. IF
!                        THE POSITION OF A LOCAL DIFFICULTY CAN BE
!                        DETERMINED (E.G. SINGULARITY,
!                        DISCONTINUITY WITHIN THE INTERVAL) ONE
!                        WILL PROBABLY GAIN FROM SPLITTING UP THE
!                        INTERVAL AT THIS POINT AND CALLING THE
!                        INTEGRATOR ON THE SUBRANGES. IF POSSIBLE,
!                        AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
!                        SHOULD BE USED, WHICH IS DESIGNED FOR
!                        HANDLING THE TYPE OF DIFFICULTY INVOLVED.
!                    = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS DETEC-
!                        TED, WHICH PREVENTS THE REQUESTED
!                        TOLERANCE FROM BEING ACHIEVED.
!                        THE ERROR MAY BE UNDER-ESTIMATED.
!                    = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS AT
!                        SOME POINTS OF THE INTEGRATION INTERVAL.
!                    = 4 THE ALGORITHM DOES NOT CONVERGE.
!                        ROUNDOFF ERROR IS DETECTED IN THE
!                        EXTRAPOLATION TABLE.
!                        IT IS PRESUMED THAT THE REQUESTED TOLERANCE
!                        CANNOT BE ACHIEVED, AND THAT THE RETURNED RESULT
!                        IS THE BEST WHICH CAN BE OBTAINED.
!                    = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR SLOWLY
!                        CONVERGENT.   IT MUST BE NOTED THAT DIVERGENCE
!                        CAN OCCUR WITH ANY OTHER VALUE OF IER.
!                    = 6 THE INPUT IS INVALID BECAUSE EPSABS OR
!                        EPSREL IS NEGATIVE. RESULT, ABSERR,
!                        LAST, RLIST(1), IORD(1), AND ELIST(1)
!                        ARE SET TO ZERO. ALIST(1) AND BLIST(1)
!                        ARE SET TO A AND B RESPECTIVELY.

!       ALIST  - REAL
!                VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
!                 LAST  ELEMENTS OF WHICH ARE THE LEFT END POINTS
!                OF THE SUBINTERVALS IN THE PARTITION OF THE
!                GIVEN INTEGRATION RANGE (A,B)

!       BLIST  - REAL
!                VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
!                 LAST  ELEMENTS OF WHICH ARE THE RIGHT END POINTS
!                OF THE SUBINTERVALS IN THE PARTITION OF THE GIVEN
!                INTEGRATION RANGE (A,B)

!       RLIST  - REAL
!                VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST `LAST'
!                ELEMENTS OF WHICH ARE THE INTEGRAL APPROXIMATIONS ON
!                THE SUBINTERVALS

!       ELIST  - REAL
!                VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
!                 LAST  ELEMENTS OF WHICH ARE THE MODULI OF THE
!                ABSOLUTE ERROR ESTIMATES ON THE SUBINTERVALS

!       IORD   - INTEGER
!                VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST K ELEMENTS
!                OF WHICH ARE POINTERS TO THE ERROR ESTIMATES OVER THE
!                SUBINTERVALS, SUCH THAT ELIST(IORD(1)), ..., ELIST(IORD(K))
!                FORM A DECREASING SEQUENCE, WITH K = LAST IF
!                LAST <= (LIMIT/2+2), AND K = LIMIT+1-LAST OTHERWISE

!       LAST   - INTEGER
!                NUMBER OF SUBINTERVALS ACTUALLY PRODUCED IN THE
!                SUBDIVISION PROCESS

!       VALP   - REAL
!       VALN     ARRAYS OF DIMENSION AT LEAST (21,LIMIT) USED TO
!                SAVE THE FUNCTIONAL VALUES

!       LP     - INTEGER
!       LN       VECTORS OF DIMENSION AT LEAST LIMIT, USED TO
!                STORE THE ACTUAL NUMBER OF FUNCTIONAL VALUES
!                SAVED IN THE CORRESPONDING COLUMN OF VALP,VALN

!***ROUTINES CALLED  F, SPMPAR, QELG, QXLQM, QPSRT, QXRRD, QXCPY

USE constants_NSWC
IMPLICIT NONE

REAL (dp), INTENT(IN)   :: a, b, epsabs, epsrel
REAL (dp), INTENT(OUT)  :: result, abserr
INTEGER, INTENT(IN)     :: limit
INTEGER, INTENT(OUT)    :: ier, last

INTERFACE
  FUNCTION f(x) RESULT(fx)
    USE constants_NSWC
    IMPLICIT NONE
    REAL (dp), INTENT(IN) :: x
    REAL (dp)             :: fx
  END FUNCTION f
END INTERFACE

REAL (dp) :: abseps, alist(limit), area, area1, area12, area2, a1, a2, b1, &
             b2, blist(limit), correc, defab1, defab2, dres, elist(limit), &
             epmach, erlarg, erlast, errbnd, errmax, error1, error2,       &
             erro12, errsum, ertest, oflow, rerr, resabs, reseps, res3la(3), &
             rlist(limit), rlist2(52), small, t, uflow, valp(21,limit),    &
             valn(21,limit), vp1(21), vp2(21), vn1(21), vn2(21), defabs
INTEGER   :: id, ierro, iord(limit), iroff1, iroff2, iroff3, jupbnd, k,  &
             ksgn, lp(limit), ln(limit), ktmin, maxerr, nres, nrmax, numrl2, &
             lp1, lp2, ln1, ln2
LOGICAL   :: extrap, noext

!            MACHINE DEPENDENT CONSTANTS
!            ---------------------------

!          EPMACH IS THE LARGEST RELATIVE SPACING.
!          UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
!          OFLOW IS THE LARGEST POSITIVE MAGNITUDE.

epmach = dpmpar(1)
uflow = dpmpar(2)
oflow = dpmpar(3)

!            TEST ON VALIDITY OF PARAMETERS
!            ------------------------------
last = 0
result = 0.0D0
abserr = 0.0D0
alist(1) = a
blist(1) = b
rlist(1) = 0.0D0
elist(1) = 0.0D0
ier = 6
IF (epsabs < 0.0D0 .OR. epsrel < 0.0D0) GO TO 999
ier = 0
rerr = MAX(epsrel, 50.0D0*epmach)

!      FIRST APPROXIMATION TO THE INTEGRAL
!      -----------------------------------

ierro = 0
lp(1) = 1
ln(1) = 1
valp(1,1) = f((a + b)*0.5D0)
valn(1,1) = valp(1,1)
CALL qxlqm(f, a, b, result, abserr, resabs, defabs, valp(1:,1), valn(1:,1), &
           lp(1), ln(1), 2, epmach, uflow, oflow)

!      TEST ON ACCURACY.

dres = ABS(result)
errbnd = MAX(epsabs,rerr*dres)
last = 1
rlist(1) = result
elist(1) = abserr
iord(1) = 1
IF (abserr <= 100.0D0*epmach*defabs .AND. abserr > errbnd) ier = 2
IF (limit == 1) ier = 1
IF (ier /= 0 .OR. (abserr <= errbnd .AND. abserr /= resabs) .OR.  &
abserr == 0.0D0) GO TO 999

!      INITIALIZATION
!      --------------

rlist2(1) = result
errmax = abserr
maxerr = 1
area = result
errsum = abserr
abserr = oflow
nrmax = 1
nres = 0
numrl2 = 2
ktmin = 0
extrap = .false.
noext = .false.
iroff1 = 0
iroff2 = 0
iroff3 = 0
ksgn = -1
IF (dres >= (1.0D0 - 50.0D0*epmach)*defabs) ksgn = 1
t = 1.0D0 + 100.0D0*epmach

!      MAIN DO-LOOP
!      ------------

DO last = 2, limit
  
!      BISECT THE SUBINTERVAL WITH THE NRMAX-TH LARGEST ERROR ESTIMATE.
  
  a1 = alist(maxerr)
  b1 = 0.5D0*(alist(maxerr) + blist(maxerr))
  a2 = b1
  b2 = blist(maxerr)
  erlast = errmax
  CALL qxrrd(f, valn(1:,maxerr), ln(maxerr), b1, a1, vn1, vp1, ln1, lp1)
  CALL qxlqm(f, a1, b1, area1, error1, resabs, defab1, vp1, vn1, lp1, ln1,  &
             2, epmach, uflow, oflow)
  CALL qxrrd(f, valp(1:,maxerr), lp(maxerr), a2, b2, vp2, vn2, lp2, ln2)
  CALL qxlqm(f, a2, b2, area2, error2, resabs, defab2, vp2, vn2, lp2, ln2,  &
             2, epmach, uflow, oflow)
  
!      IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
!      AND ERROR AND TEST FOR ACCURACY.
  
  area12 = area1 + area2
  erro12 = error1 + error2
  errsum = errsum + erro12 - errmax
  area = area + area12 - rlist(maxerr)
  IF (defab1 == error1 .OR. defab2 == error2) GO TO 15
  IF (ABS(rlist(maxerr) - area12) > 0.1D-04*ABS(area12)  &
        .OR. erro12 < 0.99D0*errmax) GO TO 10
  IF (extrap) iroff2 = iroff2 + 1
  IF (.NOT.extrap) iroff1 = iroff1 + 1
  10 IF (last > 10 .AND. erro12 > errmax) iroff3 = iroff3 + 1
  15 rlist(maxerr) = area1
  rlist(last) = area2
  errbnd = MAX(epsabs,rerr*ABS(area))
  
!      TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG.
  
  IF (iroff1 + iroff2 >= 10 .OR. iroff3 >= 20) ier = 2
  IF (iroff2 >= 5) ierro = 3
  
!      SET ERROR FLAG IN THE CASE THAT THE NUMBER OF SUBINTERVALS EQUALS LIMIT.
  
  IF (last == limit) ier = 1
  
!      SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
!      AT A POINT OF THE INTEGRATION RANGE.
  
  IF (MAX(ABS(a1),ABS(b2)) <= t*(ABS(a2) + 1.d+03*uflow)) ier = 4
  
!      APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
  
  IF (error2 > error1) GO TO 20
  alist(last) = a2
  blist(maxerr) = b1
  blist(last) = b2
  elist(maxerr) = error1
  elist(last) = error2
  CALL qxcpy(valp(1:,maxerr), vp1, lp1)
  lp(maxerr) = lp1
  CALL qxcpy(valn(1:,maxerr), vn1, ln1)
  ln(maxerr) = ln1
  CALL qxcpy(valp(1:,last), vp2, lp2)
  lp(last) = lp2
  CALL qxcpy(valn(1:,last), vn2, ln2)
  ln(last) = ln2
  GO TO 30
  
  20 alist(maxerr) = a2
  alist(last) = a1
  blist(last) = b1
  rlist(maxerr) = area2
  rlist(last) = area1
  elist(maxerr) = error2
  elist(last) = error1
  CALL qxcpy(valp(1:,maxerr), vp2, lp2)
  lp(maxerr) = lp2
  CALL qxcpy(valn(1:,maxerr), vn2, ln2)
  ln(maxerr) = ln2
  CALL qxcpy(valp(1:,last), vp1, lp1)
  lp(last) = lp1
  CALL qxcpy(valn(1:,last), vn1, ln1)
  ln(last) = ln1
  
!      CALL SUBROUTINE QPSRT TO MAINTAIN THE DESCENDING ORDERING
!      IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
!      WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT).
  
  30 CALL qpsrt(limit, last, maxerr, errmax, elist, iord, nrmax)
! ***JUMP OUT OF DO-LOOP
  IF(errsum <= errbnd) GO TO 115
! ***JUMP OUT OF DO-LOOP
  IF (ier /= 0) GO TO 100
  IF (last == 2) GO TO 80
  IF (noext) CYCLE
  erlarg = erlarg - erlast
  IF (ABS(b1 - a1) > small) erlarg = erlarg + erro12
  IF (extrap) GO TO 40
  
!      TEST WHETHER THE INTERVAL TO BE BISECTED NEXT IS THE SMALLEST INTERVAL.
  
  IF (ABS(blist(maxerr) - alist(maxerr)) > small) CYCLE
  extrap = .true.
  nrmax = 2
  
!      THE BOUND 0.3*ERTEST HAS BEEN INTRODUCED TO PERFORM A
!      MORE CAUTIOUS EXTRAPOLATION THAN IN THE ORIGINAL DQAGSE ROUTINE
  
  40 IF (ierro == 3 .OR. erlarg <= 0.3D0*ertest) GO TO 60
  
!      THE SMALLEST INTERVAL HAS THE LARGEST ERROR.
!      BEFORE BISECTING DECREASE THE SUM OF THE ERRORS OVER THE
!      LARGER INTERVALS (ERLARG) AND PERFORM EXTRAPOLATION.
  
  id = nrmax
  jupbnd = last
  IF (last > (2 + limit/2)) jupbnd = limit + 3 - last
  DO k = id, jupbnd
    maxerr = iord(nrmax)
    errmax = elist(maxerr)
! ***JUMP OUT OF DO-LOOP
    IF(ABS(blist(maxerr) - alist(maxerr)) > small) CYCLE
    nrmax = nrmax + 1
  END DO
  
!      PERFORM EXTRAPOLATION.
  
  60 numrl2 = numrl2 + 1
  rlist2(numrl2) = area
  CALL qelg (numrl2, rlist2, reseps, abseps, res3la, nres, epmach, oflow)
  ktmin = ktmin + 1
  IF (ktmin > 5 .AND. abserr < 0.1D-02*errsum) ier = 5
  IF (abseps >= abserr) GO TO 70
  ktmin = 0
  abserr = abseps
  result = reseps
  correc = erlarg
  ertest = MAX(epsabs,rerr*ABS(reseps))
! ***JUMP OUT OF DO-LOOP
  IF (abserr <= ertest) GO TO 100
  
!      PREPARE BISECTION OF THE SMALLEST INTERVAL.
  
  70 IF (numrl2 == 1) noext = .true.
  IF (ier == 5) GO TO 100
  maxerr = iord(1)
  errmax = elist(maxerr)
  nrmax = 1
  extrap = .false.
  small = small*0.5D0
  erlarg = errsum
  CYCLE
  80 small = ABS(b - a)*0.375D0
  erlarg = errsum
  ertest = errbnd
  rlist2(2) = area
END DO

!      SET FINAL RESULT AND ERROR ESTIMATE.
!      ------------------------------------

100 IF (abserr == oflow) GO TO 115
IF (ier + ierro == 0) GO TO 110
IF (ierro == 3) abserr = abserr + correc
IF (ier == 0) ier = 3
IF (result /= 0.0D0 .AND. area /= 0.0D0) GO TO 105
IF (abserr > errsum) GO TO 115
IF (area == 0.0D0) GO TO 130
GO TO 110
105 IF (abserr/ABS(result) > errsum/ABS(area)) GO TO 115

!      TEST ON DIVERGENCE.

110 IF(ksgn == (-1) .AND. MAX(ABS(result),ABS(area)) <=  &
defabs*0.1D-01) GO TO 130
IF(0.1D-01 > (result/area) .OR. (result/area) > 0.1D+03  &
   .OR. errsum > ABS(area)) ier = 6
GO TO 130

!      COMPUTE GLOBAL INTEGRAL SUM.

115 result = SUM( rlist(1:last) )
abserr = errsum
130 IF (ier > 2) ier = ier - 1
999 RETURN
END SUBROUTINE qxgse


SUBROUTINE qxcpy (a, b, l)

!  TO COPY THE REAL VECTOR B OF LENGTH L   I N T O
!          THE REAL VECTOR A OF LENGTH L

USE constants_NSWC
IMPLICIT NONE

INTEGER, INTENT(IN)   :: l
REAL (dp), DIMENSION(:), INTENT(IN)  :: b
REAL (dp), DIMENSION(:), INTENT(OUT) :: a

a(1:l) = b(1:l)
RETURN
END SUBROUTINE qxcpy


SUBROUTINE qxlqm (f, a, b, result, abserr, resabs, resasc, vr, vs, lr, ls,  &
                  key, epmach, uflow, oflow)

!    TO COMPUTE I = INTEGRAL OF F OVER (A, B), WITH ERROR ESTIMATE
!               J = INTEGRAL OF ABS(F) OVER (A,B)

!   PARAMETERS
!    ON ENTRY
!      F      - REAL
!               FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
!               FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
!               DECLARED E X T E R N A L IN THE DRIVER PROGRAM.

!      A      - REAL
!               LOWER LIMIT OF INTEGRATION

!      B      - REAL
!               UPPER LIMIT OF INTEGRATION

!      VR     - REAL
!               VECTOR OF LENGTH LR CONTAINING THE
!               SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS

!      VS     - REAL
!               VECTOR OF LENGTH LS CONTAINING THE
!               SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS

!      LR     - INTEGER
!      LS       NUMBER OF ELEMENTS IN VR,VS RESPECTIVELY

!    KEY    - INTEGER
!             KEY FOR CHOICE OF LOCAL INTEGRATION RULE
!             RMS FORMULAS ARE USED WITH
!              13 - 19               POINTS IF KEY < 1,
!              13 - 19 - (27)        POINTS IF KEY = 1,
!              13 - 19 - (27) - (41) POINTS IF KEY = 2,
!                   19 -  27  - (41) POINTS IF KEY = 3,
!                         27  -  41  POINTS IF KEY > 3.

!             (RULES) USED IF THE FUNCTION APPEARS REGULAR ENOUGH

!      EPMACH - REAL
!               THE RELATIVE PRECISION OF THE FLOATING
!               ARITHMETIC BEING USED.

!      UFLOW  - REAL
!               THE SMALLEST POSITIVE MAGNITUDE.

!      OFLOW  - REAL
!               THE LARGEST POSITIVE MAGNITUDE.

!    ON RETURN
!      RESULT - REAL
!               APPROXIMATION TO THE INTEGRAL I

!      ABSERR - REAL
!               ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
!               WHICH SHOULD NOT EXCEED ABS(I-RESULT)

!      RESABS - REAL
!               APPROXIMATION TO THE INTEGRAL J

!      RESASC - REAL
!               APPROXIMATION TO THE INTEGRAL OF ABS(F-I/(B-A)) OVER (A,B)

!      VR     - REAL
!               VECTOR OF LENGTH LR CONTAINING THE
!               SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS

!      VS     - REAL
!               VECTOR OF LENGTH LS CONTAINING THE
!               SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS

!      LR     - INTEGER
!      LS       NUMBER OF ELEMENTS IN VR,VS RESPECTIVELY

!***ROUTINES CALLED  QXRUL

USE constants_NSWC
IMPLICIT NONE

REAL (dp), INTENT(IN)                   :: a, b, epmach, oflow, uflow
INTEGER, INTENT(IN)                     :: key
REAL (dp), INTENT(OUT)                  :: result, abserr, resabs, resasc
REAL (dp), DIMENSION(:), INTENT(IN OUT) :: vr, vs
INTEGER, INTENT(IN OUT)                 :: lr, ls

REAL (dp) :: t, resg, resk, errold
INTEGER   :: k, k0, k1, k2, key1

INTERFACE
  FUNCTION f(x) RESULT(fx)
    USE constants_NSWC
    IMPLICIT NONE
    REAL (dp), INTENT(IN) :: x
    REAL (dp)             :: fx
  END FUNCTION f
END INTERFACE

key1 = MAX(key ,  0)
key1 = MIN(key1,  4)
k0   = MAX(key1-2,0)
k1   = k0 + 1
k2   = MIN(key1+1,3)

CALL qxrul (f, a, b, resg, resabs, resasc, k0, k1, vr, vs, lr, ls)
errold = oflow
t = 10.0D0*epmach
DO k = k1, k2
  CALL qxrul (f, a, b, resk, resabs, resasc, k, k1, vr, vs, lr, ls)
  result = resk
  abserr = ABS(resk - resg)
  IF (resasc /= 0.0D0 .AND. abserr /= 0.0D0)  &
  abserr = resasc*MIN(1.0D0,(200.0D0*abserr/resasc)**1.5D0)
  IF (resabs > uflow/t) abserr = MAX(t*resabs,abserr)
  resg = resk
  IF (abserr > errold*0.16D0) EXIT
  IF (abserr < 1000.0D0*epmach*resabs) CYCLE
  errold = abserr
END DO

RETURN
END SUBROUTINE qxlqm


SUBROUTINE qxrul (f, xl, xu, y, ya, ym, ke, k1, fv1, fv2, l1, l2)

!    TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR ESTIMATE
!    AND CONDITIONALLY COMPUTE
!               J = INTEGRAL OF ABS(F) OVER (A,B)
!               BY USING AN  RMS RULE
!   PARAMETERS
!    ON ENTRY
!      F      - REAL
!               FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
!               FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
!               DECLARED E X T E R N A L IN THE DRIVER PROGRAM.

!      XL     - REAL
!               LOWER LIMIT OF INTEGRATION

!      XU     - REAL
!               UPPER LIMIT OF INTEGRATION

!      KE     - INTEGER
!             KEY FOR CHOICE OF LOCAL INTEGRATION RULE
!             AN RMS RULE IS USED WITH
!                 13      POINTS IF KE  = 2,
!                 19      POINTS IF KE  = 3,
!                 27      POINTS IF KE  = 4,
!                 42      POINTS IF KE  = 5

!      K1     INTEGER
!             VALUE OF KEY FOR WHICH THE ADDITIONAL ESTIMATES
!             YA, YM ARE TO BE COMPUTED

!      FV1    - REAL
!               VECTOR CONTAINING L1
!               SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS

!      FV2    - REAL
!               VECTOR CONTAINING L2
!               SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS

!      L1     - INTEGER
!      L2       NUMBER OF ELEMENTS IN FV1,FV2  RESPECTIVELY

!    ON RETURN
!      Y      - REAL
!               APPROXIMATION TO THE INTEGRAL I
!               RESULT IS COMPUTED BY APPLYING THE REQUESTED RMS RULE

!      YA     - REAL
!               IF KEY = K1  APPROXIMATION TO THE INTEGRAL J
!               ELSE UNCHANGED

!      YM     - REAL
!               IF KEY = K1  APPROXIMATION TO THE INTEGRAL OF
!                              ABS(F-I/(XU-XL)   OVER (XL,XU)
!               ELSE UNCHANGED

!      FV1    - REAL
!               VECTOR CONTAINING L1
!               SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS

!      FV2    - REAL
!               VECTOR CONTAINING L2
!               SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS

!      L1     - INTEGER
!      L2       NUMBER OF ELEMENTS IN FV1,FV2  RESPECTIVELY

!------------------------
USE constants_NSWC
IMPLICIT NONE

REAL (dp), INTENT(IN)                   :: xl, xu
INTEGER, INTENT(IN)                     :: ke, k1
REAL (dp), INTENT(OUT)                  :: y, ya, ym
REAL (dp), DIMENSION(:), INTENT(IN OUT) :: fv1, fv2
INTEGER, INTENT(IN OUT)                 :: l1, l2

REAL (dp) :: ldl, y2, aa, bb, c
INTEGER   :: istart(4) = (/ 0, 7, 17, 31 /), length(4) = (/ 7, 10, 14, 21 /), &
             j, i, is, k, ks

INTERFACE
  FUNCTION f(x) RESULT(fx)
    USE constants_NSWC
    IMPLICIT NONE
    REAL (dp), INTENT(IN) :: x
    REAL (dp)             :: fx
  END FUNCTION f
END INTERFACE

!------------------------
REAL (dp) :: xx(21) = (/ 0.D0,                    .25000000000000000000D+00, &
                       .50000000000000000000D+00, .75000000000000000000D+00, &
                       .87500000000000000000D+00, .93750000000000000000D+00, &
                       .10000000000000000000D+01, .37500000000000000000D+00, &
                       .62500000000000000000D+00, .96875000000000000000D+00, &
                       .12500000000000000000D+00, .68750000000000000000D+00, &
                       .81250000000000000000D+00, .98437500000000000000D+00, &
                       .18750000000000000000D+00, .31250000000000000000D+00, &
                       .43750000000000000000D+00, .56250000000000000000D+00, &
                       .84375000000000000000D+00, .90625000000000000000D+00, &
                       .99218750000000000000D+00 /)
REAL (dp) :: ww(52) = (/     &
1.303262173284849021810473057638590518409112513421D-1, 2.390632866847646220320329836544615917290026806242D-1,  &
2.630626354774670227333506083741355715758124943143D-1, 2.186819313830574175167853094864355208948886875898D-1,  &
2.757897646642836865859601197607471574336674206700D-2, 1.055750100538458443365034879086669791305550493830D-1,  &
1.571194260595182254168429283636656908546309467968D-2, 1.298751627936015783241173611320651866834051160074D-1,  &
2.249996826462523640447834514709508786970828213187D-1, 1.680415725925575286319046726692683040162290325505D-1,  &
1.415567675701225879892811622832845252125600939627D-1, 1.006482260551160175038684459742336605269707889822D-1,  &
2.510604860724282479058338820428989444699235030871D-2, 9.402964360009747110031098328922608224934320397592D-3,  &
5.542699233295875168406783695143646338274805359780D-2, 9.986735247403367525720377847755415293097913496236D-2,  &
4.507523056810492466415880450799432587809828791196D-2, 6.300942249647773931746170540321811473310938661469D-2,  &
1.261383225537664703012999637242003647020326905948D-1, 1.273864433581028272878709981850307363453523117880D-1,  &
8.576500414311820514214087864326799153427368592787D-2, 7.102884842310253397447305465997026228407227220665D-2,  &
5.026383572857942403759829860675892897279675661654D-2, 4.683670010609093810432609684738393586390722052124D-3,  &
1.235837891364555000245004813294817451524633100256D-1, 1.148933497158144016800199601785309838604146040215D-1,  &
1.252575774226122633391477702593585307254527198070D-2, 1.239572396231834242194189674243818619042280816640D-1,  &
2.501306413750310579525950767549691151739047969345D-2, 4.915957918146130094258849161350510503556792927578D-2,  &
2.259167374956474713302030584548274729936249753832D-2, 6.362762978782724559269342300509058175967124446839D-2,  &
9.950065827346794643193261975720606296171462239514D-2, 7.048220002718565366098742295389607994441704889441D-2,  &
6.512297339398335645872697307762912795346716454337D-2, 3.998229150313659724790527138690215186863915308702D-2,  &
3.456512257080287509832054272964315588028252136044D-2, 2.212167975884114432760321569298651047876071264944D-3,  &
8.140326425945938045967829319725797511040878579808D-2, 6.583213447600552906273539578430361199084485578379D-2,  &
2.592913726450792546064232192976262988065252032902D-2, 1.187141856692283347609436153545356484256869129472D-1,  &
5.999947605385971985589674757013565610751028128731D-2, 5.500937980198041736910257988346101839062581489820D-2,  &
5.264422421764655969760271538981443718440340270116D-3, 1.533126874056586959338368742803997744815413565014D-2,  &
3.527159369750123100455704702965541866345781113903D-2, 5.000556431653955124212795201196389006184693561679D-2,  &
5.744164831179720106340717579281831675999717767532D-2, 1.598823797283813438301248206397233634639162043386D-2,  &
2.635660410220884993472478832884065450876913559421D-2, 1.196003937945541091670106760660561117114584656319D-2 /)
!------------------------
k = ke + 1
is = istart(k)
ks = length(k)
ldl = xu - xl
bb = ldl*0.5D0
aa = xl + bb

y = 0.0D0
DO i = 1, ks
  c = bb*xx(i)
  IF (i > l1) fv1(i) = f(aa + c)
  IF (i > l2) fv2(i) = f(aa - c)
  j = is + i
  y = y + (fv1(i) + fv2(i))*ww(j)
END DO

y2 = y
y = y*bb
IF (l1 < ks) l1 = ks
IF (l2 < ks) l2 = ks
IF (ke /= k1) RETURN

ya = 0.0D0
DO i = 1, ks
  j = is + i
  ya = ya + (ABS(fv1(i)) + ABS(fv2(i)))*ww(j)
END DO
ya = ya*ABS(bb)

y2 = y2*0.5D0
ym = 0.0D0
DO i = 1, ks
  j = is + i
  ym = ym + (ABS(fv1(i) - y2) + ABS(fv2(i) - y2))*ww(j)
END DO
ym = ym*ABS(bb)
RETURN
END SUBROUTINE qxrul


SUBROUTINE qxrrd (f, z, lz, xl, xu, r, s, lr, ls)

!    TO REORDER THE COMPUTED FUNCTIONAL VALUES BEFORE
!    THE BISECTION OF AN INTERVAL

!   PARAMETERS
!    ON ENTRY
!      F      - REAL
!               FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
!               FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
!               DECLARED E X T E R N A L IN THE DRIVER PROGRAM.

!      XL     - REAL
!               LOWER LIMIT OF INTERVAL

!      XU     - REAL
!               UPPER LIMIT OF INTERVAL

!      Z      - REAL
!               VECTOR CONTAINING LZ SAVED FUNCTIONAL VALUES

!      LZ     - INTEGER
!               NUMBER OF ELEMENTS IN LZ

!    ON RETURN
!      R      - REAL
!      S        VECTORS CONTAINING LR, LS
!               SAVED  FUNCTIONAL VALUES FOR THE NEW INTERVALS

!      LR     - INTEGER
!      LS       NUMBER OF ELEMENTES IN R,S RESPECTIVELY

!***ROUTINES CALLED  F
USE constants_NSWC
IMPLICIT NONE

REAL (dp), DIMENSION(:), INTENT(IN)  :: z
REAL (dp), INTENT(IN)                :: xl, xu
INTEGER, INTENT(IN)                  :: lz
REAL (dp), DIMENSION(:), INTENT(OUT) :: r, s
INTEGER, INTENT(OUT)                 :: lr, ls

REAL (dp) :: dlen, centr

INTERFACE
  FUNCTION f(x) RESULT(fx)
    USE constants_NSWC
    IMPLICIT NONE
    REAL (dp), INTENT(IN) :: x
    REAL (dp)             :: fx
  END FUNCTION f
END INTERFACE

dlen = (xu - xl)*0.5D0
centr = xl + dlen
r(1) =  z(3)
r(2) =  z(9)
r(3) =  z(4)
r(4) =  z(5)
r(5) =  z(6)
r(6) =  z(10)
r(7) =  z(7)
s(1) =  z(3)
s(2) =  z(8)
s(3) =  z(2)
s(7) =  z(1)
IF (lz > 11) GO TO 10

r(8) =  f(centr + dlen*0.375D0)
r(9) =  f(centr + dlen*0.625D0)
r(10) = f(centr + dlen*0.96875D0)
lr = 10
IF (lz /= 11) s(4) = f(centr - dlen*0.75D0)
IF (lz == 11) s(4) = z(11)
s(5) =  f(centr - dlen*0.875D0)
s(6) =  f(centr - dlen*0.9375D0)
s(8) =  f(centr - dlen*0.375D0)
s(9) =  f(centr - dlen*0.625D0)
s(10) = f(centr - dlen*0.96875D0)
ls = 10
RETURN

10 r(8) = z(12)
r(9) = z(13)
r(10) = z(14)
lr = 10
s(4) = z(11)
s(5) = f(centr - dlen*0.875D0)
s(6) = f(centr - dlen*0.9375D0)
IF (lz > 14) GO TO 20
s(8)  = f(centr - dlen*0.375D0)
s(9)  = f(centr - dlen*0.625D0)
s(10) = f(centr - dlen*0.96875D0)
ls = 10
RETURN

20 r(11) = z(18)
r(12) = z(19)
r(13) = z(20)
r(14) = z(21)
lr = 14
s(8) = z(16)
s(9) = z(15)
s(10) = f(centr - dlen*0.96875D0)
s(11) = z(17)
ls = 11
RETURN
END SUBROUTINE qxrrd

END MODULE adapt_quad
