function flow
clear 
clear path 
whos
global Nx Ny Nz
global lat long dz dm mdep
global delt_su su_its t_su delt
global descriptor this_path
global f deltaf Q beta r_expt r_heat H
global time rots it
global g Cp rho_bar alpha
global u v t w 
global iterations

figure;k=1;
        uxy=u(:,:,k);
        ubx=(uxy(1:(Nx-1),:)+uxy(2:Nx,:))*0.5;
        ubx(Nx,:)=(uxy(Nx,:)+uxy(1,:))*0.5;
        uxy(find(uxy==0))=uxy(find(uxy==0))*NaN;
        ubx(find(ubx==0))=ubx(find(ubx==0))*NaN;
        % V
        vxy=v(:,:,k);
        vby=(vxy(:,1:(Ny-1))+vxy(:,2:Ny))*0.5;
        vby(:,Ny)=(vxy(:,1)+vxy(:,Ny))*0.5;
        vxy(find(vxy==0))=vxy(find(vxy==0))*NaN;
        vby(find(vby==0))=vby(find(vby==0))*NaN;

        % Scale arrow
        ubx(6,6) = 0.01;
        vby(6,6) = 0.0;
%        [x,y]=meshgrid(long(2:4:Ny-2),lat(2:4:Nx-2));
        [x,y]=meshgrid(long,lat);
%        quiver(x,y,ubx((2:4:Nx-2),(2:4:Ny-2))',vby((2:4:Nx-2),(2:4:Ny-2))',2); axis square
        quiver(x,y,ubx,vby,2); axis square
        title(['speed timestep ' num2str(eval(it)) ' level ' int2str(k) ]);
        text(6,8,'1 cm/s')
        set(gca,'Xlim',[0,max(long)]);
        set(gca,'Ylim',[0,max(lat)]);
        drawnow
return