How do people deal with resized grid steps while numerically integrating using discrete Fourier Transform?

I am trying to simulate light propagation on python using FFT following the Fresnel diffraction equation given on Wikipedia: The problem with this is that the output matrix from the DFT would be given in coordinates $$(p,q)$$, which would inversely scale with my variable propagation distance $$z$$. To use the new $$(x,y)$$ coordinates to account for the values outside of the FFT I would need to factor out at least the $$z$$ value from my coordinate grid to get a 1-to-1 match, but doing so would increase my grid steps with increasing $$z$$.

What would be the normal solution to something like this? My current idea would be to convert the $$h(x,y)$$ section to $$h(p,q)$$ and then apply the evolution again(and use $$z/2$$ for each evolution), which should take me back to the same grid size. However, that would require applying FFT twice, which I would prefer to prevent if possible.

Any suggestion would be appreciated.

EDIT: Having thought about it a bit more, the two step fft might not work (I need to think about this more). The reason being that $$p$$ and $$q$$ being dependent on $$z$$ seems problematic on integration.

EDIT2: Never mind, it would be fine in this case simply because each integration is fixed with fixed distance $$z$$. Therefore I can treat it like a constant.