How to make a 3D plot of two integrals?

One option is to use SetDelay,

S1[t_?NumericQ, z_?NumericQ] := ((A*Subscript[τ, 0])/Sqrt[2*Pi])*
  NIntegrate[Exp[-(((c*ω + (-1 + c)*Subscript[ω, 0])*(2*I*z + 
            Subscript[τ, 0]^2*(c*ω + (-1 + c)*
                Subscript[ω, 0])))/(2*c^2))]*
         Exp[I*(P1 + P2*ω + P3*ω^2 + P4*ω^3)*z]*
    Exp[I*ω*t], {ω, -Infinity, Infinity}]

S2[t_?NumericQ, z_?NumericQ] := ((A*Subscript[τ, 0])/Sqrt[2*Pi])*
  NIntegrate[
   Exp[-(((c*ω + (-1 + c)*Subscript[ω, 0])*(2*I*z + 
            Subscript[τ, 
               0]^2*(c*ω + (-1 + c)*
                Subscript[ω, 0])))/(2*c^2))]*
         Exp[I*(G1 + G2*ω + G3*ω^2 + G4*ω^3)*z]*
    Exp[I*ω*t], {ω, -Infinity, Infinity}]

Plot3D[{Abs[S1[t, z]], Abs[S2[t, z]]}, {t, 0, 20}, {z, 0, 20}]

But the problem here is that the plot is not smooth (I used PlotPoints without any luck) and it is time consuming. I hope someone will come up with a better solution.

I'm on my Inactivate and Activate kick. Also, if the relationships are reasonanbly smooth and you want pictures, often creating a Table and using ListPlot can give you what you need with plenty of control and less computational effort.

As @zhk does, convert to numerical integration. First I broke out the integrand and simplified it so I could see what was going on.

S1integrand = Exp[-(((c*ω + (-1 + c)*Subscript[ω, 0])*(2*I*z + 
       Subscript[τ, 
          0]^2*(c*ω + (-1 + c)*
           Subscript[ω, 0])))/(2*c^2))]* Exp[I*(P1 + P2*ω + P3*ω^2 
       + P4*ω^3)*z]* Exp[I*ω*t] // FullSimplify;

S2integrand =  Exp[-(((c*ω + (-1 + c)*Subscript[ω, 0])*(2*I*z + 
       Subscript[τ, 
          0]^2*(c*ω + (-1 + c)*
           Subscript[ω, 0])))/(2*c^2))]*  Exp[I*(G1 + G2*ω + G3*ω^2  
           + G4*ω^3)*z]*  Exp[I*ω*t] // FullSimplify

Then create the S1 and S2 terms using NIntegrate and a specified method (standard struggled for t close to 20), but inactivate them.

S1 = (((A*Subscript[τ, 0])/Sqrt[2*Pi])*
     NIntegrate[S1integrand, {ω, -Infinity, Infinity}, 
     Method -> "LocalAdaptive"]) // Inactivate

S2 = (((A*Subscript[τ, 0])/Sqrt[2*Pi])*
     NIntegrate[\!\(TraditionalForm\`S2integrand\), {ω, \
     -Infinity, Infinity}, Method -> "LocalAdaptive"]) // Inactivate

Make tables, activating the terms once t and z are set (vice creating a function that only takes in Numeric quantities).

S1tab = Flatten[Table[{t, z, Abs[S1]}, {t, 0, 20}, {z, 0, 20}] // Activate, 1];
S2tab = Flatten[Table[{t, z, Abs[S2]}, {t, 0, 20}, {z, 0, 20}] // Activate, 1];

It took approximately 20 seconds per table. Plot them.

ListPlot3D[{S1tab, S2tab}]