    !Module of generally useful routines and definitions
    !Antony Lewis, http://cosmologist.info/

    !April 2006: fix to TList_RealArr_Thin
    !March 2008: fix to Ranges
    !June 2010: fixed bug in  DeleteFile gradually using up file units

    module Ranges
    !A collection of ranges, consisting of sections of minimum step size
    implicit none

    integer, parameter :: Max_Ranges = 100
    double precision, parameter :: RangeTol = 0.1d0
    !fraction of bin width we are prepared for merged bin widths to increase by

    Type Region
        integer start_index
        integer steps
        logical :: IsLog
        double precision Low, High
        double precision delta
        double precision delta_max, delta_min !for log spacing, the non-log max and min step size
    end Type Region

    Type Regions

        integer count
        integer npoints
        double precision Lowest, Highest
        Type(Region) :: R(Max_ranges)
        logical :: has_dpoints
        double precision, dimension(:), pointer :: points => null(), dpoints => null()
        !dpoints is (points(i+1)-points(i-1))/2

    end Type Regions

    contains

    subroutine Ranges_Init(R)
    Type(Regions) R

    call Ranges_Free(R)

    end  subroutine Ranges_Init

    subroutine Ranges_Free(R)
    Type(Regions) R
    integer status

    deallocate(R%points,stat = status)
    deallocate(R%dpoints,stat = status)
    call Ranges_Nullify(R)

    end  subroutine Ranges_Free


    subroutine Ranges_Nullify(R)
    Type(Regions) R

    nullify(R%points)
    nullify(R%dpoints)
    R%count = 0
    R%npoints = 0
    R%has_dpoints = .false.


    end subroutine Ranges_Nullify

    subroutine Ranges_Assign(R,Rin)
    Type(Regions) R, Rin

    call Ranges_Init(R)
    R = Rin
    nullify(R%points,R%dpoints)
    if (associated(Rin%points)) then
        call Ranges_GetArray(R, associated(Rin%dpoints))
    end if

    end subroutine Ranges_Assign

    function Ranges_IndexOf(Reg, tau) result(pointstep)
    Type(Regions), intent(in), target :: Reg
    Type(Region), pointer :: AReg
    double precision :: tau
    integer pointstep
    integer i


    pointstep=0
    do i=1,Reg%count
        AReg => Reg%R(i)

        if (tau < AReg%High .and. tau >= AReg%Low) then
            if (AReg%IsLog) then
                pointstep = AReg%start_index + int( log(tau/AReg%Low)/AReg%delta)
            else
                pointstep = AReg%start_index + int(( tau - AReg%Low)/AReg%delta)
            end if
            return
        end if

    end do

    if (tau >= Reg%Highest) then
        pointstep = Reg%npoints
    else
        write (*,*) 'Ranges_IndexOf: value out of range'
        stop
    end if

    end function Ranges_IndexOf


    subroutine Ranges_GetArray(Reg, want_dpoints)
    Type(Regions), target :: Reg
    Type(Region), pointer :: AReg
    logical, intent(in), optional :: want_dpoints
    integer status,i,j,ix


    if (present(want_dpoints)) then
        Reg%has_dpoints = want_dpoints
    else
        Reg%has_dpoints = .true.
    end if

    deallocate(Reg%points,stat = status)
    allocate(Reg%points(Reg%npoints))

    ix=0
    do i=1, Reg%count
        AReg => Reg%R(i)
        do j = 0, AReg%steps-1
            ix=ix+1
            if (AReg%IsLog) then
                Reg%points(ix) = AReg%Low*exp(j*AReg%delta)
            else
                Reg%points(ix) = AReg%Low + AReg%delta*j
            end if
        end do
    end do
    ix =ix+1
    Reg%points(ix) = Reg%Highest
    if (ix /= Reg%npoints) stop 'Ranges_GetArray: ERROR'

    if (Reg%has_dpoints) call Ranges_Getdpoints(Reg)

    end subroutine Ranges_GetArray


    subroutine Ranges_Getdpoints(Reg, half_ends)
    Type(Regions), target :: Reg
    logical, intent(in), optional :: half_ends
    integer i, status
    logical halfs

    if (present(half_ends)) then
        halfs = half_ends
    else
        halfs = .true.
    end if

    deallocate(Reg%dpoints,stat = status)
    allocate(Reg%dpoints(Reg%npoints))

    do i=2, Reg%npoints-1
        Reg%dpoints(i) = (Reg%points(i+1) - Reg%points(i-1))/2
    end do
    if (halfs) then
        Reg%dpoints(1) = (Reg%points(2) - Reg%points(1))/2
        Reg%dpoints(Reg%npoints) = (Reg%points(Reg%npoints) - Reg%points(Reg%npoints-1))/2
    else
        Reg%dpoints(1) = (Reg%points(2) - Reg%points(1))
        Reg%dpoints(Reg%npoints) = (Reg%points(Reg%npoints) - Reg%points(Reg%npoints-1))
    end if
    end subroutine Ranges_Getdpoints


    subroutine Ranges_Add_delta(Reg, t_start, t_end, t_approx_delta, IsLog)
    Type(Regions), target :: Reg
    logical, intent(in), optional :: IsLog
    double precision, intent(in) :: t_start, t_end, t_approx_delta
    integer n
    logical :: WantLog

    if (present(IsLog)) then
        WantLog = IsLog
    else
        WantLog = .false.
    end if

    if (t_end <= t_start) &
        stop 'Ranges_Add_delta: end must be larger than start'
    if (t_approx_delta <=0) stop 'Ranges_Add_delta: delta must be > 0'

    if (WantLog) then
        n  = max(1,int(log(t_end/t_start)/t_approx_delta + 1.d0 - RangeTol))
    else
        n  = max(1,int((t_end-t_start)/t_approx_delta + 1.d0 - RangeTol))
    end if
    call Ranges_Add(Reg,t_start, t_end, n, WantLog)

    end subroutine Ranges_Add_delta


    subroutine Ranges_Add(Reg, t_start, t_end, nstep, IsLog)
    Type(Regions), target :: Reg
    logical, intent(in), optional :: IsLog
    double precision, intent(in) :: t_start, t_end
    integer, intent(in) :: nstep
    Type(Region), pointer :: AReg, LastReg
    Type(Region), target :: NewRegions(Max_Ranges)
    double precision EndPoints(0:Max_Ranges*2)
    integer ixin, nreg, ix, i,j, nsteps
    double precision delta
    logical WantLog
    double precision min_request, max_request, min_log_step, max_log_step, diff, max_delta
    double precision RequestDelta(Max_Ranges)

    if (present(IsLog)) then
        WantLog = IsLog
    else
        WantLog = .false.
    end if

    if (WantLog) then
        delta = log(t_end/t_start) / nstep
    else
        delta = (t_end - t_start) / nstep
    end if

    if (t_end <= t_start) stop 'Ranges_Add: end must be larger than start'
    if (nstep <=0) stop 'Ranges_Add: nstep must be > 0'
    if (Reg%Count>= Max_Ranges) stop 'Ranges_Add: Increase Max_Ranges'

    !avoid IBM compiler bug, from Angel de Vicente
    !    if (Reg%count > 0) NewRegions(1:Reg%count) = Reg%R(1:Reg%count)
    if (Reg%count > 0) THEN
        DO i=1,Reg%count
            NewRegions(i) = Reg%R(i)
        END DO
    END IF
    nreg = Reg%count + 1
    AReg=> NewRegions(nreg)
    AReg%Low = t_start
    AReg%High = t_end
    AReg%delta = delta
    AReg%steps = nstep
    AReg%IsLog = WantLog

    !Get end point in order
    ix = 0
    do i=1, nreg

        AReg => NewRegions(i)
        if (ix==0) then
            ix = 1
            EndPoints(ix) = AReg%Low
            ix = 2
            EndPoints(ix) = AReg%High
        else
            ixin = ix
            do j=1,ixin
                if (AReg%Low < EndPoints(j)) then
                    EndPoints(j+1:ix+1) = EndPoints(j:ix)
                    EndPoints(j) = AReg%Low
                    ix=ix+1
                    exit
                end if
            end do
            if (ixin == ix) then
                ix = ix+1
                EndPoints(ix) = AReg%Low
                ix = ix+1
                EndPoints(ix) = AReg%High
            else
                ixin = ix
                do j=1,ixin
                    if (AReg%High < EndPoints(j)) then
                        EndPoints(j+1:ix+1) = EndPoints(j:ix)
                        EndPoints(j) = AReg%High
                        ix=ix+1
                        exit
                    end if
                end do
                if (ixin == ix) then
                    ix = ix+1
                    EndPoints(ix) = AReg%High
                end if

            end if
        end if

    end do

    !remove duplicate points
    ixin = ix
    ix = 1
    do i=2, ixin
        if (EndPoints(i) /= EndPoints(ix)) then
            ix=ix+1
            EndPoints(ix) = EndPoints(i)
        end if
    end do


    !ix is the number of end points
    Reg%Lowest = EndPoints(1)
    Reg%Highest = EndPoints(ix)
    Reg%count = 0

    max_delta = Reg%Highest - Reg%Lowest

    do i=1, ix - 1
        AReg => Reg%R(i)
        AReg%Low = EndPoints(i)
        AReg%High = EndPoints(i+1)

        !          max_delta = EndPoints(i+1) - EndPoints(i)
        delta = max_delta
        AReg%IsLog = .false.

        do j=1, nreg
            if (AReg%Low >= NewRegions(j)%Low .and. Areg%Low < NewRegions(j)%High) then
                if (NewRegions(j)%IsLog) then
                    if (AReg%IsLog) then
                        delta = min(delta,NewRegions(j)%delta)
                    else
                        min_log_step = AReg%Low*(exp(NewRegions(j)%delta)-1)
                        if (min_log_step < delta) then
                            max_log_step = AReg%High*(1-exp(-NewRegions(j)%delta))
                            if  (delta < max_log_step) then
                                delta = min_log_step
                            else
                                AReg%IsLog = .true.
                                delta = NewRegions(j)%delta
                            end if
                        end if
                    end if
                else !NewRegion is not log
                    if (AReg%IsLog) then
                        max_log_step = AReg%High*(1-exp(-delta))
                        if (NewRegions(j)%delta < max_log_step) then
                            min_log_step = AReg%Low*(exp(delta)-1)
                            if (min_log_step <  NewRegions(j)%delta) then
                                AReg%IsLog = .false.
                                delta =  min_log_step
                            else
                                delta = - log(1- NewRegions(j)%delta/AReg%High)
                            end if
                        end if
                    else
                        delta = min(delta, NewRegions(j)%delta)
                    end if
                end if
            end if
        end do

        if (AReg%IsLog) then
            Diff = log(AReg%High/AReg%Low)
        else
            Diff = AReg%High - AReg%Low
        endif
        if (delta >= Diff) then
            AReg%delta = Diff
            AReg%steps = 1
        else
            AReg%steps  = max(1,int(Diff/delta + 1.d0 - RangeTol))
            AReg%delta = Diff / AReg%steps
        end if

        Reg%count = Reg%count + 1
        RequestDelta(Reg%Count) = delta

        if (AReg%IsLog) then
            if (AReg%steps ==1) then
                AReg%Delta_min = AReg%High - AReg%Low
                AReg%Delta_max = AReg%Delta_min
            else
                AReg%Delta_min = AReg%Low*(exp(AReg%delta)-1)
                AReg%Delta_max = AReg%High*(1-exp(-AReg%delta))
            end if
        else
            AReg%Delta_max = AReg%delta
            AReg%Delta_min = AReg%delta
        end if
    end do


    !Get rid of tiny regions
    ix = Reg%Count
    do i=ix, 1, -1
        AReg => Reg%R(i)
        if (AReg%steps ==1) then
            Diff = AReg%High - AReg%Low
            if (AReg%IsLog) then
                min_request = AReg%Low*(exp(RequestDelta(i))-1)
                max_request = AReg%High*(1-exp(-RequestDelta(i)))
            else
                min_request = RequestDelta(i)
                max_request = min_request
            end if
            if (i/= Reg%Count) then  !from i/= ix Mar08
                LastReg => Reg%R(i+1)
                if (RequestDelta(i) >= AReg%delta .and. Diff <= LastReg%Delta_min &
                    .and. LastReg%Delta_min <= max_request) then

                LastReg%Low = AReg%Low
                if (Diff > LastReg%Delta_min*RangeTol) then
                    LastReg%steps =  LastReg%steps + 1
                end if
                if (LastReg%IsLog) then
                    LastReg%delta = log(LastReg%High/LastReg%Low) / LastReg%steps
                else
                    LastReg%delta = (LastReg%High -LastReg%Low) / LastReg%steps
                end if
                Reg%R(i:Reg%Count-1) = Reg%R(i+1:Reg%Count)
                Reg%Count = Reg%Count -1
                cycle
                end if
            end if
            if (i/=1) then
                LastReg => Reg%R(i-1)
                if (RequestDelta(i) >= AReg%delta .and. Diff <= LastReg%Delta_max &
                    .and. LastReg%Delta_max <= min_request) then
                LastReg%High = AReg%High
                !AlMat08 LastReg%Low = AReg%Low
                if (Diff > LastReg%Delta_max*RangeTol) then
                    LastReg%steps =  LastReg%steps + 1
                end if
                if (LastReg%IsLog) then
                    LastReg%delta = log(LastReg%High/LastReg%Low) / LastReg%steps
                else
                    LastReg%delta = (LastReg%High -LastReg%Low) / LastReg%steps
                end if
                Reg%R(i:Reg%Count-1) = Reg%R(i+1:Reg%Count)
                Reg%Count = Reg%Count -1
                end if
            end if
        end if
    end do


    !Set up start indices and get total number of steps
    nsteps = 1
    do i = 1, Reg%Count
        AReg => Reg%R(i)
        AReg%Start_index = nsteps
        nsteps = nsteps + AReg%steps
        if (AReg%IsLog) then
            if (AReg%steps ==1) then
                AReg%Delta_min = AReg%High - AReg%Low
                AReg%Delta_max = AReg%Delta_min
            else
                AReg%Delta_min = AReg%Low*(exp(AReg%delta)-1)
                AReg%Delta_max = AReg%High*(1-exp(-AReg%delta))
            end if
        else
            AReg%Delta_max = AReg%delta
            AReg%Delta_min = AReg%delta
        end if
    end do

    Reg%npoints = nsteps

    end subroutine Ranges_Add


    subroutine Ranges_Write(Reg)
    Type(Regions), intent(in), target :: Reg
    Type(Region), pointer :: AReg
    integer i

    do i=1,Reg%count
        AReg => Reg%R(i)
        if (AReg%IsLog) then
            Write (*,'("Range ",I3,":", 3E14.4," log")') i, AReg%Low, AReg%High, AReg%delta
        else
            Write (*,'("Range ",I3,":", 3E14.4," linear")') i, AReg%Low, AReg%High, AReg%delta
        end if
    end do
    end subroutine Ranges_Write


    end module Ranges


    module Lists
    !Currently implements lists of strings and lists of arrays of reals
    implicit none

    type real_pointer
        real, dimension(:), pointer :: p
    end type real_pointer

    type double_pointer
        double precision, dimension(:), pointer :: p
    end type double_pointer

    type String_pointer
        character, dimension(:), pointer :: p
    end type String_pointer


    Type TList_RealArr
        integer Count
        integer Delta
        integer Capacity
        type(Real_Pointer), dimension(:), pointer :: Items
    end Type TList_RealArr

    Type TStringList
        integer Count
        integer Delta
        integer Capacity
        type(String_Pointer), dimension(:), pointer :: Items
    end Type TStringList

    contains

    subroutine TList_RealArr_Init(L)
    Type (TList_RealArr) :: L

    L%Count = 0
    L%Capacity = 0
    L%Delta = 1024
    nullify(L%items)

    end subroutine TList_RealArr_Init

    subroutine TList_RealArr_Clear(L)
    Type (TList_RealArr) :: L
    integer i, status

    do i=L%Count,1,-1
        deallocate (L%Items(i)%P, stat = status)
    end do
    deallocate (L%Items, stat = status)
    nullify(L%Items)
    L%Count = 0
    L%Capacity = 0

    end subroutine TList_RealArr_Clear


    subroutine TList_RealArr_Add(L, P)
    Type (TList_RealArr) :: L
    real, intent(in) :: P(:)
    integer s

    if (L%Count == L%Capacity) call TList_RealArr_SetCapacity(L, L%Capacity + L%Delta)
    s = size(P)
    L%Count = L%Count + 1
    allocate(L%Items(L%Count)%P(s))
    L%Items(L%Count)%P = P

    end subroutine TList_RealArr_Add

    subroutine TList_RealArr_SetCapacity(L, C)
    Type (TList_RealArr) :: L
    integer C
    type(Real_Pointer), dimension(:), pointer :: TmpItems

    if (L%Count > 0) then
        if (C < L%Count) stop 'TList_RealArr_SetCapacity: smaller than Count'
        allocate(TmpItems(L%Count))
        TmpItems = L%Items(1:L%Count)
        deallocate(L%Items)
        allocate(L%Items(C))
        L%Items(1:L%Count) = TmpItems
        deallocate(TmpItems)
    else
        allocate(L%Items(C))
    end if
    L%Capacity = C
    end subroutine TList_RealArr_SetCapacity

    subroutine TList_RealArr_Delete(L, i)
    Type (TList_RealArr) :: L
    integer, intent(in) :: i
    integer status

    deallocate(L%items(i)%P, stat = status)
    if (L%Count > 1) L%Items(i:L%Count-1) = L%Items(i+1:L%Count)
    L%Count = L%Count -1

    end subroutine TList_RealArr_Delete

    subroutine TList_RealArr_SaveBinary(L,fid)
    Type (TList_RealArr) :: L
    integer, intent(in) :: fid
    integer i

    write (fid) L%Count
    do i=1,L%Count
        write(fid) size(L%Items(i)%P)
        write(fid) L%Items(i)%P
    end do

    end subroutine TList_RealArr_SaveBinary

    subroutine TList_RealArr_ReadBinary(L,fid)
    Type (TList_RealArr) :: L
    integer, intent(in) :: fid
    integer num,i,sz

    call TList_RealArr_Clear(L)
    read (fid) num
    call TList_RealArr_SetCapacity(L, num)
    do i=1,num
        read(fid) sz
        allocate(L%Items(i)%P(sz))
        read(fid) L%Items(i)%P
    end do
    L%Count = num

    end subroutine TList_RealArr_ReadBinary


    subroutine TList_RealArr_Thin(L, i)
    Type (TList_RealArr) :: L
    integer, intent(in) :: i
    integer newCount
    type(Real_Pointer), dimension(:), pointer :: TmpItems

    if (L%Count > 1) then
        newCount = (L%Count-1)/i+1
        allocate(TmpItems(newCount))
        TmpItems = L%Items(1:L%Count:i)
        deallocate(L%Items)
        L%Capacity = newCount
        allocate(L%Items(L%Capacity))
        L%Items = TmpItems
        L%Count = newCount
        deallocate(TmpItems)
    end if
    end subroutine TList_RealArr_Thin

    subroutine TList_RealArr_ConfidVal(L, ix, limfrac, ix1, ix2, Lower, Upper)
    !Taking the ix'th entry in each array to be a sample, value for which
    !limfrac of the items between ix1 and ix2 (inc) are above or below
    !e.g. if limfrac = 0.05 get two tail 90% confidence limits
    Type (TList_RealArr) :: L
    integer, intent(IN) :: ix
    real, intent(IN) :: limfrac
    real, intent(OUT), optional :: Lower, Upper
    integer, intent(IN), optional :: ix1,ix2
    integer b,t,samps
    real pos, d
    type(Real_Pointer), dimension(:), pointer :: SortItems

    b=1
    t=L%Count
    if (present(ix1)) b = ix1
    if (present(ix2)) t = ix2
    samps = t - b + 1

    allocate(SortItems(samps))
    SortItems = L%Items(b:t)
    call QuickSortArr_Real(SortItems, 1, samps, ix)
    if (present(Lower)) then
        pos = (samps-1)*limfrac + 1
        b = max(int(pos),1)
        Lower = SortItems(b)%P(ix)
        if (b < samps .and. pos>b) then
            d = pos - b
            Lower = Lower*(1 - d) + d * SortItems(b+1)%P(ix)
        end if
    end if
    if (present(Upper)) then
        pos = (samps-1)*(1.-limfrac) + 1
        b = max(int(pos),1)
        Upper = SortItems(b)%P(ix)
        if (b < samps .and. pos>b) then
            d = pos - b
            Upper = Upper*(1 - d) + d * SortItems(b+1)%P(ix)
        end if
    end if

    deallocate(SortItems)

    end subroutine TList_RealArr_ConfidVal

    subroutine TStringList_Init(L)
    Type (TStringList) :: L

    L%Count = 0
    L%Capacity = 0
    L%Delta = 128
    nullify(L%items)

    end subroutine TStringList_Init

    subroutine TStringList_Clear(L)
    Type (TStringList) :: L
    integer i, status

    do i=L%Count,1,-1
        deallocate (L%Items(i)%P, stat = status)
    end do
    deallocate (L%Items, stat = status)
    call TStringList_Init(L)

    end subroutine TStringList_Clear

    subroutine TStringList_SetFromString(L, S, valid_chars_in)
    Type (TStringList) :: L
    character(Len=*), intent(in) :: S
    character(Len=*), intent(in), optional :: valid_chars_in
    character(LEN=1024) item
    integer i,j
    character(LEN=256) valid_chars

    if (present(valid_chars_in)) then
        valid_chars = valid_chars_in
    else
        valid_chars='abcdefghijklmopqrstuvwxyzABCDEFGHIJKLMOPQRSTUVWXYZ0123456789_-.'
    endif

    call TStringList_Clear(L)
    item =''
    j=0
    do i=1, len_trim(S)
        if (verify(S(i:i),trim(valid_chars)) == 0) then
            j=j+1
            item(j:j) = S(i:i)
        else
            if (trim(S(i:i))/='') then
                write (*,*) 'Invalid character in: '//trim(S)
            end if
            if (j>0) call TStringList_Add(L, item(1:j))
            j=0
        end if
    end do
    if (j>0) call TStringList_Add(L, item(1:j))

    end subroutine TStringList_SetFromString



    subroutine TStringList_Add(L, P)
    Type (TStringList) :: L
    character(LEN=*), intent(in) :: P
    integer s,i

    if (L%Count == L%Capacity) call TStringList_SetCapacity(L, L%Capacity + L%Delta)
    s = len_trim(P)
    L%Count = L%Count + 1
    allocate(L%Items(L%Count)%P(s))
    do i=1,s
        L%Items(L%Count)%P(i) = P(i:i)
    end do
    end subroutine TStringList_Add

    function TStringList_Item(L, i) result(S)
    Type (TStringList) :: L
    integer, intent(in) :: i
    integer j
    character(LEN=1024) S

    S=''
    if (i<=L%Count .and. i>0) then
        do j=1,size(L%Items(i)%P)
            S(j:j)=L%Items(i)%P(j)
        end do
    end if
    end function TStringList_Item

    subroutine TStringList_SetCapacity(L, C)
    Type (TStringList) :: L
    integer C
    type(String_Pointer), dimension(:), pointer :: TmpItems

    if (L%Count > 0) then
        if (C < L%Count) stop 'TStringList_SetCapacity: smaller than Count'
        allocate(TmpItems(L%Count))
        TmpItems = L%Items(1:L%Count)
        deallocate(L%Items)
        allocate(L%Items(C))
        L%Items(1:L%Count) = TmpItems
        deallocate(TmpItems)
    else
        allocate(L%Items(C))
    end if
    L%Capacity = C

    end subroutine TStringList_SetCapacity

    subroutine TStringList_Delete(L, i)
    Type (TStringList) :: L
    integer, intent(in) :: i
    integer status

    deallocate(L%items(i)%P, stat = status)
    if (L%Count > 1) L%Items(i:L%Count-1) = L%Items(i+1:L%Count)
    L%Count = L%Count -1

    end subroutine TStringList_Delete

    function TStringList_IndexOf(L, S)
    Type (TStringList) :: L
    character(LEN=*), intent(in) :: S
    integer TStringList_IndexOf, i, j,slen

    slen = len_trim(S)
    do i=1,L%Count
        if ( size(L%Items(i)%P)==slen) then
            !Yes, comparing strings and pointer strings really is this horrible...
            j=1
            do while (L%Items(i)%P(j)==S(j:j))
                j=j+1
                if (j>slen) then
                    TStringList_IndexOf = i
                    return
                end if
            end do
        end if
    end do
    TStringList_IndexOf=-1

    end function TStringList_IndexOf


    recursive subroutine QuickSortArr_Real(Arr, Lin, R, index)
    !Sorts an array of pointers to arrays of reals by the value of the index'th entry
    integer, intent(in) :: Lin, R, index
