matrix - Parallelization Fortran 90 -

i started learn how parallelize programs in fortran , trying include in program.

the problem using lapack (to calculate eigenvalues , eigenvectors) , apparently creating conflict parallelization routine, because depending on matrix size calculates eigenvalues (wrong values small sizes, , doesn't compute them @ matrices equal or bigger 10x10).

when run program 1 thread works correctly.

program nuevo  use omp_lib      implicit none   integer ::i,j,k,l,kk,jj,red,ok,info,lwork,stat1    integer :: m,n,q,rea,i0,j0,chunk,chunksize   !!!! rea= number of realizations >1   parameter(m=10,n=10,rea=12,lwork=2*m*n,ok=m*n,chunksize=rea/2)     complex*8 t,aa,tt,ww !! caution dimension of variables   integer, dimension(1)::  seed, old   integer, parameter ::  ikind=selected_real_kind(p=15)   real (kind=ikind) ::  pi   real ( kind = 8 ) wtime   real*8 d, dd(m,n),zz, x(2*m*n*rea), xx(m*n,2*rea), ss, hh, harvest(2*n*m*rea), w   complex*16 a(m*n,m*n), work(3*n-2,3*n-2)   real*8 rwork(2*m*n), b(m*n), f   !!!!!!!!!!!!random numbers!!!!!!!!!!!!!!!!!!!      seed(1) = 67937     call random_seed     call random_seed(size=k)     call random_seed(get=old(1:k))     call random_number(harvest)         x=harvest   !!!!!!!!!!!random matrix!!!!!!!!!!!!!!!!!!!!!!        xx=reshape( x, (/ m*n, 2*rea /) ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! wtime = omp_get_wtime ( ) !call omp_set_num_threads(1)    t=1.0d0                  !!  hoppings    chunk=chunksize  !$omp parallel !$omp schedule(static,chunk) ordered     jj=0,2*(rea-1),2  !!for realizations   w=0.1d0              !! desorder ratio   ww=n   q=m*n   aa=(0.0d0,1.0d0) !! complex number    a=0d0    b=0d0    f=2/ww     j=1,(m*n)     if (xx(j,jj+1)<0.5)        ss=1.0d0      else        ss=-1.0d0      end if     a(j,j)=ss*0.5d0*w*xx(j,jj+2)        !!  diagonal terms slighty random modifications      end      !!!!!!!!!!!!!!!!!!!!!!!!!!!     !!!!!!  matrix    !!!!!!!!!     !!!!!!!!!!!!!!!!!!!!!!!!!!!    !! in x direction  (i odd)          i=1,m,2             j=1,n-2,2                a((j*m)+i,((j+1)*m)+i)=t             end          end   !!  in x direction  (i even)          i=2,m,2             j=1,n-1,2                a(((j-1)*m)+i,(j*m)+i)=t             end          end    !!  in y direction             j=1,n             i=1,m-1             a(((j-1)*m)+i,((j-1)*m)+i+1)=t*exp(pi*j*f*aa) !!  phase factor due magnetic field              end          end !! periodic boundaries in x direction            i=1,m,2            a(i,(i+m*n-m))=t          end  !!    periodic boundaries in y direction              j=1,n                a(((j-1)*m)+1,j*m)=t*(exp(-pi*j*f*aa))              end  !!=========conjugate===============           i0=1,m*n             j0=1,m*n                a(j0,i0)=dconjg(a(i0,j0))             end          end  call zheev('v', 'u', ok, a, ok, b, work, lwork, rwork, info) !!lapack routine     i=1,m*n         write(*,*) b(i)         write(*,*) ' '     end  end     !! realizations   !$omp end !$omp end parallel  end program nuevo