%% usage
load gfs_20161003_1200.mat 
p_vort=p(6:end);

lev=find(p==850);
pcolor(lon,lat,squeeze(phi(lev,:,:)));shading('flat')
axis([-72-6/cos(15*pi/180),-72,12,18])
%axis([-85,-65,10,30])
hold on
plot(lonm,latm,'w')
quiver(lon,lat,squeeze(u(lev,:,:)),squeeze(v(lev,:,:)),3,'w')
hold off
colorbar
return

filename="http://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/Global_0p25deg/GFS_Global_0p25deg_20161003_1200.grib2"
ncf=netcdf(filename,'r')
%ncvar(ncf)
lat=ncf{'lat'}(:);
lon=ncf{'lon'}(:);
p=ncf{'isobaric'}(:)/100;
p4=ncf{'isobaric4'}(:)/100;
latr=241:361;
lonr=1081:1201;
lat=lat(latr);
lon-lon(lonr)-360;

t=1
zeta_a=squeeze(ncf{"Absolute_vorticity_isobaric"}(t,:,latr,lonr));
phi=squeeze(ncf{"Geopotential_height_isobaric"}(t,:,latr,lonr));
slp=squeeze(ncf{"Pressure_reduced_to_MSL_msl"}(t,latr,lonr));
T0=squeeze(ncf{"Temperature_surface"}(t,latr,lonr));
T=squeeze(ncf{"Temperature_isobaric"}(t,:,latr,lonr));
omega=squeeze(ncf{"Vertical_velocity_pressure_isobaric"}(t,:,latr,lonr));
u=squeeze(ncf{"u-component_of_wind_isobaric"}(t,:,latr,lonr));
v=squeeze(ncf{"v-component_of_wind_isobaric"}(t,:,latr,lonr));

return

