% Longitude limits of plot
% Can set lower bound to negative value for a plot that
% spans the prime meridian
lons = [120, 180];

% Latitude limits of plot
lats = [20, 50];

% Density of wind vectors
% must be a positive integer
% smaller = more dense
skip = 1;

% Scale of wind vectors
% larger = longer
scale = 1.5;

% Load near-surface currents
lat = double(ncread('world_oscar_vel_clim_1993_2012.nc', 'latitude'));
lon = double(ncread('world_oscar_vel_clim_1993_2012.nc', 'longitude'));
u = double(ncread('world_oscar_vel_clim_1993_2012.nc', 'u'));
v = double(ncread('world_oscar_vel_clim_1993_2012.nc', 'v'));

% Sparsify data based on skip parameter
u = u(1:skip:end,1:skip:end);
v = v(1:skip:end, 1:skip:end);
lat = lat(1:skip:end);
lon = lon(1:skip:end);

% Load sea surface height data
% (already time-averaged)
toplat = double(ncread( ...
    'hawaii_soest_a695_f334_f99c_5420_a7bb_a872.nc', 'latitude'));
toplon = double(ncread( ...
    'hawaii_soest_a695_f334_f99c_5420_a7bb_a872.nc', 'longitude'));
mdot = double(ncread( ...
    'hawaii_soest_a695_f334_f99c_5420_a7bb_a872.nc', 'mdot'));

% Set up map
figure('Position', [0, 0, 800, 400])
h = worldmap(lats, lons);
axis off;
setm(h, 'MapProjection','pcarree', 'MLabelParallel', 'south', ...
    'MLineLocation', 2, 'PLineLocation', 2, 'GLineStyle', '-');

% Plot sea surface height
contourfm(toplat, toplon, mdot', -100:20:100);
colormap(redblue);
caxis([-100, 100])
h = colorbar('northoutside');
ylabel(h, 'Sea surface height anomaly (cm)')

% Plot wind vectors
[LON, LAT] = meshgrid(lon, lat);
quiverm(LAT, LON, v', u', 'k', scale);

% Plot coastlines
load coastlines;
plotm(coastlat, coastlon, 'k');
tightmap;

% Save figure
exportgraphics(gca,'surface_currents.png')

% Generate a red-blue colormap
function c = redblue(m)
if nargin < 1, m = size(get(gcf,'colormap'),1); end
if (mod(m,2) == 0)
    m1 = m*0.5;
    r = (0:m1-1)'/max(m1-1,1);
    g = r;
    r = [r; ones(m1,1)];
    g = [g; flipud(g)];
    b = flipud(r);
else
    m1 = floor(m*0.5);
    r = (0:m1-1)'/max(m1,1);
    g = r;
    r = [r; ones(m1+1,1)];
    g = [g; 1; flipud(g)];
    b = flipud(r);
end
c = [r g b]; 
end