How to calculate the gradient in cylindrical coordinates in more efficient way?

It is not clear if $k_z$ depends on z or not. In both cases, you can try writing in the following way:

f[r_, [Theta]_, z_] := 
 BesselJ[0, r Sqrt[[Omega]^2/c^2 - k[z]^2]] S[k[z]] Exp[I k[z] z];
Grad[f[r, [Theta], z], {r, [Theta], z}, "Cylindrical"]

with $k[z]$ taking into account for the dependance.

The output is the gradient in cylindrical coordinates

{-E^(I k z) BesselJ[1, r Sqrt[[Omega]^2/9 - k[z]^2]] Sqrt[[Omega]^2/9 - k[z]^2]S[k[z]],
 0, 
 (E^(I z k[z])
   r BesselJ[1, r Sqrt[[Omega]^2/9 - k[z]^2]] k[z] S[k[z]]
 Derivative[1][k][z])/Sqrt[[Omega]^2/9 - k[z]^2] + BesselJ[0, r Sqrt[[Omega]^2/9-k[z]^2]]
 S[k[z]] (I k[z] + I z Derivative[1][k][z]) + 
 E^(I z k[z]) BesselJ[0, r Sqrt[[Omega]^2/9 - k[z]^2]] Derivative[1][k][z] Derivative[1][S][k[z]]}

and it is quite complicated because k depends on z. If k does not depend on z, the output is:

{-E^(I k z) BesselJ[1, r Sqrt[[Omega]^2/9 - k^2]] Sqrt[[Omega]^2/9 - k^2] S[k],
 0, 
 I E^(I k z) k BesselJ[0, r Sqrt[[Omega]^2/9 - k^2]] S[k]}