It is a tricky question. `CountourPlot[F[x,y]==0,...]` finds points where `F[x,y]>0` and `F[x,y]<0`. Then by dichotomy it finds points where `F[x,y]` is approximately zero. In your case `Abs[f[x,y]]-0.001` is almost always positive, so the algorithm fails. 

It is nontrivial (at least at the first sight) to find points when both the real and the imaginary part of the complex function is zero. Especially if the function has a quickly oscillating phase.

For this moment I found the following workaround:

It finds points where derivatives of `Abs[f[x,y]]` have a singularities or are zero.

    f[a_, a0_, k_, K0_] := 
      a^2 Sech[(a a0)/2]^2 (-2 I (1 + E^(2 I a0 k)) k + 
          a (-1 + E^(2 I a0 k)) Tanh[(a a0)/2]) + 
       2 k (I E^(I a0 k) ((a - k) (a + k) Cos[a0 k] + (a^2 + k^2) Cos[
               a0 K0]) + a (-1 + E^(2 I a0 k)) k Tanh[(a a0)/2]);
    With[{a0 = 10., a = 1.4}, 
     ContourPlot[{Abs[f[a, a0, y, x - 0.01]] == Abs[f[a, a0, y, x + 0.01]], 
       Abs[f[a, a0, y - 0.01, x]] == Abs[f[a, a0, y + 0.01, x]]}, 
        {x, -\[Pi]/a0, \[Pi]/a0}, {y, 0, 2}]]

![enter image description here][1]

In plot only those points are meaningful where both colors match.