Function Brm_Dmlin, a, b, maxfcn, rerr, eph, gamlow, gambrk, gamhigh, $
    p, q, z, EFD, ier

l = Where(a)

;   Create arrays for integral sum and error flags.

intsum = DblArr(N_Elements(a))
ier = IntArr(N_Elements(a))

;   Maximum possible order of the Gaussian quadrature integration is 2^nlim.

nlim = 12

FOR ires = 1, nlim DO BEGIN

    npoint = 2l^ires
    IF (npoint GT maxfcn) THEN GOTO, noconv

    eph1 = eph(l)

    ;   Create arrays for abscissas and weights.

    x = DblArr(N_Elements(l),npoint)
    w = DblArr(N_Elements(l),npoint)

    ;   Use Brm_GauLeg to perform the outer integration.


    Brm_GauLeg, a(l), b, npoint, x, w


    lastsum = intsum

    intsum(l) = (w*Brm_FThin(x, eph1, gamlow, gambrk, gamhigh, p, q, z, EFD)) $
        # (Make_Array(npoint, Value=1))

    ;   Use last two calculations to check for convergence.

    l1 = Abs(intsum-lastsum)
    l2 = rerr*Abs(intsum)
    l = Where(l1 GT l2)

    IF (l(0) EQ -1) THEN GOTO, conv

ENDFOR

;   ier = 120 indicates no convergence.

noconv: ier(l) = 120

;   ier = 0 indicates convergence.

conv:

Brm_Dmlin = intsum

Return, Brm_Dmlin
Return, ier

EN