#ifdef __GFORTRAN__
    type(real_pointer), dimension(:) :: Arr
#else
    type(real_pointer), dimension(*) :: Arr
#endif
    integer I, J, L
    real P
    type(real_pointer) :: temp

    L = Lin
    do

        I = L
        J = R
        P = Arr((L + R)/2)%p(index)

        do
            do while (Arr(I)%p(index) <  P)
                I = I + 1
            end do

            do while (Arr(J)%p(index) > P)
                J = J - 1
            end do

            if (I <= J) then

                Temp%p => Arr(I)%p
                Arr(I)%p => Arr(J)%p
                Arr(J)%p => Temp%p
                I = I + 1
                J = J - 1
            end if
            if (I > J) exit

        end do
        if (L < J) call QuickSortArr_Real(Arr, L, J, index);
        L = I
        if (I >= R) exit
    end do

    end subroutine QuickSortArr_Real



    recursive subroutine QuickSortArr(Arr, Lin, R, index)
    !Sorts an array of pointers to arrays of reals by the value of the index'th entry
    integer, intent(in) :: Lin, R, index
#ifdef __GFORTRAN__
    type(double_pointer), dimension(:) :: Arr
#else
    type(double_pointer), dimension(*) :: Arr
#endif
    integer I, J, L
    double precision P
    type(double_pointer) :: temp

    L = Lin
    do

        I = L
        J = R
        P = Arr((L + R)/2)%p(index)

        do
            do while (Arr(I)%p(index) <  P)
                I = I + 1
            end do

            do while (Arr(J)%p(index) > P)
                J = J - 1
            end do

            if (I <= J) then

                Temp%p => Arr(I)%p
                Arr(I)%p => Arr(J)%p
                Arr(J)%p => Temp%p
                I = I + 1
                J = J - 1
            end if
            if (I > J) exit

        end do
        if (L < J) call QuickSortArr(Arr, L, J, index);
        L = I
        if (I >= R) exit
    end do

    end subroutine QuickSortArr


    end module Lists

    module AMLutils
    use Lists

