generate_bernoulli_mask Subroutine

  subroutine generate_bernoulli_mask(this)
    !! Generate Bernoulli mask
    implicit none

    ! Arguments
    class(dropblock3d_layer_type), intent(inout) :: this
    !! Instance of the 3D dropblock layer

    ! Local variables
    real(real32), allocatable, dimension(:,:,:) :: mask_real
    !! Real mask
    integer :: i, j, k
    !! Loop indices
    integer, dimension(2) :: ilim, jlim, klim
    !! Limits for mask


    ! Generate Bernoulli mask
    !---------------------------------------------------------------------------
    ! assume random number already seeded and don't need to again
    allocate(mask_real(size(this%mask,1), size(this%mask,2), size(this%mask,3)))
    call random_number(mask_real)  ! Generate random values in [0..1]

    this%mask = .true. ! 1 = keep

    !! Apply threshold to create binary mask
    !---------------------------------------------------------------------------
    do k = 1 + this%half, size(this%mask, dim=3) - this%half
       do j = 1 + this%half, size(this%mask, dim=2) - this%half
          do i = 1 + this%half, size(this%mask, dim=1) - this%half
             if(mask_real(i, j, k).lt.this%gamma)then
                ilim(:) = [ &
                     max(i - this%half, lbound(this%mask,1)), &
                     min(i + this%half, ubound(this%mask,1)) ]
                jlim(:) = [ &
                     max(j - this%half, lbound(this%mask,2)), &
                     min(j + this%half, ubound(this%mask,2)) ]
                klim(:) = [ &
                     max(k - this%half, lbound(this%mask,3)), &
                     min(k + this%half, ubound(this%mask,3)) ]
                this%mask( &
                     ilim(1):ilim(2), &
                     jlim(1):jlim(2), &
                     klim(1):klim(2)) = .false. ! 0 = drop
             end if
          end do
       end do
    end do

  end subroutine generate_bernoulli_mask