w: vorticity field
function vm_ic(nx,ny,x,y,w)
xc1 = pi-pi/4.0 # horizontal location of left vortex
yc1 = pi # vertical location of left vortex
xc2 = pi+pi/4.0 # horizontal location of right vortex
yc2 = pi # vertical location of right~vortex