How to make a random function in fortran to generate the same random distribution into array?

I think my code below it's not exactly give me the same random distribution.

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

According to my code, I try to random value into l1(1,:,:,:), l1(2,:,:,:), l2(1,:,:,:) and l2(2,:,:,:) where l(t, x, y, z) is (3+1)-dimension array

Why i use trigonometry function for my random function? because i want a circular randomization. If i plot distribution of l1(1,:,:,:) vs l2(1,:,:,:) or l1(2,:,:,:) vs l2(2,:,:,:) you will see that a shape of plotting is circle.

So, let's get into my problem. Why i tell you that this's not exactly give me a same distribution? because i was tried to measure a variance of them and i got

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

where (variance_l1_t2 - variance_l1_t1) = (variance_l2_t1 - variance_l1_t2) = (variance_l2_t2 - variance_l2_t1) = 0.00166

That's quite weird. In actaully i should get almost the same variance value of l1(1,:,:,:), l1(2,:,:,:), l2(1,:,:,:) and l2(2,:,:,:)

How to solve this problem?