如何在fortran中创建随机函数以将相同的随机分布生成数组？

我认为下面的代码并不能完全给我相同的随机分布。

subroutine trig_random_value()
    implicit none
    integer :: t, z, y, x
    real(real64) :: theta, r
    real(real64), parameter :: PI=4.D0*DATAN(1.D0)

    integer, dimension(12) :: date_time
    integer, dimension(12) :: seed

    call date_and_time(values=date_time)
    call random_seed
    seed = date_time(6) * date_time(7) + date_time(8)
    call random_seed(put = seed)

    do z = 1, z_size
        do y = 1, y_size
            do x = 1, x_size

                theta = rand()*2*PI
                r = 0.1*rand()
                l1(1, z, y, x) = r*cos(theta)
                l2(1, z, y, x) = r*sin(theta)

                theta = rand()*2*PI
                r = 0.1*rand()
                l1(2, z, y, x) = r*cos(theta)
                l2(2, z, y, x) = r*sin(theta)

            end do
        end do
    end do

    return
end subroutine trig_random_value

根据我的代码，我尝试随机值l1(1,:,:,:)，l1(2,:,:,:)，l2(1,:,:,:)和l2(2,:,:,:)其中l（T，X，Y，Z）为（3 + 1） -尺寸数组

为什么我对随机函数使用三角函数？因为我要循环随机化。如果我绘制l1（1，：，：，）vs l2（1，：，：，）或l1（2，：，：，）vs l2（2，：，：:)的分布图，我会获得半径为0.1的圆形分布

那么，为什么我要告诉你，这并不能完全给我相同的分配？因为我试图测量它们的方差，我得到了

 variance_l1_t1 =   1.6670507752921395E-003
 variance_l1_t2 =   3.3313151655785292E-003
 variance_l2_t1 =   4.9965623815717321E-003
 variance_l2_t2 =   6.6641054728288360E-003

请注意，（variance_l1_t2-方差_l1_t1）=（variance_l2_t1-方差_l1_t2）=（variance_l2_t2-方差_l2_t1）= 0.00166

那是一个很奇怪的结果。在实际上我应该得到几乎相同的方差值l1(1,:,:,:)，l1(2,:,:,:)，l2(1,:,:,:)和l2(2,:,:,:)如果这个功能，如果良好的随机函数。可能是我做错了。

如何解决这个问题呢？

根据要求提供的其他信息：

real(real64) function find_variance(l)
    implicit none
    real(real64), dimension(z_size, y_size, x_size), intent(in) :: l
    integer :: z, y, x
    real(real64) :: l_avg = 0
    real(real64) :: sum_val = 0

    do z = 1, z_size
        do y = 1, y_size
            do x = 1, x_size

                l_avg = l_avg + l(z, y, x)

            end do
        end do
    end do

    l_avg = l_avg/(z_size*y_size*x_size)

    do z = 1, z_size
        do y = 1, y_size
            do x = 1, x_size

                sum_val = sum_val + (l(z , y, x) - l_avg)**2

            end do
        end do
    end do

    find_variance = sum_val/(z_size*y_size*x_size)

    return
end function find_variance