For those that may have an interest in random number generation in SF, I've posted a basic program in SF that uses random numbers.
Even though it doesn't address the original post here, it might be instructive for those interested. Although the code used has been modified by me, it comes from various miscellaneous sources and is not original. I've named it WH_RAND after the  WICHMAN-HILL method.

!***********************************************************************
!*  WH_RAND: UNIFORM RANDOM NUMBER GENERATOR VERSION 3.0               *
!*                                                                     *
!*  RANF IS A UNIFORM RANDOM NUMBER ROUTINE, ORIGIN IS UNKOWN.         *
!*  RDUNWH IS USED IN PLACE OF RANF AS A UNIFORM RANDOM NUMBER ROUTINE.*
!*                                                                     *
!*  MODIFIED WITH RANDOM SEEDS RANSEED() TO INITIALIZE THE RANDOM SEED *
!*  USING THE SYSTEM CLOCK TO ALLOW FOR NEWLY GENERATED NUMBERS AFTER  *
!*  EACH EXECUTION OF THE ROUTINE.                                     *
!*                                                                     *
!*  REVISED AND MODIFIED BY F.W. PERRELLA, PH.D.  FEBRUARY 29, 2024    *
!***********************************************************************
!
!-----------------------------------------------------------------------
!---------- MODULES ----------------------------------------------------
!-----------------------------------------------------------------------
MODULE NRTYPES
!--- KIND TYPES FOR 64- AND 32-BIT SIGNED INTEGERS
    INTEGER, PARAMETER :: I32 = SELECTED_INT_KIND(9)
    !INTEGER, PARAMETER :: I64 = SELECTED_INT_KIND(18)

    INTEGER, PARAMETER :: SP = SELECTED_REAL_KIND(6)
    !INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(15)
END MODULE NRTYPES

MODULE IRANDOM
    USE NRTYPES
    INTEGER(I32) :: K2
    INTEGER(I32) :: IRAND(250)
    INTEGER(I32) :: IX,IY,IZ
END MODULE IRANDOM

