Azimuthal average in Fortran? Find indexes in Fortran?

I am working on an eigenvalue problem in fortran. I have used Lapack to solve the problem and get the eigenvalues and eigenvectors. This is done for $201\times101$ wavenumbers, only half the wavespace due to symmetry, and for each grid-point in a large domain(in the ocean). I am searching for the maximum eigenvalue for each gridpoint, and I would like to not just pick the absolute maximum in the $201\times101$ matrix of eigenvalues, but perform an azimuthal average in wavespace - and then pick the maximum average. I am struggeling with seeing how to do this.

At first I coded it like this:

! Wavenumber domain

dx=4000.
pi = 4.*atan(1.)

DO m=1,ktot
kx(m) = -(2.*pi)/(dx) + ((m-1)*2.*pi)/(100.*dx)
END DO

DO l=1,ltot
ly(l)= ((l-1)*2.*pi)/(100.*dx)
END DO


DO m=1,ktot
DO l=1,ltot
END DO
END DO


! Azimuthal average (omegai are the eigenvalues I have found, a ktot*ltot large matrix)

DO i=1,ltot-1
length=count(ind)
WHERE (ind) average_omegai = sum(omegai)/length
END DO


But it seems I am summing a horizontal chunk for all wavenumbers in k-direction, between ly(i) and ly(i+1). I rather need to make a semicircle in wavespace to sum all omegai-values in between it. Can anyone help with this?

(BTW: I also tried

 DO m=1,ktot-1
DO l=1,ltot-1
length=count(ind)
WHERE (ind) average_omegai = sum(omegai)/length
END DO
END DO


but I right now I am blind to what I am doing wrong! My problem stated differently would be: is there another way in Fortran to find indexes in a matrix which fulfill some criteria? The find-function in Matlab would do the work, but I have not found an equivalent intrinsic function in Fortran? )