ghost points
w[1,:] = w[nx+1,:]
w[:,1] = w[:,ny+1]
for i = 1:nx for j = 1:ny
sf[i,j] = wf[i,j]/k2[i,j] # streamfunction in Fourier space
end~end
s[2:nx+1, 2:ny+1] = real(ifft(sf)) # convert streamfunction to physical space