MODULE WH_RAND
    USE NRTYPES
    IMPLICIT NONE

    PUBLIC :: RANSEED, RANSET, RANF, RDUNWH

    INTERFACE RANSEED
        MODULE PROCEDURE RANSEED
    END INTERFACE RANSEED

    INTERFACE RANSET
        MODULE PROCEDURE RANSET
    END INTERFACE RANSET

    INTERFACE RANF
        MODULE PROCEDURE RANF
    END INTERFACE RANF

    INTERFACE RDUNWH
        MODULE PROCEDURE RDUNWH
    END INTERFACE RDUNWH
    PRIVATE

    !-------------------------------------------------------------------
    CONTAINS
    !-------------------------------------------------------------------

    SUBROUTINE RANSEED(ISEED)
        IMPLICIT NONE
        INTEGER(I32), INTENT(OUT) :: ISEED
        INTEGER(I32)              :: I
        REAL(SP)                  :: A, AC, AM, R1
        PARAMETER(A = 4.5875, AC = 45381.6693, AM = 214748.3648)

        CALL INIT_RANDOM_SEED(ISEED)

        I = 0
        DO WHILE ((ISEED == 0).AND.(I < 10))
            I = I + 1
            CALL INIT_RANDOM_SEED(ISEED)
        END DO

        IF(ISEED .EQ. 0) ISEED = 10000
        IF(ISEED .LT. 0) ISEED = -ISEED

        R1 = 1.0 + MODULO(A*REAL(ISEED) + AC, AM)
        ISEED = INT(R1)

        RETURN
    END SUBROUTINE RANSEED

    !--- INITIALIZE A RANDOM SEED FROM THE SYSTEM CLOCK AT EVERY RUN.
    SUBROUTINE INIT_RANDOM_SEED(ISEED)
        IMPLICIT NONE
        INTEGER(I32), INTENT(OUT) :: ISEED
        INTEGER(I32)              :: I, N, ICLOCK
        INTEGER(I32), DIMENSION(:), ALLOCATABLE :: RSEED

        CALL RANDOM_SEED(SIZE = N)
        ALLOCATE(RSEED(N))

        CALL SYSTEM_CLOCK(COUNT = ICLOCK)

        IF (ICLOCK < 0) ICLOCK = 30000
        RSEED = ICLOCK + 37 * (/ (I - 1, I = 1, N) /)

        CALL RANDOM_SEED(PUT = RSEED)

        ISEED = INT( RSEED(N) )

        DEALLOCATE(RSEED)

        RETURN
    END SUBROUTINE INIT_RANDOM_SEED

    SUBROUTINE RANSET(ISD)
        USE IRANDOM

        IMPLICIT NONE
        INTEGER(I32), INTENT(IN OUT) :: ISD
        INTEGER(I32)                 :: I, J, K, JRAND, J1

        ISD = MODULO(ISD, 65536)

        DO I = 1, 250
            DO J = 1, 5
                ISD = MODULO(ISD*8973 + 1, 65536) - 32768
            END DO

            JRAND = ISD*65536
            J1 = MODULO(I, 2)
            DO K = 1, 5
                ISD = MODULO(ISD*8973 + J1, 65536)
            END DO

            IRAND(I) = JRAND + ISD
        END DO

        K2 = 1

        RETURN
    END SUBROUTINE RANSET

    FUNCTION RANF()
        USE IRANDOM

        IMPLICIT NONE
        REAL(SP)     :: RANF
        INTEGER(I32) :: L1, IRANF

        L1 = MODULO(K2 + 146, 250) + 1
        IRANF = IEOR(IRAND(K2), IRAND(L1))

        RANF = REAL(IRANF)*2.328306E-10 + 0.5
        IRAND(K2) = IRANF
        K2 = MODULO(K2, 250) + 1

        RETURN
    END FUNCTION RANF

    FUNCTION RDUNWH (IX,IY,IZ)
        !--- GENERATES UNIFORM PSEUDO-RANDOM DEVIATES IN THE
        !--- INTERVAL [0,1) USING THE WICHMAN-HILL METHOD.
        !--- AN EXACT VALUE OF ZERO IS POSSIBLE BUT NOT LIKELY.
        !--- THE CYCLE LENGTH EXCEEDS 6.95E12.
        !---
        !--- CURRENT VERSION COMPLETED MAY 15, 1987
        !---
        !--- REFERENCE:
        !--- GRIFFITHS, P. AND HILL, I.D. (EDS.),
        !--- 'APPLIED STATISTICS ALGORITHMS',
        !--- THE ROYAL STATISTICAL SOCIETY/ELLIS HARWOOD LIMITED,
        !--- ALGORITHM AS 183, PP. 238-242 (1985).
        !---
        !--- THE INTEGERS IX, IY, AND IZ REQUIRE INTEGER
        !--- ARITHMETIC UP TO 5,212,632
        !--- BEFORE THE FIRST CALL TO THIS ROUTINE THE INTEGERS
        !--- IX, IY, AND IZ SHOULD BE SET TO VALUES WITHIN THE
        !--- RANGE [1,30000]. AFTER THAT THEIR VALUES SHOULD NOT
        !--- BE CHANGED EXCEPT BY THIS ROUTINE.
        IMPLICIT NONE
        REAL(SP)                     :: RDUNWH, SUMX
        INTEGER(I32), INTENT(IN OUT) :: IX, IY, IZ

        IX = MOD(171*IX,30269)
        IY = MOD(172*IY,30307)
        IZ = MOD(170*IZ,30323)
        SUMX = REAL(IX)/30269.0+REAL(IY)/30307.0+REAL(IZ)/30323.0

        RDUNWH = AMOD(SUMX,1.0)

        RETURN
    END FUNCTION RDUNWH

