function [UP, VP, psi1]=sGstep(u0, v0, psi0, h)

% Input u0, v0, psi0, h
% Output UP, VP, psi1
UR = u0([ 2:end 1 ]);
UL = u0([ end 1:end-1 ]);
VR = v0([ 2:end 1 ]);
VL = v0([ end 1:end-1 ]);
pR = psi0([ 2:end 1 ]);
pL = psi0([ end 1:end-1 ]);
fR = -sin(pL);
fL = -sin(pR);

UP = (UR+UL)/2 + (VR-VL)/2 + (fR-fL)*h/2;
VP = (UR-UL)/2 + (VR+VL)/2 + (fR+fL)*h/2;

p1 = pL + (UL + UP + VL + VP)*h/2;
p2 = pR + (-UR - UP + VR + VP)*h/2;
psi1 = (p1+p2)/2;
