Fast way to calculate source response using Green's function
Performance wise

In general:

You use numpy, but you write it almost like in fortran. Python with numpy is good for scientific programing and computing as long as you don't do many loops


read this first A beginners guide to using Python for performance computing
if you really need to do tight loops than use
cython.
But most of the time you can avoid that by expressing the loop as some functional operation on the whole array 


most numpy functions can operator on whole array in one pass, use it.


always prefer y[:]=np.sqrt(x[:]) and avoid




          `for i in xrange(N): 
                y[i]=np.sqrt(x[i])`



operation like dot,convolve,einsum,ufunc.reduce are very potent to express most sciencetific algorithms in more abstract way
use advanced array slicing


you can express most inner if by boolen mask arrays or integer index array


e.g. f[ f&gt;15.0 ] = 15.0 clamp all elements in array f which are &gt;15.0 to 15.0; you can also store the boolean mask mask = (f&gt;15.0) and than use it like f[mask]+=g[mask]. This way you can express branching as fully functional programs/expressions with arrays.
Do not construct new np.array (np.zeros, np.ones etc.) too often (i.e. inside tight loops). Optimal performance is obtained if you prepare all arrays you need at the beggining. Avoid frequent conversion between list and array




To address your code example in particular:


Green's functions are basically convolutions. I'm pretty sure you can express it using e.g. scipy.ndimage.filters.convolve
if your convolution kernel is large (i.e. pixels interact more than with few neighbors) than it is often much faster to do it in Fourier space (convolution transforms as multiplication) and using np.fftn with O(nlog(N)) cost.
cannot def gw(xx, yy, tt, mm, nn, ss): operate on the whole array,
rather than individual numbers?
if denom &lt; 0: can be expressed using boolean mask array like mask = denom &gt;0; denom[mask] = np.sqrt(denom[maks])
Never do this if you can:


  for k in range (Nt):
     for j in range (Ny):
        for i in range (Nx):


the operations seems to be possible to rewrite that operate on the whole x,y space at once   


This is how you do convolution-like operations in numpy (from here)


           # The actual iteration
           u[1:-1, 1:-1] = ((u[0:-2, 1:-1] + u[2:, 1:-1])*dy2 +
                            (u[1:-1,0:-2] + u[1:-1, 2:])*dx2)*dnr_inv



this is so horrible:


green = np.array([gw(x[i],y[j],t[k],i_grid,j_grid,k_grid) for k_grid in t for j_grid in y for i_grid in x])



list comprehesion is relatively fast, but still much slower than numpy array operations (which are implemented in C).
do not create temporary list and convert it to temporary array, you loose lot time doing that.
why you cannot just prealocate Green[Nx,Ny,:] and  Gauss[Nx,Ny,:] and use it as a whole?