END MODULE WH_RAND
!-----------------------------------------------------------------------
!---------- END MODULES ------------------------------------------------
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!---------- TEST PROGRAM -----------------------------------------------
!-----------------------------------------------------------------------
PROGRAM RANDOM_TEST
    USE NRTYPES
    USE WH_RAND

    IMPLICIT NONE
    CHARACTER(1) :: CR
    INTEGER(I32) :: I, J, K, CNT, N, ISEED, IX,IY,IZ, IXYZ
    INTEGER(I32) :: NBINS, BIN, YBIN, IFACTOR, ICYCLE, GCOUNT
    REAL(SP)     :: YMIN, YMAX, YMEAN, YSIGMA
    REAL(SP)     :: SUMN, SUMSQ, COVCNT, COVSCORE
    REAL(SP)     :: Y, RANGEX, W, GMEAN, GSIGMA, GSCORE
    PARAMETER(IFACTOR = 100, ICYCLE = 5)

    ALLOCATABLE :: Y(:), IXYZ(:), BIN(:), YBIN(:), RANGEX(:)
    ALLOCATE(Y(1), IXYZ(3))
    ALLOCATE(BIN(IFACTOR), YBIN(IFACTOR), RANGEX(IFACTOR))

    !--- INITIALIZE SEED VALUE
    CALL RANSEED(ISEED)

    DO I = 1, 3
        !--- RESET SEED VALUE
        CALL RANSET(ISEED)

        !--- CHECK THAT THE VALUES IN IN RANGE.
        IF (ISEED > 30000) THEN
            ISEED = 30000*INT(FLOAT(ISEED - 30000)/FLOAT(ISEED))
        END IF
        IF (ISEED == 0) ISEED = 13849 / I
        IXYZ(I) = ISEED
    END DO

    !--- SEEDS FOR RDUNWH() RANDOM NUMBER GENERATOR
    IX = IXYZ(1)
    IY = IXYZ(2)
    IZ = IXYZ(3)

    !--- GENERATE N*N = 10,000 UNIFORM RANDOM NUMBERS
    N = 100