#ifdef DECONLY
    !Comment out if linking to LAPACK/MKL separetly
    !CXML only has LAPACK 2.0
    include 'CXML_INCLUDE.F90'
#endif


#ifdef NAGF95
    use F90_UNIX
#endif

    implicit none

#ifndef NAGF95
#ifndef GFC
#ifndef __INTEL_COMPILER_BUILD_DATE
#ifndef __GFORTRAN__
    integer iargc
    external iargc
#endif        
#endif
#endif
#endif


#ifdef MPI
    include "mpif.h"
#endif

    integer :: Feedback = 1
    integer, parameter :: tmp_file_unit = 50


    double precision, parameter :: pi=3.14159265358979323846264338328d0, &
        twopi=2*pi, fourpi=4*pi
    double precision, parameter :: root2 = 1.41421356237309504880168872421d0, sqrt2 = root2
    double precision, parameter :: log2 = 0.693147180559945309417232121458d0

    real, parameter :: pi_r = 3.141592653, twopi_r = 2*pi_r, fourpi_r = twopi_r*2

    logical :: flush_write = .true.
    !True means no data lost on crashes, but may make it slower

    integer, parameter :: file_units_start = 20
    integer, parameter :: file_units_end = 100

    logical file_units(file_units_start:file_units_end)

    INTERFACE CONCAT
    module procedure concat_s, concat_s_n
    END INTERFACE

    INTERFACE RealToStr
    module procedure SingleToStr, DoubleToStr
    END INTERFACE RealToStr

    contains

    function new_file_unit()
    integer i, new_file_unit
    logical, save :: file_units_inited = .false.
    logical notfree

    if (.not. file_units_inited) then
        file_units = .false.
        file_units_inited = .true.
    end if

    do i=file_units_start, file_units_end
        if (.not. file_units(i) .and. i/=tmp_file_unit) then
            inquire(i,opened=notfree)
            if (notfree) cycle
            file_units(i)=.true.
            new_file_unit = i
            return
        end if
    end do

    call mpiStop('No unused file unit numbers')

    end function new_file_unit


    subroutine CloseFile(i)
    integer, intent(in) :: i

    close(i)
    file_units(i) = .false.

    end subroutine CloseFile

    subroutine ClearFileUnit(i)
    integer, intent(in) :: i

    file_units(i) = .false.

    end subroutine ClearFileUnit

    function GetParamCount()
    integer GetParamCount

    GetParamCount = iargc()

    end function GetParamCount

    function GetMpiRank()
    integer GetMpiRank
#ifdef MPI 
    integer ierror
    call mpi_comm_rank(mpi_comm_world,GetMPIrank,ierror)
#else
    GetMpiRank=0
#endif    

    end function GetMpiRank

    function IsMainMPI()
    logical IsMainMPI

    IsMainMPI =  GetMpiRank() == 0

    end function IsMainMPI

    subroutine MpiStop(Msg)
    character(LEN=*), intent(in), optional :: Msg
    integer i
#ifdef MPI 
    integer ierror, MpiRank
#endif

    if (present(Msg)) write(*,*) trim(Msg)

#ifdef MPI
    call mpi_comm_rank(mpi_comm_world,MPIrank,ierror)
    write (*,*) 'MpiStop: ', MpiRank
    call MPI_ABORT(MPI_COMM_WORLD,i)
#endif
    i=1     !put breakpoint on this line to debug
    error stop

    end subroutine MpiStop

    subroutine MpiStat(MpiID, MpiSize)
    implicit none
    integer MpiID,MpiSize
#ifdef MPI  
    integer ierror
    call mpi_comm_rank(mpi_comm_world,MpiID,ierror)
    if (ierror/=MPI_SUCCESS) stop 'MpiStat: MPI rank'
    call mpi_comm_size(mpi_comm_world,MpiSize,ierror)
#else
    MpiID=0
    MpiSize=1
#endif
    end subroutine MpiStat

    subroutine MpiQuietWait
    !Set MPI thread to sleep, e.g. so can run openmp on cpu instead
#ifdef MPI  
    integer flag, ierr, STATUS(MPI_STATUS_SIZE)
    integer i, MpiId, MpiSize

    call MpiStat(MpiID, MpiSize)
    if (MpiID/=0) then
        do
            call MPI_IPROBE(0,0,MPI_COMM_WORLD,flag, MPI_STATUS_IGNORE,ierr)
            if (flag/=0) then
                call MPI_RECV(i,1,MPI_INTEGER, 0,0,MPI_COMM_WORLD,status,ierr)
                exit
            end if
            call sleep(1)
        end do
    end if
#endif
    end subroutine

    subroutine MpiWakeQuietWait
#ifdef MPI  
    integer j, MpiId, MpiSize, ierr,r

    call MpiStat(MpiID, MpiSize)
    if (MpiID==0) then
        do j=1, MpiSize-1
            call MPI_ISSEND(MpiId,1,MPI_INTEGER, j,0,MPI_COMM_WORLD,r,ierr)
        end do
    end if
#endif
    end subroutine MpiWakeQuietWait

#ifdef __GFORTRAN__

    ! ===========================================================
    function iargc ()
    ! ===========================================================
    integer iargc
    ! ===========================================================

    iargc=command_argument_count()
    end function iargc

    ! ===========================================================
    subroutine getarg(num, res)
    ! ===========================================================
    integer, intent(in) :: num
    character(len=*), intent(out) :: res
    integer l, err
    ! ===========================================================
    call get_command_argument(num,res,l,err)
    end subroutine getarg

