LAPACK: ZHEEV routine failing for large matrices

I have a Fortran 90 code which is finding the eigenvalues of a spin-chain Hamiltonian after block-diagonalizing. I diagonalize each block as they are generated, and this seems to work just fine until the size of the block gets too large, at which point I get the error.

Exception thrown at 0x05E8A247 (mkl_avx2.dll) in spinchain.exe: 0xC0000005: Access violation writing location 0x054E8000.

where I am using the Intel Fortran compiler with MKL libraries in Visual Studio 2017. From what I can tell it appears that when the order of the block matrix is larger than around 60 I get this error. Running in Release mode (as opposed to Debug) the code will run all the way through without any complaints, but I assume the output is incorrect.

I wrote a short code to create an N x N hermitian matrix and run it through the zheev routine. For simplicity I chose the matrix of all ones, which has N-1 zero eigenvalues and a single eigenvalue N. I simply loop through values of N, generating the matrix and diagonalizing it with zheev(), printing N, the max eigenvalue, and the sum of all eigenvalues. I find that everything is fine up to N=51, at which point the maximum eigenvalue returns 52 and then I get an exception when I try to deallocate the arrays with the error Critical error detected c0000374.

Hitting "continue" the exception turns into

Unhandled exception at 0x77A5A879 (ntdll.dll) in Eigtest.exe: 0xC0000374: A heap has been corrupted (parameters: 0x77A95910).

Code below:

program Eigenvalues
    implicit none
    integer, parameter :: dp = kind(0.d0)
    integer :: N, lwork, info, i
    real(dp), allocatable :: lambda(:), work(:), rwork(:)
    complex(dp), allocatable :: H(:,:)

    do N = 1,70

        lwork = 3*N
        allocate(H(N,N))
        allocate(lambda(N))
        allocate(work(lwork))
        allocate(rwork(lwork))

        H = dcmplx(1.0_dp,0)

        call zheev('N','U',N,H,N,lambda,work,lwork,rwork,info)

        if (info == 0) then
                write(*,'(I5, F16.12, F16.12)'), N, lambda(N), sum(lambda)
        else
            print*, "diagonalization failed: info = ", info
            read*, i
            stop
        end if

        deallocate(H,lambda,work,rwork)

    end do

    read*, i

end program

I also wrote another version of the code where I fix N as a parameter at the beginning and define the arrays without allocatable, which will run without any exceptions but gives the wrong eigenvalues for N=51 and most values above N=51.

Forgive me if my code is bad, I'm a physicist, not a programmer. Comments or suggestions appreciated.

According to the man page of zheev, the array work is expected to be complex*16, so we probably need something like

real(dp), allocatable :: lambda(:), rwork(:)
complex(dp), allocatable :: H(:,:), work(:)

(this seems to be working with gfortran on my computer).