A part from efficiency, I noticed that with the definition below, lmax >= 64 works:
Clear[field]; field[θ_, ϕ_] := Chop@ Total[ Table[ alms[l, m] SphericalHarmonicY[l, m, θ, ϕ], {l, 0, lmax}, {m, -l, l} ] , 2 ]; nn = 4.; dat = ParallelTable[ field[θ, ϕ], {θ, 0, Pi, Pi/nn}, {ϕ, 0., 2 Pi, 2 Pi/nn/2} ]; Round[Re[dat], 0.01] { {0.51, 0.51, 0.51, 0.51, 0.51, 0.51, 0.51, 0.51, 0.51},
{0.02, -1.54, -1.48, -0.47, -2.29, 0.61, 2.05, 1.82, 0.02},
{0.24, -1.61, 0.44, 0.32, -0.55, 0.65, -1.12, -0.08, 0.24},
{-0.34, 0.63, -0.54, 2.08, -0.72, 1.09, -2.09, -1.3, -0.34},
{0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3} }
which looks better and produces a reasonable plot.
So, I tested also:
Clear[field]; field[θ_, ϕ_] := Sum[ alms[l, m] SphericalHarmonicY[l, m, θ, ϕ], {l, 0, lmax}, {m, -l, l} ]; and this too seems to me to work.
Consequently, it appears to me that Compile over Sum is doing some inappropriate manipulation.