1   CONTINUE
    GMEAN = 0.0
    GSIGMA = 0.0
    GSCORE = 0.0
    GCOUNT = 0.0

    IF (ALLOCATED(Y)) DEALLOCATE (Y)
    ALLOCATE(Y(N*N))

    WRITE(*,"(/1X,A)") 'WH_RAND UNIFORM RANDOM NUMBER GENERATOR'

    WRITE(*,"(/1X,A,3I6)") 'INITIAL RDUNWH SEEDS IX,IY,IZ =', IXYZ(1:3)
    WRITE(*,"(A)") REPEAT('~',40)

    WRITE(*,"(/A)") REPEAT('-',72)

    !--- CALCULATE ICYCLE MEAN VALUES
    DO K = 1, ICYCLE
        !--- RESET SEED VALUES
        J = RDUNWH(IX,IY,IZ)

        !--- STORE NEW SEEDS
        IXYZ(1) = IX
        IXYZ(2) = IY
        IXYZ(3) = IZ

        CNT = 0
        DO J = 1, N
            DO I = 1, N
                CNT = CNT + 1
                !--- GENERATE AND SAVE RANDOM DATA Y
                !--- Y(CNT) = RANF()
                Y(CNT) = RDUNWH(IX,IY,IZ)
            END DO
        END DO

        !--- CALCULATE THE MINUMUM AND MAXIMUM VALUES
        YMIN = MINVAL(Y)
        YMAX = MAXVAL(Y)

        !--- CALCULATE THE MEAN AND STANDARD DEVIATION
        SUMN = SUM( Y(:) )
        SUMSQ = SUM( Y(:)**2 )

        YMEAN = SUMN / FLOAT(CNT)

        !--- YSIGMA = SQRT( SUM( (Y(:) - YMEAN)**2 ) / CNT )
        YSIGMA = SQRT( SUMSQ/FLOAT(CNT) - YMEAN**2 )

        !--- The covariance score is computed over each pair of the
        !--- data set Y, so that Y1 is paired with Y2, Y3 with Y4, etc.
        COVCNT = 0.0
        DO I = 1, CNT-1
            COVCNT = COVCNT + (Y(I+1)-YMEAN)*(Y(I)-YMEAN)
        END DO
        COVSCORE = COVCNT / FLOAT(CNT-1)

        WRITE(*,"(1X,A,1X,I6,1X,A,F10.4,2X,A,F10.4,2X,A,F10.4)") &
            "COUNT", CNT, "MEAN=", YMEAN, " SIGMA=", YSIGMA, &
            " COV SCORE=", COVSCORE

        IF (ICYCLE > 1) THEN
            !--- GRAND MEAN, GRAND SIGMA, GRAND SCORE
            GMEAN = GMEAN + YMEAN
            GSIGMA = GSIGMA + YSIGMA
            GSCORE = GSCORE + COVSCORE
            GCOUNT = GCOUNT + CNT
        END IF
    END DO

    WRITE(*,"(A)") REPEAT('-',72)

    IF (ICYCLE > 1) THEN
        !--- GRAND AVERAGES
        GMEAN = GMEAN / FLOAT(ICYCLE)
        GSIGMA = GSIGMA / FLOAT(ICYCLE)
        GSCORE = GSCORE / FLOAT(ICYCLE)

        WRITE(*,"(1X,A,1X,I6,1X,A,F10.4,2X,A,F10.4,2X,A,F10.6)") &
            "COUNT", GCOUNT, "MEAN=", GMEAN, &
            " SIGMA=", GSIGMA, " COV SCORE=", GSCORE
    END IF

    !-------------------------------------------------------------------
    !--- THIS SECTION USES A SCORE ARRAY AND A RANGE ARRAY, COUNTS
    !--- THE NUMBER OF EACH SCORES IN EACH RANGE, AND DISPLAYS A
    !--- HORIZONTAL BAR HISTOGRAM.
    !-------------------------------------------------------------------
    !--- CALCULATE THE NUMBER OF BINS FOR Y VALUES
    W = 5.0 * (FLOAT(IFACTOR)*YSIGMA + 0.10) * (N*N)**(-1.0/3.0)

    !--- ADJUST VALUES BY IFACTOR FOR PLOTTING
    Y(:) = FLOAT(IFACTOR)*Y(:)
    YMIN = FLOAT(IFACTOR)*YMIN
    YMAX = FLOAT(IFACTOR)*YMAX
    NBINS = INT((YMAX - YMIN) / W)

    WRITE(*,"(/A,F8.3,2X,A,I3)") 'W=', W, '  NBINS=', NBINS

    !--- CALCULATE THE RANGE OF BINS STARTING AT Y MINIMUM
    RANGEX(:) = YMIN
    DO I = 1, NBINS
        IF ( I .EQ. 1 ) THEN
        RANGEX(I) = RANGEX(I) + W
        ELSE
            RANGEX(I) = RANGEX(I-1) + W
        END IF
    END DO

    WRITE(*,"(A)") REPEAT('~',40)

    !--- STORE THE Y VALUES IN THE APPROPRIATE BIN
    WRITE(*,"(/1X,A)") 'NUMERIC HISTOGRAM'
    DO I = 1, NBINS
        YBIN(I) = COUNT( Y(:) .LE. RANGEX(I) )
        IF( I .GT. 1 ) THEN
        BIN(I) = YBIN(I) - YBIN(I-1)
        ELSE
            BIN(I) = YBIN(I)
        ENDIF
        WRITE(*,'(F10.2,I10)') RANGEX(I), BIN(I)
    END DO

    WRITE(*,"(A)") REPEAT('~',40)

    !--- DISPLAY HORIZONTAL BAR HISTOGRAM
    WRITE(*,"(/1X,A)") 'HISTOGRAM PLOT'
    DO I = 1, NBINS
        J = INT( FLOAT(BIN(I))/ FLOAT(NBINS) / 1.0 )
        WRITE(*, "(F4.1, A, A)") &
        REAL(I), ": ", REPEAT("=", J)
    END DO

    WRITE(*,"(A)") REPEAT('~',40)

    !...END RANDOM NUMBER ROUTINE
    WRITE(*,'(/A$)') "PRESS '1' TO CONTINUE OR 'ENTER' TO QUIT:"
    READ(*,'(A)') CR
    IF (TRIM(CR) == "1") THEN
        GO TO 1
    ELSE
        STOP 'PROGRAM ENDED...'
    END IF

END PROGRAM RANDOM_TEST
!-----------------------------------------------------------------------
!---------- END TEST PROGRAM -------------------------------------------
!-----------------------------------------------------------------------

/*

   We use the xoshiro256** generator, a fast high-quality generator
   that:

   - passes TestU1 without any failures

   - provides a "jump" function making it easy to provide many
     independent parallel streams.

   - Long period of 2**256 - 1

   A description can be found at

   http://prng.di.unimi.it/

   or

   https://arxiv.org/abs/1805.01407

   The paper includes public domain source code which is the basis for
   the implementation below.

*/

Just to be clear, that is the algorithm the current runtime library is using.  This algorithm could always change if something better comes along.