w0 = 2; zR = Pi*w0^2; rmax = 30; zmax = 100;
plot = Plot3D[
   Exp[-2 r^2/(w0^2 (1 + (z/zR)^2))]/(1 + (z/zR)^2),
   {r, -rmax, rmax}, {z, -zmax, zmax}, PlotRange -> {0, 1},
   PlotPoints -> 101, MaxRecursion -> 4, 
   Method -> {Refinement -> {ControlValue -> 0.01} },
   MeshFunctions -> {#1 &, #2 &, #3 &, Exp[-2 #1^2/(w0^2 (1 + (#2/zR)^2))] &}, 
   MeshStyle -> {Black, Black, Gray, White},
   Mesh -> {Range[-Floor[rmax, 10], Floor[rmax, 10], 10],
     Range[-Floor[zmax, 10], Floor[zmax, 10], 10], 
     Range[0.05, 0.95, 0.05], {1/2, 1/E^2} },
   ColorFunction -> Function[{r, z, i}, Lighter[ColorData["Rainbow"][Sqrt[i] ], 0.3] ],
   BoxRatios -> {rmax, zmax, 30},
   ViewPoint -> {1.2, -3, 0.7},
   AxesEdge -> { {-1, -1}, {1, -1}, {-1, -1} }, 
   AxesLabel -> {"r/\[Lambda]", "z/\[Lambda]", "I/\!\(\*SubscriptBox[\(I\), \(0\)]\)"},
   BaseStyle -> 13, 
   Epilog -> 
    Text["\!\(\*SubscriptBox[\(w\), \(0\)]\) = 2\[Lambda]", {-15, -3}, BaseStyle -> Gray] ]
Export["Gaussian-beam_intensity_surfaceplot.png", 
  ImageResize[Image[plot, ImageResolution -> 400], 2000, Resampling -> "Linear"] ];