As @mecej4 observed, the code can be simplified if one is only generating uniform variates. I have implemented pure subroutines (not functions) to generate single values or sequences of random 32-bit integers or uniform variates, with an `integer :: intent(in out)`

argument representing the state:

```
module ziggurat_pure_mod
! Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating random variables', J. Statist. Software, v5(8).
implicit none
private
public :: dp, random_int_32, random_uni
interface random_int_32
module procedure random_int_32_scalar, random_int_32_vec
end interface random_int_32
interface random_uni
module procedure random_uni_scalar,random_uni_vec
end interface random_uni
integer, parameter :: dp = selected_real_kind(12, 60)
contains
!
pure elemental subroutine random_int_32_scalar(jsr,iran)
! generate random 32-bit integers
integer, intent(in out) :: jsr ! state of RNG
integer, intent(out) :: iran ! random integer
integer :: jz
jz = jsr
jsr = ieor(jsr, ishft(jsr, 13))
jsr = ieor(jsr, ishft(jsr, -17))
jsr = ieor(jsr, ishft(jsr, 5))
iran = jz + jsr
end subroutine random_int_32_scalar
!
pure subroutine random_int_32_vec(jsr,iran)
! generate random 32-bit integers
integer, intent(in out) :: jsr ! state of RNG
integer, intent(out) :: iran(:) ! random integers
integer :: i
do i=1,size(iran)
call random_int_32_scalar(jsr,iran(i))
end do
end subroutine random_int_32_vec
!
pure elemental subroutine random_uni_scalar(jsr,xran)
integer , intent(in out) :: jsr ! state of RNG
real(kind=dp), intent(out) :: xran ! random uniform variate
integer :: iran
call random_int_32(jsr,iran)
xran = 0.2328306e-9_dp*iran + 0.5_dp
end subroutine random_uni_scalar
!
pure subroutine random_uni_vec(jsr,xran)
integer , intent(in out) :: jsr ! state of RNG
real(kind=dp), intent(out) :: xran(:) ! random uniform variates
integer :: i
do i=1,size(xran)
call random_uni_scalar(jsr,xran(i))
end do
end subroutine random_uni_vec
end module ziggurat_pure_mod
```

```
program xzig_pure
! 10/27/2021 12:54 AM driver for random_int_32 and random_uni
use ziggurat_pure_mod, only: dp, random_int_32, random_uni
implicit none
integer, parameter :: seed = 123, nran = 10000000
integer :: iran(nran),state
real(kind=dp) :: xran(nran)
state = seed
call random_int_32(state,iran)
call random_uni(state,xran)
write (*,"(2a12)") "min","max"
write (*,"(2i12)") minval(iran),maxval(iran)
write (*,"(/,4a10)") "min","max","mean","mean_sq"
write (*,"(4f10.6)") minval(xran),maxval(xran),sum(xran)/nran,sum(xran**2)/nran
end program xzig_pure
```

Compiling with `gfortran ziggurat_pure.f90 xzig_pure.f90`

or the analogous, gfortran, Intel Fortran, flang, and g95 all give output

```
min max
-2147483306 2147482694
min max mean mean_sq
0.000000 1.000000 0.500141 0.333496
```