the general idea is that solve for u v velocity from d(PV)/dt=0 for the nodes on each PV contour. PV (potential vorticity = relative vorticity + beta term) follows the material conservation law.

is Log(r) still the Green function of:

laplace(stream function) = Pv - beta*y  ??

if so,  the stream function solution at point (x,y) would be the integration for (Pv - beta*y)*log(r) dx'dy' over all the domain of (x' y') coordinate system. Here r2 = (x-x')2+(y-y')2.

the stream function derivatives, with respect to (x,y), from the first term Pv *log(r) dx'dy' could be simplified by contour dynamics, since d(Pv)/dx or d(Pv)/dy exist only on the contour of PV(Zabusky 1978). However now the question is how to get d(integrate beta*y*log(r) dx'dy' over all the domain)/dx as a function of (x,y) ??