#endif


    function GetParam(i)

    character(LEN=512) GetParam
    integer, intent(in) :: i

    if (iargc() < i) then
        GetParam = ''
    else
        call getarg(i,GetParam)
    end if
    end function GetParam

    function concat_s(S1,S2,S3,S4,S5,S6,S7,S8) result(concat)
    character(LEN=*), intent(in) :: S1, S2
    character(LEN=*), intent(in) , optional :: S3, S4, S5, S6,S7,S8
    character(LEN = 1000) concat

    concat = trim(S1) // S2
    if (present(S3)) then
        concat = trim(concat) // S3
        if (present(S4)) then
            concat = trim(concat) // S4
            if (present(S5)) then
                concat = trim(concat) // S5
                if (present(S6)) then
                    concat = trim(concat) // S6
                    if (present(S7)) then
                        concat = trim(concat) // S7
                        if (present(S8)) then
                            concat = trim(concat) // S8
                        end if
                    end if
                end if
            end if
        end if
    end if

    end function concat_s

    function concat_s_n(SS1,N2,SS3,N4,SS5,N6,SS7,N8,SS9,N10,SS11) result(concat)
    character(LEN=*), intent(in) :: SS1
    integer, intent(in) :: N2
    character(LEN=*), intent(in) , optional :: SS3, SS5, SS7, SS9,SS11
    integer, intent(in), optional ::N4,N6,N8, N10
    character(LEN = 1000) concat

    concat = trim(SS1) //trim(IntToStr(N2))
    if (present(SS3)) then
        concat = trim(concat) // SS3
        if (present(N4)) then
            concat = trim(concat) // trim(IntToStr(N4))
            if (present(SS5)) then
                concat = trim(concat) // SS5
                if (present(N6)) then
                    concat = trim(concat) // trim(intToStr(N6))
                    if (present(SS7)) then
                        concat = trim(concat) // SS7
                        if (present(N8)) then
                            concat = trim(concat) // trim(intToStr(N8))
                            if (present(SS9)) then
                                concat = trim(concat) // SS9
                                if (present(N10)) then
                                    concat = trim(concat) // trim(intToStr(N10))
                                    if (present(SS11)) then
                                        concat = trim(concat) // SS11
                                    end if
                                end if
                            end if
                        end if
                    end if
                end if
            end if
        end if
    end if

    end  function concat_s_n

    subroutine Exchange(i1,i2)
    integer i1,i2,tmp

    tmp=i1
    i1=i2
    i2=tmp

    end subroutine Exchange

    subroutine WriteS(S)
    character(LEN=*), intent(in) :: S

    write (*,*) trim(S)

    end subroutine WriteS

    subroutine StringReplace(FindS, RepS, S)
    character(LEN=*), intent(in) :: FindS, RepS
    character(LEN=*), intent(inout) :: S
    integer i

    i = index(S,FindS)
    if (i>0) then
        S = S(1:i-1)//trim(RepS)//S(i+len_trim(FindS):len_trim(S))
    end if

    end subroutine StringReplace

    function numcat(S, num)
    character(LEN=*) S
    character(LEN=1024) numcat, numstr
    integer num

    write (numstr, *) num
    numcat = trim(S) // trim(adjustl(numstr))
    !OK, so can probably do with with a format statement too...
    end function numcat

    function LogicalToint(B)
    integer LogicalToint
    logical, intent(in) :: B

    if (B) then
        LogicalToInt=1
    else
        LogicalToint=0
    end if

    end function LogicalToInt

    function IntToLogical(I)
    integer, intent(in) :: I
    logical IntToLogical

    IntToLogical = I /= 0

    end  function IntToLogical

    function IntToStr(I, minlen)
    integer , intent(in) :: I
    character(LEN=30) IntToStr
    integer, intent(in), optional :: minlen
    integer n
    character (LEN=20) :: form

    if (present(minlen)) then
        n = minlen
        if (I<0) n=n+1
        form = concat('(I',n,'.',minlen,')')
        write (IntToStr,form) i
    else
        write (IntToStr,*) i
        IntToStr = adjustl(IntToStr)
    end if


    end function IntToStr

    function StrToInt(S)
    integer :: StrToInt
    character(LEN=30), intent(in) :: S

    read (S,*) StrToInt
    end function StrToInt

    function DoubleToStr(R, figs)
    double precision, intent(in) :: R
    integer, intent(in), optional :: figs
    character(LEN=30) DoubleToStr

    DoubleToStr = SingleToStr(real(R),figs)

    end function DoubleToStr

    function SingleToStr(R, figs)
    real, intent(in) :: R
    integer, intent(in), optional :: figs
    character(LEN=30) SingleToStr

    if (abs(R)>=0.001 .or. R==0.) then
        write (SingleToStr,'(f12.6)') R

        SingleToStr = adjustl(SingleToStr)
        if (present(figs)) then
            SingleToStr = SingleToStr(1:figs)
        else
            SingleToStr = SingleToStr(1:6)
        end if

    else
        if (present(figs)) then
            write (SingleToStr,trim(numcat('(E',figs))//'.2)') R
        else
            write (SingleToStr,'(G9.2)') R
        end if
        SingleToStr = adjustl(SingleToStr)
    end if

    end function SingleToStr

    subroutine WriteFormatInts(unit, formatst, i1,i2,i3,i4)
    integer, intent(in) :: unit
    character(LEN=*), intent(in) :: formatst
    integer, intent(in) :: i1
    integer, intent(in),optional :: i2,i3,i4
    character(LEN=1024*16) S

    S = formatst
    call StringReplace('%u', IntToStr(i1), S)
    if (present(i2)) call StringReplace('%u', IntToStr(i2), S)
    if (present(i3)) call StringReplace('%u', IntToStr(i3), S)
    if (present(i4)) call StringReplace('%u', IntToStr(i4), S)

    write(unit,'(a)') trim(S)
    end subroutine WriteFormatInts

    function IndexOf(aval,arr, n)
    integer, intent(in) :: n, arr(n), aval
    integer IndexOf, i

    do i=1,n
        if (arr(i)==aval) then
            IndexOf= i
            return
        end if
    end do
    IndexOf = 0

    end function IndexOf

    function MaxIndex(arr, n)
    integer, intent(in) :: n
    real, intent(in) :: arr(n)
    integer locs(1:1), MaxIndex

    locs = maxloc(arr(1:n))
    MaxIndex = locs(1)

    end function MaxIndex


    function MinIndex(arr, n)
    integer, intent(in) :: n
    real, intent(in) :: arr(n)
    integer locs(1:1), MinIndex

    locs = minloc(arr(1:n))
    MinIndex = locs(1)

    end function MinIndex


    subroutine TList_RealArr_SaveToFile(L,fname)
    character(LEN=*), intent(IN) :: fname
    Type (TList_RealArr) :: L
    character(LEN=20) aform
    integer i
    integer :: Plen = -1
    integer :: file_id

    file_id = new_file_unit()
    call CreateTxtFile(fname,file_id)
    do i=1, L%Count
        if (PLen /= size(L%Items(i)%P)) then
            PLen = size(L%Items(i)%P)
            aform = '('//trim(IntToStr(PLen))//'E16.8)'
        end if
        write (file_id,aform) L%Items(i)%P
    end do
    call CloseFile(file_id)

    end subroutine TList_RealArr_SaveToFile


    function ExtractFilePath(aname)
    character(LEN=*), intent(IN) :: aname
    character(LEN=1024) ExtractFilePath
    integer len, i

    len = len_trim(aname)
    do i = len, 1, -1
        if (aname(i:i)=='/') then
            ExtractFilePath = aname(1:i)
            return
        end if
    end do
    ExtractFilePath = ''

    end function ExtractFilePath

    function ExtractFileExt(aname)
    character(LEN=*), intent(IN) :: aname
    character(LEN=120) ExtractFileExt
    integer len, i

    len = len_trim(aname)
    do i = len, 1, -1
        if (aname(i:i)=='/') then
            ExtractFileExt = ''
            return
        else if (aname(i:i)=='.') then
            ExtractFileExt= aname(i:len)
            return
        end if
    end do
    ExtractFileExt = ''

    end function ExtractFileExt


    function ExtractFileName(aname)
    character(LEN=*), intent(IN) :: aname
    character(LEN=1024) ExtractFileName
    integer len, i

    len = len_trim(aname)
    do i = len, 1, -1
        if (aname(i:i)=='/') then
            ExtractFileName = aname(i+1:len)
            return
        end if
    end do
    ExtractFileName = aname

    end function ExtractFileName

    function ChangeFileExt(aname,ext)
    character(LEN=*), intent(IN) :: aname,ext
    character(LEN=1024) ChangeFileExt
    integer len, i

    len = len_trim(aname)
    do i = len, 1, -1
        if (aname(i:i)=='.') then
            ChangeFileExt = aname(1:i) // trim(ext)
            return
        end if
    end do
    ChangeFileExt = trim(aname) // '.' // trim(ext)

    end function ChangeFileExt


    function CheckTrailingSlash(aname)
    character(LEN=*), intent(in) :: aname
    character(LEN=1024) CheckTrailingSlash
    integer len

    len = len_trim(aname)
#ifdef IBMXL
    if (aname(len:len) /= '\\' .and. aname(len:len) /= '/') then
#else
#ifdef ESCAPEBACKSLASH
        if (aname(len:len) /= '\\' .and. aname(len:len) /= '/') then
#else
            if (aname(len:len) /= '\' .and. aname(len:len) /= '/') then
#endif
#endif
                CheckTrailingSlash = trim(aname)//'/'
            else
                CheckTrailingSlash = aname
            end if


        end  function CheckTrailingSlash


        subroutine DeleteFile(aname)
        character(LEN=*), intent(IN) :: aname
        integer file_id

        if (FileExists(aname)) then
            file_id = new_file_unit()
            open(unit = file_id, file = aname, err = 2)
            close(unit = file_id, status = 'DELETE')
2           file_units(file_id) = .false.
        end if

        end subroutine DeleteFile


        subroutine FlushFile(aunit)
#ifdef __INTEL_COMPILER_BUILD_DATE
        USE IFPORT
#endif
        integer, intent(IN) :: aunit


#ifdef IBMXL
        call flush_(aunit)
#else
        call flush(aunit)
#endif

        end subroutine FlushFile


        function FileExists(aname)
        character(LEN=*), intent(IN) :: aname
        logical FileExists

        inquire(file=aname, exist = FileExists)

        end function FileExists

        subroutine OpenFile(aname, aunit,mode)
        character(LEN=*), intent(IN) :: aname,mode
        integer, intent(in) :: aunit

        open(unit=aunit,file=aname,form=mode,status='old', action='read', err=500)
        return

500     call MpiStop('File not found: '//trim(aname))


        end subroutine OpenFile


        subroutine OpenTxtFile(aname, aunit)
        character(LEN=*), intent(IN) :: aname
        integer, intent(in) :: aunit

        call OpenFile(aname,aunit,'formatted')

        end subroutine OpenTxtFile

        subroutine CreateOpenTxtFile(aname, aunit, append)
        character(LEN=*), intent(IN) :: aname
        integer, intent(in) :: aunit
        logical, optional, intent(in) :: append
        logical A

        if (present(append)) then
            A=append
        else
            A = .false.
        endif

        call CreateOpenFile(aname,aunit,'formatted',A)

        end subroutine CreateOpenTxtFile


        subroutine CreateTxtFile(aname, aunit)
        character(LEN=*), intent(IN) :: aname
        integer, intent(in) :: aunit

        call CreateFile(aname,aunit,'formatted')

        end subroutine CreateTxtFile


        subroutine CreateFile(aname, aunit,mode)
        character(LEN=*), intent(IN) :: aname,mode
        integer, intent(in) :: aunit

        open(unit=aunit,file=aname,form=mode,status='replace', err=500)

        return

500     call MpiStop('Error creating file '//trim(aname))


        end subroutine CreateFile

        subroutine CreateOpenFile(aname, aunit,mode, append)
        character(LEN=*), intent(IN) :: aname,mode
        integer, intent(in) :: aunit
        logical, optional, intent(in) :: append
        logical A

        if (present(append)) then
            A=append
        else
            A = .false.
        endif

        if (A) then
            open(unit=aunit,file=aname,form=mode,status='unknown', err=500, position='append')
        else
            open(unit=aunit,file=aname,form=mode,status='replace', err=500)
        end if

        return

500     call MpiStop('Error creatinging or opening '//trim(aname))


        end subroutine CreateOpenFile

        function TxtNumberColumns(InLine) result(n)
        character(LEN=*) :: InLine
        integer n,i
        logical isNum

        n=0
        isNum=.false.
        do i=1, len_trim(InLIne)
            if (verify(InLine(i:i),'-+eE.0123456789') == 0) then
                if (.not. IsNum) n=n+1
                IsNum=.true.
            else
                IsNum=.false.
            end if
        end do

        end function TxtNumberColumns

        function TxtColumns(InLine) result(n)
        character(LEN=*) :: InLine
        integer n,i
        logical isNum

        n=0
        isNum=.false.
        do i=1, len_trim(InLine)
            if (InLine(i:i) > char(32)) then
                if (.not. IsNum) n=n+1
                IsNum=.true.
            else
                IsNum=.false.
            end if
        end do

        end function TxtColumns

        function FileColumns(aunit) result(n)
        integer, intent(in) :: aunit
        integer n
        character(LEN=4096*32) :: InLine

        n=0
        read(aunit,'(a)', end = 10) InLine
        n = TxtNumberColumns(InLine)
10      rewind aunit

        end function FileColumns

        function FileLines(aunit) result(n)
        integer, intent(in) :: aunit
        integer n
        character(LEN=4096) :: InLine

        n=0
        do

            read(aunit,'(a)', end = 200) InLine
            n = n+1
        end do

200     rewind aunit


        end function FileLines


        function TopCommentLine(aname) result(res)
        character(LEN=*), intent(IN) :: aname
        integer file_id
        character(LEN=1024) :: InLine, res

        res = ''
        file_id = new_file_unit()
        call OpenTxtFile(aname, file_id)
        InLine=''
        do while (InLine /= '')
            read(file_id,'(a)', end = 10) InLine
        end do
        If (InLIne(1:1)=='#') then
            res = InLine
        end if

10      call CloseFile(file_id)

        end function TopCommentLine


        function TxtFileColumns(aname) result(n)
        character(LEN=*), intent(IN) :: aname
        integer n, file_id


        file_id = new_file_unit()

        call OpenTxtFile(aname, file_id)
        n = FileColumns(file_id)
        call CloseFile(file_id)

        end function TxtFileColumns


        function LastFileLine(aname)
        character(LEN=*), intent(IN) :: aname
        character(LEN = 5000) LastFileLine, InLine
        integer  file_id

        file_id = new_file_unit()

        InLine = ''
        call OpenTxtFile(aname,file_id)
        do
            read(file_id,'(a)', end = 200) InLine
        end do

200     call CloseFile(file_id)

        LastFileLine = InLine

        end function LastFileLine

        subroutine writeArrayLine(unit, arr)
        real, intent(in) :: arr(:)
        integer, intent(in) :: unit

        write(unit,concat('(',size(arr),'E16.6)')) arr

        end subroutine writeArrayLine


        subroutine spline_real(x,y,n,d2)
        integer, intent(in) :: n
        integer, parameter :: dp=KIND(1.0)
        real(dp), intent(in) :: x(n), y(n)
        real(dp), intent(out) :: d2(n)
        real(dp), dimension(:), allocatable :: u
        integer i
        real(dp) xp,sig,xxdiv,d1l,d1r

        allocate(u(1:n-1))

        d2(1)=0._dp
        u(1)=0._dp

        d1r= (y(2)-y(1))/(x(2)-x(1))
        do i=2,n-1
            d1l=d1r
            d1r=(y(i+1)-y(i))/(x(i+1)-x(i))
            xxdiv=1._dp/(x(i+1)-x(i-1))
            sig=(x(i)-x(i-1))*xxdiv
            xp=1._dp/(sig*d2(i-1)+2._dp)
            d2(i)=(sig-1._dp)*xp
            u(i)=(6._dp*(d1r-d1l)*xxdiv-sig*u(i-1))*xp
        end do

        d2(n)=0._dp
        do i=n-1,1,-1
            d2(i)=d2(i)*d2(i+1)+u(i)
        end do

        deallocate(u)
        end subroutine spline_real

        subroutine spline_double(x,y,n,d2)
        integer, intent(in) :: n
        integer, parameter :: dp=KIND(1.d0)
        real(dp), intent(in) :: x(n), y(n)
        real(dp), intent(out) :: d2(n)
        real(dp), dimension(:), allocatable :: u
        integer i
        real(dp) xp,sig,xxdiv,d1l,d1r

        allocate(u(1:n-1))

        d2(1)=0._dp
        u(1)=0._dp

        d1r= (y(2)-y(1))/(x(2)-x(1))
        do i=2,n-1
            d1l=d1r
            d1r=(y(i+1)-y(i))/(x(i+1)-x(i))
            xxdiv=1._dp/(x(i+1)-x(i-1))
            sig=(x(i)-x(i-1))*xxdiv
            xp=1._dp/(sig*d2(i-1)+2._dp)
            d2(i)=(sig-1._dp)*xp
            u(i)=(6._dp*(d1r-d1l)*xxdiv-sig*u(i-1))*xp
        end do

        d2(n)=0._dp
        do i=n-1,1,-1
            d2(i)=d2(i)*d2(i+1)+u(i)
        end do

        deallocate(u)
        end subroutine spline_double


        function DLGAMMA(x)
        !Use Stirling generalization for large x
        !See e.g. http://en.wikipedia.org/wiki/Stirling's_approximation
        !Is accurate to at least 10 decimals, worse just about 30
        double precision :: x
        double precision:: DLGAMMA !approx log gamma
        double precision, parameter :: const = .91893853320467274180d0 !log(2pi)/2

        if (x<32.d0) then
            DLGAMMA = log(GAMMA(x))
        else
            DLGAMMA = (x-0.5d0)*log(x) - x + const +  &
                1/12.d0/(1+x)*(1+1/(x+2)*(1+59.d0/30/(x+3)*(1+2.9491525423728813559d0/(x+4))))
        end if
        end function DLGAMMA


        function LogGamma(x)
        real LogGamma
        real, intent(in) :: x
        integer i, j
        real r

        i = nint(x*2)
        if (abs(i-x*2) > 1e-4) call MpiStop('LogGamma function for half integral only')
        if (mod(i,2) == 0) then
            r=0
            do j = 2, i/2-1
                r = r + log(real(j))
            end do
            LogGamma = r
        else
            r = log(pi)/2
            do j = 1, i-2 , 2
                r = r+ log(j/2.0)
            end do
            LogGamma = r
        end if

        end function LogGamma

        DOUBLE PRECISION FUNCTION GAMMA(X)
        !----------------------------------------------------------------------
        !
        ! This routine calculates the GAMMA function for a real argument X.
        !   Computation is based on an algorithm outlined in reference 1.
        !   The program uses rational functions that approximate the GAMMA
        !   function to at least 20 significant decimal digits.  Coefficients
        !   for the approximation over the interval (1,2) are unpublished.
        !   Those for the approximation for X .GE. 12 are from reference 2.
        !   The accuracy achieved depends on the arithmetic system, the
        !   compiler, the intrinsic functions, and proper selection of the
        !   machine-dependent constants.
        !*******************************************************************
        !
        ! Explanation of machine-dependent constants
        !
        ! beta   - radix for the floating-point representation
        ! maxexp - the smallest positive power of beta that overflows
        ! XBIG   - the largest argument for which GAMMA(X) is representable
        !          in the machine, i.e., the solution to the equation
        !                  GAMMA(XBIG) = beta**maxexp
        ! XINF   - the largest machine representable floating-point number;
        !          approximately beta**maxexp
        ! EPS    - the smallest positive floating-point number such that
        !          1.0+EPS .GT. 1.0
        ! XMININ - the smallest positive floating-point number such that
        !          1/XMININ is machine representable
        !
        !     Approximate values for some important machines are:
        !
        !                            beta       maxexp        XBIG
        !
        ! CRAY-1         (S.P.)        2         8191        966.961
        ! Cyber 180/855
        !   under NOS    (S.P.)        2         1070        177.803
        ! IEEE (IBM/XT,
        !   SUN, etc.)   (S.P.)        2          128        35.040
        ! IEEE (IBM/XT,
        !   SUN, etc.)   (D.P.)        2         1024        171.624
        ! IBM 3033       (D.P.)       16           63        57.574
        ! VAX D-Format   (D.P.)        2          127        34.844
        ! VAX G-Format   (D.P.)        2         1023        171.489
        !
        !                            XINF         EPS        XMININ
        !
        ! CRAY-1         (S.P.)   5.45E+2465   7.11E-15    1.84E-2466
        ! Cyber 180/855
        !   under NOS    (S.P.)   1.26E+322    3.55E-15    3.14E-294
        ! IEEE (IBM/XT,
        !   SUN, etc.)   (S.P.)   3.40E+38     1.19E-7     1.18E-38
        ! IEEE (IBM/XT,
        !   SUN, etc.)   (D.P.)   1.79D+308    2.22D-16    2.23D-308
        ! IBM 3033       (D.P.)   7.23D+75     2.22D-16    1.39D-76
        ! VAX D-Format   (D.P.)   1.70D+38     1.39D-17    5.88D-39
        ! VAX G-Format   (D.P.)   8.98D+307    1.11D-16    1.12D-308
        !
        !*******************************************************************
        !*******************************************************************
        !
        ! Error returns
        !
        !  The program returns the value XINF for singularities or
        !     when overflow would occur.  The computation is believed
        !     to be free of underflow and overflow.
        !
        !
        !  Intrinsic functions required are:
        !
        !     INT, DBLE, EXP, LOG, REAL, SIN
        !
        !
        ! References: "An Overview of Software Development for Special
        !              Functions", W. J. Cody, Lecture Notes in Mathemati,
        !              506, Numerical Analysis Dundee, 1975, G. A. Watson
        !              (ed.), Springer Verlag, Berlin, 1976.
        !
        !              Computer Approximations, Hart, Et. Al., Wiley and
        !              sons, New York, 1968.
        !
        !  Latest modification: October 12, 1989
        !
        !  Authors: W. J. Cody and L. Stoltz
        !           Applied Mathemati Division
        !           Argonne National Laboratory
        !           Argonne, IL 60439
        !
        !----------------------------------------------------------------------
        INTEGER I,N
        LOGICAL PARITY
        DOUBLE PRECISION C,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI,SUM,TWELVE, &
            TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO
        DIMENSION C(7),P(8),Q(8)
        !----------------------------------------------------------------------
        !  Mathematical constants
        !----------------------------------------------------------------------
        DATA ONE,HALF,TWELVE,TWO,ZERO/1.0D0,0.5D0,12.0D0,2.0D0,0.0D0/, &
            SQRTPI/0.9189385332046727417803297D0/, &
            PI/3.1415926535897932384626434D0/
        !----------------------------------------------------------------------
        !  Machine dependent parameters
        !----------------------------------------------------------------------
        DATA XBIG,XMININ,EPS/35.040D0,1.18D-38,1.19D-7/, &
            XINF/3.4E38/
        !----------------------------------------------------------------------
        !  Numerator and denominator coefficients for rational minimax
        !     approximation over (1,2).
        !----------------------------------------------------------------------
        DATA P/-1.71618513886549492533811E+0,2.47656508055759199108314E+1, &
            -3.79804256470945635097577E+2,6.29331155312818442661052E+2, &
            8.66966202790413211295064E+2,-3.14512729688483675254357E+4, &
            -3.61444134186911729807069E+4,6.64561438202405440627855E+4/
        DATA Q/-3.08402300119738975254353E+1,3.15350626979604161529144E+2, &
            -1.01515636749021914166146E+3,-3.10777167157231109440444E+3, &
            2.25381184209801510330112E+4,4.75584627752788110767815E+3, &
            -1.34659959864969306392456E+5,-1.15132259675553483497211E+5/
        !----------------------------------------------------------------------
        ! Coefficients for minimax approximation over (12, INF).
        !----------------------------------------------------------------------
        DATA C/-1.910444077728D-03,8.4171387781295D-04, &
            -5.952379913043012D-04,7.93650793500350248D-04, &
            -2.777777777777681622553D-03,8.333333333333333331554247D-02, &
            5.7083835261D-03/
        !----------------------------------------------------------------------
        !  Statement functions for conversion between integer and float
        !----------------------------------------------------------------------
        PARITY = .FALSE.
        FACT = ONE
        N = 0
        Y = X
        IF (Y .LE. ZERO) THEN
            !----------------------------------------------------------------------
            !  Argument is negative
            !----------------------------------------------------------------------
            Y = -X
            Y1 = AINT(Y)
            RES = Y - Y1
            IF (RES .NE. ZERO) THEN
                IF (Y1 .NE. AINT(Y1*HALF)*TWO) PARITY = .TRUE.
                FACT = -PI / SIN(PI*RES)
                Y = Y + ONE
            ELSE
                RES = XINF
                GO TO 900
            END IF
        END IF
        !----------------------------------------------------------------------
        !  Argument is positive
        !----------------------------------------------------------------------
        IF (Y .LT. EPS) THEN
            !----------------------------------------------------------------------
            !  Argument .LT. EPS
            !----------------------------------------------------------------------
            IF (Y .GE. XMININ) THEN
                RES = ONE / Y
            ELSE
                RES = XINF
                GO TO 900
            END IF
        ELSE IF (Y .LT. TWELVE) THEN
            Y1 = Y
            IF (Y .LT. ONE) THEN
                !----------------------------------------------------------------------
                !  0.0 .LT. argument .LT. 1.0
                !----------------------------------------------------------------------
                Z = Y
                Y = Y + ONE
            ELSE
                !----------------------------------------------------------------------
                !  1.0 .LT. argument .LT. 12.0, reduce argument if necessary
                !----------------------------------------------------------------------
                N = INT(Y) - 1
                Y = Y - REAL(N)
                Z = Y - ONE
            END IF
            !----------------------------------------------------------------------
            !  Evaluate approximation for 1.0 .LT. argument .LT. 2.0
            !----------------------------------------------------------------------
            XNUM = ZERO
            XDEN = ONE
            DO 260 I = 1, 8
                XNUM = (XNUM + P(I)) * Z
                XDEN = XDEN * Z + Q(I)
260         CONTINUE
            RES = XNUM / XDEN + ONE
            IF (Y1 .LT. Y) THEN
                !----------------------------------------------------------------------
                !  Adjust result for case  0.0 .LT. argument .LT. 1.0
                !----------------------------------------------------------------------
                RES = RES / Y1
            ELSE IF (Y1 .GT. Y) THEN
                !----------------------------------------------------------------------
                !  Adjust result for case  2.0 .LT. argument .LT. 12.0
                !----------------------------------------------------------------------
                DO 290 I = 1, N
                    RES = RES * Y
                    Y = Y + ONE
290             CONTINUE
            END IF
        ELSE
            !----------------------------------------------------------------------
            !  Evaluate for argument .GE. 12.0,
            !----------------------------------------------------------------------
            IF (Y .LE. XBIG) THEN
                YSQ = Y * Y
                SUM = C(7)
                DO 350 I = 1, 6
                    SUM = SUM / YSQ + C(I)
350             CONTINUE
                SUM = SUM/Y - Y + SQRTPI
                SUM = SUM + (Y-HALF)*LOG(Y)
                RES = EXP(SUM)
            ELSE
                RES = XINF
                GO TO 900
            END IF
        END IF
        !----------------------------------------------------------------------
        !  Final adjustments and return
        !----------------------------------------------------------------------
        IF (PARITY) RES = -RES
        IF (FACT .NE. ONE) RES = FACT / RES
900     GAMMA = RES

        END FUNCTION GAMMA

        subroutine SetIdlePriority
#ifdef RUNIDLE
        USE DFWIN
        Integer dwPriority
        Integer CheckPriority

        dwPriority = 64 ! idle priority
        CheckPriority = SetPriorityClass(GetCurrentProcess(), dwPriority)
#endif
        end subroutine SetIdlePriority


        subroutine GetThreeJs(thrcof,l2in,l3in,m2in,m3in)
        !Recursive evaluation of 3j symbols. Does minimal error checking on input parameters.
        implicit none
        integer, parameter :: dl = KIND(1.d0)
        integer, intent(in) :: l2in,l3in, m2in,m3in
        real(dl), dimension(*) :: thrcof
        INTEGER, PARAMETER :: i8 = selected_int_kind(18)
        integer(i8) :: l2,l3,m2,m3
        integer(i8) :: l1, m1, l1min,l1max, lmatch, nfin, a1, a2

        real(dl) :: newfac, oldfac, sumfor, c1,c2,c1old, dv, denom, x, sum1, sumuni
        real(dl) :: x1,x2,x3, y,y1,y2,y3,sum2,sumbac, ratio,cnorm, sign1, thresh
        integer i,ier, index, nlim, sign2
        integer nfinp1,nfinp2,nfinp3, lstep, nstep2,n
        real(dl), parameter :: zero = 0._dl, one = 1._dl
        real(dl), parameter ::  tiny = 1.0d-30, srtiny=1.0d-15, huge = 1.d30, srhuge = 1.d15

        ! routine to generate set of 3j-coeffs (l1,l2,l3\\ m1,m2,m3)

        ! by recursion from l1min = max(abs(l2-l3),abs(m1))
        !                to l1max = l2+l3
        ! the resulting 3j-coeffs are stored as thrcof(l1-l1min+1)

        ! to achieve the numerical stability, the recursion will proceed
        ! simultaneously forwards and backwards, starting from l1min and l1max
        ! respectively.
        !
        ! lmatch is the l1-value at which forward and backward recursion are matched.
        !
        ! ndim is the length of the array thrcof
        !
        ! ier = -1 for all 3j vanish(l2-abs(m2)<0, l3-abs(m3)<0 or not integer)
        ! ier = -2 if possible 3j's exceed ndim
        ! ier >= 0 otherwise

        l2=l2in
        l3=l3in
        m2=m2in
        m3=m3in
        newfac = 0
        lmatch = 0
        m1 = -(m2+m3)

        ! check relative magnitude of l and m values
        ier = 0

        if (l2 < abs(m2) .or. l3 < m3) then
            ier = -1
            call MpiStop('error ier = -1')
            return
        end if

        ! limits for l1
        l1min = max(abs(l2-l3),abs(m1))
        l1max = l2+l3

        if (l1min >= l1max) then
            if (l1min/=l1max) then
                ier = -1
                call MpiStop('error ier = -1')
                return
            end if

            ! reached if l1 can take only one value, i.e.l1min=l1max
            thrcof(1) = (-1)**abs(l2+m2-l3+m3)/sqrt(real(l1min+l2+l3+1,dl))
            return

        end if

        nfin = l1max-l1min+1

        ! starting forward recursion from l1min taking nstep1 steps
        l1 = l1min
        thrcof(1) = srtiny
        sum1 = (2*l1 + 1)*tiny

        lstep = 1

30      lstep = lstep+1
        l1 = l1+1

        oldfac = newfac
        a1 = (l1+l2+l3+1)*(l1-l2+l3)*(l1+l2-l3)
        a2 = (l1+m1)*(l1-m1)*(-l1+l2+l3+1)
        newfac = sqrt(a2*real(a1,dl))
        if (l1 == 1) then
            !IF L1 = 1  (L1-1) HAS TO BE FACTORED OUT OF DV, HENCE
            c1 = -(2*l1-1)*l1*(m3-m2)/newfac
        else

            dv = -l2*(l2+1)*m1 + l3*(l3+1)*m1 + l1*(l1-1)*(m3-m2)
            denom = (l1-1)*newfac

            if (lstep > 2) c1old = abs(c1)
            c1 = -(2*l1-1)*dv/denom

        end if

        if (lstep<= 2) then

            ! if l1=l1min+1 the third term in the recursion eqn vanishes, hence
            x = srtiny*c1
            thrcof(2) = x
            sum1 = sum1+tiny*(2*l1+1)*c1*c1
            if(lstep==nfin) then
                sumuni=sum1
                go to 230
            end if
            goto 30

        end if

        c2 = -l1*oldfac/denom

        ! recursion to the next 3j-coeff x
        x = c1*thrcof(lstep-1) + c2*thrcof(lstep-2)
        thrcof(lstep) = x
        sumfor = sum1
        sum1 = sum1 + (2*l1+1)*x*x
        if (lstep/=nfin) then

            ! see if last unnormalised 3j-coeff exceeds srhuge
            if (abs(x) >= srhuge) then

                ! REACHED IF LAST 3J-COEFFICIENT LARGER THAN SRHUGE
                ! SO THAT THE RECURSION SERIES THRCOF(1), ... , THRCOF(LSTEP)
                ! HAS TO BE RESCALED TO PREVENT OVERFLOW

                ier = ier+1
                do i = 1, lstep
                    if (abs(thrcof(i)) < srtiny) thrcof(i)= zero
                    thrcof(i) = thrcof(i)/srhuge
                end do

                sum1 = sum1/huge
                sumfor = sumfor/huge
                x = x/srhuge

            end if

            ! as long as abs(c1) is decreasing, the recursion proceeds towards increasing
            ! 3j-valuse and so is numerically stable. Once an increase of abs(c1) is
            ! detected, the recursion direction is reversed.

            if (c1old > abs(c1)) goto 30

        end if !lstep/=nfin

        ! keep three 3j-coeffs around lmatch for comparison with backward recursion

        lmatch = l1-1
        x1 = x
        x2 = thrcof(lstep-1)
        x3 = thrcof(lstep-2)
        nstep2 = nfin-lstep+3

        ! --------------------------------------------------------------------------
        !
        ! starting backward recursion from l1max taking nstep2 stpes, so that
        ! forward and backward recursion overlap at 3 points
        ! l1 = lmatch-1, lmatch, lmatch+1

        nfinp1 = nfin+1
        nfinp2 = nfin+2
        nfinp3 = nfin+3
        l1 = l1max
        thrcof(nfin) = srtiny
        sum2 = tiny*(2*l1+1)

        l1 = l1+2
        lstep=1

        do
            lstep = lstep + 1
            l1= l1-1

            oldfac = newfac
            a1 = (l1+l2+l3)*(l1-l2+l3-1)*(l1+l2-l3-1)
            a2 = (l1+m1-1)*(l1-m1-1)*(-l1+l2+l3+2)
            newfac = sqrt(a1*real(a2,dl))

            dv = -l2*(l2+1)*m1 + l3*(l3+1)*m1 +l1*(l1-1)*(m3-m2)

            denom = l1*newfac
            c1 = -(2*l1-1)*dv/denom
            if (lstep <= 2) then

                ! if l2=l2max+1, the third term in the recursion vanishes

                y = srtiny*c1
                thrcof(nfin-1) = y
                sumbac = sum2
                sum2 = sum2 + tiny*(2*l1-3)*c1*c1

                cycle

            end if

            c2 = -(l1-1)*oldfac/denom

            ! recursion to the next 3j-coeff y
            y = c1*thrcof(nfinp2-lstep)+c2*thrcof(nfinp3-lstep)

            if (lstep==nstep2) exit

            thrcof(nfinp1-lstep) = y
            sumbac = sum2
            sum2 = sum2+(2*l1-3)*y*y

            ! see if last unnormalised 3j-coeff exceeds srhuge
            if (abs(y) >= srhuge) then

                ! reached if 3j-coeff larger than srhuge so that the recursion series
                ! thrcof(nfin),..., thrcof(nfin-lstep+1) has to be rescaled to prevent overflow

                ier=ier+1
                do i = 1, lstep
                    index=nfin-i+1
                    if (abs(thrcof(index)) < srtiny) thrcof(index)=zero
                    thrcof(index) = thrcof(index)/srhuge
                end do

                sum2=sum2/huge
                sumbac=sumbac/huge

            end if

        end do

        ! the forward recursion 3j-coeffs x1, x2, x3 are to be matched with the
        ! corresponding backward recursion vals y1, y2, y3

        y3 = y
        y2 = thrcof(nfinp2-lstep)
        y1 = thrcof(nfinp3-lstep)

        ! determine now ratio such that yi=ratio*xi (i=1,2,3) holds with minimal error

        ratio = (x1*y1+x2*y2+x3*y3)/(x1*x1+x2*x2+x3*x3)
        nlim = nfin-nstep2+1

        if (abs(ratio) >= 1) then

            thrcof(1:nlim) = ratio*thrcof(1:nlim)
            sumuni = ratio*ratio*sumfor + sumbac

        else

            nlim = nlim+1
            ratio = 1/ratio
            do n = nlim, nfin
                thrcof(n) = ratio*thrcof(n)
            end do
            sumuni = sumfor + ratio*ratio*sumbac

        end if
        ! normalise 3j-coeffs

230     cnorm = 1/sqrt(sumuni)

        ! sign convention for last 3j-coeff determines overall phase

        sign1 = sign(one,thrcof(nfin))
        sign2 = (-1)**(abs(l2+m2-l3+m3))
        if (sign1*sign2 <= 0) then
            cnorm = -cnorm
        end if
        if (abs(cnorm) >= one) then
            thrcof(1:nfin) = cnorm*thrcof(1:nfin)
            return
        end if

        thresh = tiny/abs(cnorm)

        do n = 1, nfin
            if (abs(thrcof(n)) < thresh) thrcof(n) = zero
            thrcof(n) = cnorm*thrcof(n)
        end do
        return

        end subroutine GetThreeJs



    end module AMLutils


#ifdef ZIGGURAT
    MODULE Ziggurat
    ! Marsaglia & Tsang generator for random normals & random exponentials.
    ! Translated from C by Alan Miller (amiller@bigpond.net.au)

    ! Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating
    ! random variables', J. Statist. Software, v5(8).

    ! This is an electronic journal which can be downloaded from:
    ! http://www.jstatsoft.org/v05/i08

    ! N.B. It is assumed that all integers are 32-bit.
    ! N.B. The value of M2 has been halved to compensate for the lack of
    !      unsigned integers in Fortran.

    ! Latest version - 1 January 2001
    !
    ! AL: useful material at http://en.wikipedia.org/wiki/Ziggurat_algorithm
    IMPLICIT NONE

    PRIVATE

    INTEGER,  PARAMETER  ::  DP=SELECTED_REAL_KIND( 12, 60 )
    REAL(DP), PARAMETER  ::  m1=2147483648.0_DP,   m2=2147483648.0_DP,      &
        half=0.5_DP
    REAL(DP)             ::  dn=3.442619855899_DP, tn=3.442619855899_DP,    &
        vn=0.00991256303526217_DP,                     &
        q,                    de=7.697117470131487_DP, &
        te=7.697117470131487_DP,                       &
        ve=0.003949659822581572_DP
    INTEGER,  SAVE       ::  iz, jz, jsr=123456789, kn(0:127),              &
        ke(0:255), hz
    REAL(DP), SAVE       ::  wn(0:127), fn(0:127), we(0:255), fe(0:255)
    LOGICAL,  SAVE       ::  initialized=.FALSE.

    PUBLIC  :: zigset, shr3, uni, rnor, rexp


    CONTAINS


    SUBROUTINE zigset( jsrseed )

    INTEGER, INTENT(IN)  :: jsrseed

    INTEGER  :: i

    !  Set the seed
    jsr = jsrseed

    !  Tables for RNOR
    q = vn*EXP(half*dn*dn)
    kn(0) = (dn/q)*m1
    kn(1) = 0
    wn(0) = q/m1
    wn(127) = dn/m1
    fn(0) = 1.0_DP
    fn(127) = EXP( -half*dn*dn )
    DO  i = 126, 1, -1
        dn = SQRT( -2.0_DP * LOG( vn/dn + EXP( -half*dn*dn ) ) )
        kn(i+1) = (dn/tn)*m1
        tn = dn
        fn(i) = EXP(-half*dn*dn)
        wn(i) = dn/m1
    END DO

    !  Tables for REXP
    q = ve*EXP( de )
    ke(0) = (de/q)*m2
    ke(1) = 0
    we(0) = q/m2
    we(255) = de/m2
    fe(0) = 1.0_DP
    fe(255) = EXP( -de )
    DO  i = 254, 1, -1
        de = -LOG( ve/de + EXP( -de ) )
        ke(i+1) = m2 * (de/te)
        te = de
        fe(i) = EXP( -de )
        we(i) = de/m2
    END DO
    initialized = .TRUE.
    RETURN
    END SUBROUTINE zigset



    !  Generate random 32-bit integers
    FUNCTION shr3( ) RESULT( ival )
    INTEGER  ::  ival

    jz = jsr
    jsr = IEOR( jsr, ISHFT( jsr,  13 ) )
    jsr = IEOR( jsr, ISHFT( jsr, -17 ) )
    jsr = IEOR( jsr, ISHFT( jsr,   5 ) )
    ival = jz + jsr
    RETURN
    END FUNCTION shr3



    !  Generate uniformly distributed random numbers
    FUNCTION uni( ) RESULT( fn_val )
    REAL(DP)  ::  fn_val

    fn_val = half + 0.2328306e-9_DP * shr3( )
    RETURN
    END FUNCTION uni



    !  Generate random normals
    FUNCTION rnor( ) RESULT( fn_val )
    REAL(DP)             ::  fn_val

    REAL(DP), PARAMETER  ::  r = 3.442620_DP
    REAL(DP)             ::  x, y

    IF( .NOT. initialized ) CALL zigset( jsr )
    hz = shr3( )
    iz = IAND( hz, 127 )
    IF( ABS( hz ) < kn(iz) ) THEN
        fn_val = hz * wn(iz)
    ELSE
        DO
            IF( iz == 0 ) THEN
                DO
                    x = -0.2904764_DP* LOG( uni( ) )
                    y = -LOG( uni( ) )
                    IF( y+y >= x*x ) EXIT
                END DO
                fn_val = r+x
                IF( hz <= 0 ) fn_val = -fn_val
                RETURN
            END IF
            x = hz * wn(iz)
            IF( fn(iz) + uni( )*(fn(iz-1)-fn(iz)) < EXP(-half*x*x) ) THEN
                fn_val = x
                RETURN
            END IF
            hz = shr3( )
            iz = IAND( hz, 127 )
            IF( ABS( hz ) < kn(iz) ) THEN
                fn_val = hz * wn(iz)
                RETURN
            END IF
        END DO
    END IF
    RETURN
    END FUNCTION rnor



    !  Generate random exponentials
    FUNCTION rexp( ) RESULT( fn_val )
    REAL(DP)  ::  fn_val

    REAL(DP)  ::  x

    IF( .NOT. initialized ) CALL Zigset( jsr )
    jz = shr3( )
    iz = IAND( jz, 255 )
    IF( ABS( jz ) < ke(iz) ) THEN
        fn_val = ABS(jz) * we(iz)
        RETURN
    END IF
    DO
        IF( iz == 0 ) THEN
            fn_val = 7.69711 - LOG( uni( ) )
            RETURN
        END IF
        x = ABS( jz ) * we(iz)
        IF( fe(iz) + uni( )*(fe(iz-1) - fe(iz)) < EXP( -x ) ) THEN
            fn_val = x
            RETURN
        END IF
        jz = shr3( )
        iz = IAND( jz, 255 )
        IF( ABS( jz ) < ke(iz) ) THEN
            fn_val = ABS( jz ) * we(iz)
            RETURN
        END IF
    END DO
    RETURN
    END FUNCTION rexp

    END MODULE ziggurat
#endif 



    module Random
    integer :: rand_inst = 0
    logical, parameter :: use_ziggurat = .false.
    integer, parameter :: krand = KIND(1.d0)
    !Ziggurat is significantly (3-4x) faster, see Wikipedia for details
    !Have seem some suspicious things, though couldn't replicate; may be OK..

    INTERFACE RandRotation
    MODULE PROCEDURE RandRotationS, RandRotationD
    END INTERFACE

    contains

    subroutine initRandom(i, i2)
    use AMLUtils
#ifdef ZIGGURAT
    use Ziggurat
#endif
    implicit none
    integer, optional, intent(IN) :: i
    integer, optional, intent(IN) :: i2
    integer seed_in,kl,ij
    character(len=10) :: fred
    real(krand) :: klr

    if (present(i)) then
        seed_in = i
    else
        seed_in = -1
    end if
    if (seed_in /=-1) then
        if (present(i2)) then
            kl=i2
            if (i2 > 30081) call MpiStop('initRandom:second seed too large')
        else
            kl = 9373
        end if
        ij = i
    else
        call system_clock(count=ij)
        ij = mod(ij + rand_inst*100, 31328)
        call date_and_time(time=fred)
        read (fred,'(e10.3)') klr
        kl = mod(int(klr*1000), 30081)
    end if

    if (Feedback > 0 ) write(*,'(" Random seeds:",1I6,",",1I6," rand_inst:",1I4)') ij,kl,rand_inst
    call rmarin(ij,kl)
#ifdef ZIGGURAT
    if (use_ziggurat) call zigset(ij)
#endif
    end subroutine initRandom

    subroutine RandIndices(indices, nmax, n)
    use AMLUtils
    integer, intent(in) :: nmax, n
    integer indices(n),i, ix
    integer tmp(nmax)

    if (n> nmax) call MpiStop('Error in RandIndices, n > nmax')
    do i=1, nmax
        tmp(i)=i
    end do
    do i=1, n
        ix = int(ranmar()*(nmax +1 -i)) + 1
        indices(i) = tmp(ix)
        tmp(ix) = tmp(nmax+1-i)
    end do

    end subroutine RandIndices

    subroutine RandRotationS(R, N)
    !this is most certainly not the world's most efficient or robust random rotation generator
    integer, intent(in) :: N
    real R(N,N), vec(N), norm
    integer i,j

    do j = 1, N
        do
            do i = 1, N
                vec(i) = Gaussian1()
            end do
            do i = 1, j-1
                vec = vec - sum(vec*R(i,:))*R(i,:)
            end do
            norm = sum(vec**2)
            if (norm > 1e-3) exit
        end do
        R(j,:) = vec / sqrt(norm)
    end do

    end subroutine RandRotationS


    subroutine RandRotationD(R, N)
    !this is most certainly not the world's most efficient or robust random rotation generator
    integer, intent(in) :: N
    double precision R(N,N), vec(N), norm
    integer i,j

    do j = 1, N
        do
            do i = 1, N
                vec(i) = Gaussian1()
            end do
            do i = 1, j-1
                vec = vec - sum(vec*R(i,:))*R(i,:)
            end do
            norm = sum(vec**2)
            if (norm > 1e-3) exit
        end do
        R(j,:) = vec / sqrt(norm)
    end do

    end subroutine RandRotationD


    double precision function GAUSSIAN1()
#ifdef ZIGGURAT
    use Ziggurat
#endif
    implicit none
    double precision R, V1, V2, FAC
    integer, save :: iset = 0
    double precision, save :: gset

    if (use_ziggurat) then
#ifdef ZIGGURAT
        Gaussian1 = rnor( )
#endif
    else
        !Box muller
        if (ISET==0) then
            R=2
            do while (R >= 1.d0)
                V1=2.d0*ranmar()-1.d0
                V2=2.d0*ranmar()-1.d0
                R=V1**2+V2**2
            end do
            FAC=sqrt(-2.d0*log(R)/R)
            GSET=V1*FAC
            GAUSSIAN1=V2*FAC
            ISET=1
        else
            GAUSSIAN1=GSET
            ISET=0
        endif
    end if
    end function GAUSSIAN1


    double precision function CAUCHY1()
    implicit none

    Cauchy1 = Gaussian1()/max(1d-15,abs(Gaussian1()))

    end function CAUCHY1


    real FUNCTION RANDEXP1()
    !
    !     Random-number generator for the exponential distribution
    !     Algorithm EA from J. H. Ahrens and U. Dieter,
    !     Communications of the ACM, 31 (1988) 1330--1337.
    !     Coded by K. G. Hamilton, December 1996, with corrections.
    !
    real u, up, g, y

    real, parameter ::   alog2= 0.6931471805599453
    real, parameter ::      a = 5.7133631526454228
    real, parameter ::      b = 3.4142135623730950
    real, parameter ::     c = -1.6734053240284925
    real, parameter ::      p = 0.9802581434685472
    real, parameter ::     aa = 5.6005707569738080
    real, parameter ::     bb = 3.3468106480569850
    real, parameter ::     hh = 0.0026106723602095
    real, parameter ::     dd = 0.0857864376269050

    u = ranmar()
    do while (u.le.0)                 ! Comment out this block
        u = ranmar()                    ! if your RNG can never
    enddo                             ! return exact zero
    g = c
    u = u+u
    do while (u.lt.1.0)
        g = g + alog2
        u = u+u
    enddo
    u = u-1.0
    if (u.le.p) then
        randexp1 = g + aa/(bb-u)
        return
    endif
    do
        u = ranmar()
        y = a/(b-u)
        up = ranmar()
        if ((up*hh+dd)*(b-u)**2 .le. exp(-(y+c))) then
            randexp1 = g+y
            return
        endif
    enddo

    end function randexp1


    ! This random number generator originally appeared in ''Toward a Universal
    ! Random Number Generator'' by George Marsaglia and Arif Zaman.
    ! Florida State University Report: FSU-SCRI-87-50 (1987)
    !
    ! It was later modified by F. James and published in ''A Review of Pseudo-
    ! random Number Generators''
    !
    ! THIS IS THE BEST KNOWN RANDOM NUMBER GENERATOR AVAILABLE.
    !    (However, a newly discovered technique can yield
    !        a period of 10^600. But that is still in the development stage.)
    !
    ! It passes ALL of the tests for random number generators and has a period
    !   of 2^144, is completely portable (gives bit identical results on all
    !   machines with at least 24-bit mantissas in the floating point
    !   representation).
    !
    ! The algorithm is a combination of a Fibonacci sequence (with lags of 97
    !   and 33, and operation "subtraction plus one, modulo one") and an
    !   "arithmetic sequence" (using subtraction).
    !
    ! On a Vax 11/780, this random number generator can produce a number in
    !    13 microseconds.
    !========================================================================
    !
    !      PROGRAM TstRAN
    !     INTEGER IJ, KL, I
    ! Thee are the seeds needed to produce the test case results
    !      IJ = 1802
    !      KL = 9373
    !
    !
    ! Do the initialization
    !      call rmarin(ij,kl)
    !
    ! Generate 20000 random numbers
    !      do 10 I = 1, 20000
    !         x = RANMAR()
    !10    continue
    !
    ! If the random number generator is working properly, the next six random
    !    numbers should be:
    !          6533892.0  14220222.0  7275067.0
    !    6172232.0  8354498.0   10633180.0
    !
    !
    !
    !      write(6,20) (4096.0*4096.0*RANMAR(), I=1,6)
    !20    format (3f12.1)
    !      end
    !
    subroutine RMARIN(IJ,KL)
    ! This is the initialization routine for the random number generator RANMAR()
    ! NOTE: The seed variables can have values between:    0 <= IJ <= 31328
    !                                                      0 <= KL <= 30081
    !The random number sequences created by these two seeds are of sufficient
    ! length to complete an entire calculation with. For example, if sveral
    ! different groups are working on different parts of the same calculation,
    ! each group could be assigned its own IJ seed. This would leave each group
    ! with 30000 choices for the second seed. That is to say, this random
    ! number generator can create 900 million different subsequences -- with
    ! each subsequence having a length of approximately 10^30.
    !
    ! Use IJ = 1802 & KL = 9373 to test the random number generator. The
    ! subroutine RANMAR should be used to generate 20000 random numbers.
    ! Then display the next six random numbers generated multiplied by 4096*4096
    ! If the random number generator is working properly, the random numbers
    !    should be:
    !           6533892.0  14220222.0  7275067.0
    !           6172232.0  8354498.0   10633180.0
    double precision U(97), C, CD, CM, S, T
    integer I97, J97,i,j,k,l,m
    integer ij,kl
    integer ii,jj


    !      INTEGER IRM(103)

    common /RASET1/ U, C, CD, CM, I97, J97
    if( IJ .lt. 0  .or.  IJ .gt. 31328  .or. &
        KL .lt. 0  .or.  KL .gt. 30081 ) then
    print '(A)', ' The first random number seed must have a value  between 0 and 31328'
    print '(A)',' The second seed must have a value between 0 and   30081'
    stop
    endif
    I = mod(IJ/177, 177) + 2
    J = mod(IJ    , 177) + 2
    K = mod(KL/169, 178) + 1
    L = mod(KL,     169)
    do 2 II = 1, 97
        S = 0.0
        T = 0.5
        do 3 JJ = 1, 24
            M = mod(mod(I*J, 179)*K, 179)
            I = J
            J = K
            K = M
            L = mod(53*L+1, 169)
            if (mod(L*M, 64) .ge. 32) then
                S = S + T
            endif
            T = 0.5 * T
3       continue
        U(II) = S
2   continue
    C = 362436.0 / 16777216.0
    CD = 7654321.0 / 16777216.0
    CM = 16777213.0 /16777216.0
    I97 = 97
    J97 = 33

    end subroutine RMARIN

    double precision function RANMAR()
    ! This is the random number generator proposed by George Marsaglia in
    ! Florida State University Report: FSU-SCRI-87-50
    ! It was slightly modified by F. James to produce an array of pseudorandom
    ! numbers.
    double precision U(97), C, CD, CM
    integer I97, J97
    double precision uni

    common /RASET1/ U, C, CD, CM, I97, J97
    !      INTEGER IVEC
    UNI = U(I97) - U(J97)
    if( UNI .lt. 0.0 ) UNI = UNI + 1.0
    U(I97) = UNI
    I97 = I97 - 1
    if(I97 .eq. 0) I97 = 97
    J97 = J97 - 1
    if(J97 .eq. 0) J97 = 97
    C = C - CD
    if( C .lt. 0.d0 ) C = C + CM
    UNI = UNI - C
    if( UNI .lt. 0.d0 ) UNI = UNI + 1.0 ! bug?
    RANMAR = UNI

    end function RANMAR


    end module Random


