Try approximately 250 lines of code:
subroutine rand_flat_01S(rand_val)
!------------------------------------------------------------------
! RAND_FLAT_01S
!
! Generate a double-precision uniformly distributed random variate
! in the interval (0, 1).
!------------------------------------------------------------------
use, intrinsic :: iso_fortran_env, only : real32, real64
real(real32) :: rand_val
real(real64) :: temp
call park_miller_uniformD(temp)
rand_val = realDtoS(temp)
end subroutine rand_flat_01S
subroutine rand_flat_01D(rand_val)
!------------------------------------------------------------------
! RAND_FLAT_01D
!
! Generate a double-precision uniformly distributed random variate
! in the interval (0, 1).
!------------------------------------------------------------------
use, intrinsic :: iso_fortran_env, only : real64
real(real64) :: rand_val
call park_miller_uniformD(rand_val)
end subroutine rand_flat_01D
subroutine rand_flat_abS(A, B, randAB)
!------------------------------------------------------------------
! RAND_FLAT_ABS
!
! Generate a single-precision uniformly distributed random variate
! in the interval [A, B].
!------------------------------------------------------------------
use, intrinsic :: iso_fortran_env, only : real32
implicit none
real(real32), intent(in) :: A, B
real(real32), intent(out) :: randAB
real(real32) :: rand01
call rand_flat_01S(rand01)
randAB = A + (B - A) * rand01
end subroutine rand_flat_abS
subroutine rand_flat_abD(A, B, randAB)
!------------------------------------------------------------------
! RAND_FLAT_ABD
!
! Generate a double-precision uniformly distributed random variate
! in the interval [A, B].
!------------------------------------------------------------------
use, intrinsic :: iso_fortran_env, only : real64
implicit none
real(real64), intent(in) :: A, B
real(real64), intent(out) :: randAB
real(real64) :: rand01
call rand_flat_01D(rand01)
randAB = A + (B - A) * rand01
end subroutine rand_flat_abD
subroutine rand_gaussS(randval)
!------------------------------------------------------------------
! RAND_GAUSSS
!
! Generate a single-precision standard Gaussian (mean 0, std 1)
! random variate using the Box-Muller transform.
!------------------------------------------------------------------
use, intrinsic :: iso_fortran_env, only : real32, real64
real(real32) :: randval
real(real64) :: temp
call rand_gaussD(temp)
randval = realDtoS(temp)
end subroutine rand_gaussS
subroutine rand_gaussD(randval)
!------------------------------------------------------------------
! RAND_GAUSSD
!
! Generate a double-precision standard Gaussian (mean 0, std 1)
! random variate using the Box-Muller transform.
!------------------------------------------------------------------
use, intrinsic :: iso_fortran_env, only : real64
real(real64) :: randval
real(real64) :: u1, u2, eps, twopi
twopi = 2.0_real64*4.0_real64*atan(1.0_real64)
eps = epsilon(1.0_real64)
do
call park_miller_uniformD(u1)
if (u1 .lt. eps ) cycle
call park_miller_uniformD(u2)
randval = sqrt(-2.0_real64 * log(u1)) * cos(twopi * u2)
exit
end do
end subroutine rand_gaussD
subroutine rand_gauss_mu_sigmaS(mu, sigma, randX)
!------------------------------------------------------------------
! RAND_GAUSS_MU_SIGMAS
!
! Generate a single-precision Gaussian (normal) random variate with
! mean MU and standard deviation SIGMA.
!------------------------------------------------------------------
use, intrinsic :: iso_fortran_env, only : real32
implicit none
real(real32), intent(in) :: mu
real(real32), intent(in) :: sigma
real(real32), intent(out) :: randX
real(real32) :: z
call rand_gaussS(z)
randX = mu + sigma * z
end subroutine rand_gauss_mu_sigmaS
subroutine rand_gauss_mu_sigmaD(mu, sigma, randX)
!------------------------------------------------------------------
! RAND_GAUSS_MU_SIGMAD
!
! Generate a double-precision Gaussian (normal) random variate with
! mean MU and standard deviation SIGMA.
!------------------------------------------------------------------
use, intrinsic :: iso_fortran_env, only : real64
implicit none
real(real64), intent(in) :: mu
real(real64), intent(in) :: sigma
real(real64), intent(out) :: randX
real(real64) :: z
call rand_gauss(z)
randX = mu + sigma * z
end subroutine rand_gauss_mu_sigmaD
subroutine rand_gauss_mu_sigma_trunc_abS(mu, sigma, a, b, randX)
!------------------------------------------------------------------
! RAND_GAUSS_MU_SIGMA_TRUNC_ABDS
!
! Generate a single-precision Gaussian random variate with mean
! MU and standard deviation SIGMA, truncated to the interval [A,B].
!
! Rejection sampling is used; if no accepted sample is found after
! MAX_ITER attempts, a uniform fallback in [A,B] is returned.
!------------------------------------------------------------------
use, intrinsic :: iso_fortran_env, only : real32
implicit none
real(real32), intent(in) :: mu, sigma, a, b
real(real32), intent(out) :: randX
real(real32) :: x, lo, hi
integer :: iter
integer, parameter :: max_iter = 100000
lo = min(a, b)
hi = max(a, b)
do iter = 1, max_iter
call rand_gauss_mu_sigma(mu, sigma, x)
if (x >= lo .and. x <= hi) then
randX = x
return
end if
end do
! Fallback: uniform in [a,b]
call rand_flat_abS(lo, hi, randX)
end subroutine rand_gauss_mu_sigma_trunc_abS
subroutine rand_gauss_mu_sigma_trunc_abD(mu, sigma, a, b, randX)
!------------------------------------------------------------------
! RAND_GAUSS_MU_SIGMA_TRUNC_ABD
!
! Generate a double-precision Gaussian random variate with mean
! MU and standard deviation SIGMA, truncated to the interval [A,B].
!
! Rejection sampling is used; if no accepted sample is found after
! MAX_ITER attempts, a uniform fallback in [A,B] is returned.
!------------------------------------------------------------------
use, intrinsic :: iso_fortran_env, only : real64
implicit none
real(real64), intent(in) :: mu, sigma, a, b
real(real64), intent(out) :: randX
real(real64) :: x, lo, hi
integer :: iter
integer, parameter :: max_iter = 100000
lo = min(a, b)
hi = max(a, b)
do iter = 1, max_iter
call rand_gauss_mu_sigma(mu, sigma, x)
if (x >= lo .and. x <= hi) then
randX = x
return
end if
end do
print*, iter
! Fallback: uniform in [a,b]
call rand_flat_abD(lo, hi, randX)
end subroutine rand_gauss_mu_sigma_trunc_abD
subroutine park_miller_uniformD(r)
!-----------------------------------------------------------------------
! Function : park_miller_uniformD
! Returns : A reproducible pseudo-random number in the interval [0, 1)
!
! Description:
! Stand-alone version of the Park-Miller "Minimal Standard" RNG
! with Bays-Durham shuffle. This version returns REAL(REAL64)
! values, with fixed internal seed and no arguments.
!
! Properties:
! - Returns REAL(REAL64) in [0, 1)
! - Deterministic and platform-independent
! - Self-contained and reproducible
! - Overall cycle length before repeating is after 2.1 billion calls
!
! KSS
! 05/07/2025
!-----------------------------------------------------------------------
use, intrinsic :: iso_fortran_env, only: real64, int32
implicit none
real(real64), intent(out) :: r
! integer constants (32-bit arithmetic)
integer(int32), parameter :: ia = 16807 ! ia - 7**5 (primitive root of M_31)
integer(int32), parameter :: im = 2147483647 ! im = 2**31 - 1 (Mersenne prime, M_31)
integer(int32), parameter :: iq = 127773
integer(int32), parameter :: ir = 2836
integer(int32), parameter :: ntab = 32
integer(int32), parameter :: ndiv = 1 + (im - 1) / ntab
! real constants in real64 precision
real(real64), parameter :: am = 1.0_real64 / real(im, real64)
real(real64), parameter :: rnmx = 1.0_real64 - epsilon(1.0_real64)
! rng state
integer(int32), save :: seed = 1
integer(int32), save :: iv(ntab) = 0
integer(int32), save :: iy = 0
integer(int32) :: j, k
! initialization
if (iy .eq. 0) then
seed = max(seed, 1)
do j = ntab+8, 1, -1
k = seed / iq
seed = ia * (seed - k * iq) - ir * k
if (seed .lt. 0) seed = seed + im
if (j .le. ntab) iv(j) = seed
end do
iy = iv(1)
end if
! main generation logic
k = seed / iq
seed = ia * (seed - k * iq) - ir * k
if (seed .lt. 0) seed = seed + im
j = 1 + iy / ndiv
if (j .gt. ntab) j = ntab
iy = iv(j)
iv(j) = seed
r = am * real(iy, real64)
if (r .ge. 1.0_real64) r = rnmx
end subroutine park_miller_uniformD