% example code for plotting:
% (1) daily maps in the NetCDF file "RB_2025_maps.nc"
% and
% (2) daily time series in the NetCDF file "RB_2025_timeseries.nc"

clear; clc; close all;
% use "ncdisp" to look into files
fname1 = "RB2025_maps.nc";
fname2 = "RB2025_timeseries.nc";
ncdisp(fname1);
ncdisp(fname2);

%% ---------------------- (1) Plotting daily maps -------------------------
% load in time
t = ncread(fname1,'t');
% turn it into matlab datetime
dt = datetime(t,'ConvertFrom','datenum');

% load in sea surface chlorophyll (chl) and corresponding lat/lon
latchl = ncread(fname1,'latchl');
lonchl = ncread(fname1,'lonchl');
chl = ncread(fname1,'chl');
% load in curl-driven vertical velocity (wcurl) and corresponding lat/lon
latflux = ncread(fname1,'latflux');
lonflux = ncread(fname1,'lonflux');
wcurl = ncread(fname1,'wcurl');

%% -------- define a day to look for ---------
dtuse = datetime(2019,8,1); % August 1, 2019 for example
dtind = find(dt == dtuse);
% plot chl for that day
figure('position',[1 1 1000 800]);
pcolor(lonchl,latchl,chl(:,:,dtind)); 
shading flat; colormap('turbo'); clim([0 6]); 
text(-119.8,34.6,string(dtuse),'fontsize',30);
cbar = colorbar('northoutside'); cbar.Label.String = 'SSChl (mg/m^3)';
axis equal; xlim([-120.75 -119.25]); ylim([33.75 34.75]);
xlabel('longitude (°E)'); ylabel('latitude (°N)');
set(gca,'fontsize',25,'linew',2,'layer','top');
% plot wcurl as well
figure('position',[1 1 1000 800]);
pcolor(lonflux,latflux,wcurl(:,:,dtind).*86400); % *86400 to convert m/s to m/day
shading flat; colormap(cmocean('balance')); clim([-15 15]); 
text(-119.8,34.6,string(dtuse),'fontsize',30);
cbar = colorbar('northoutside'); cbar.Label.String = 'w_c_u_r_l (m/day)';
axis equal; xlim([-120.75 -119.25]); ylim([33.75 34.75]);
xlabel('longitude (°E)'); ylabel('latitude (°N)');
set(gca,'fontsize',25,'linew',2,'layer','top');

%% ------------------------------------------------------------------------
close all; clearvars -except fname1 fname2;

%% ------------------- (2) Plotting daily time series ---------------------
% load in time
t = ncread(fname2,'t');
% turn it into matlab datetime
dt = datetime(t,'ConvertFrom','datenum');

% load in spatial mean sea surface chlorophyll (chl) time series
chl = ncread(fname2,'chl');
% load in spatial mean curl-driven vertical velocity (Wcurl) times series
Wcurl = ncread(fname2,'Wcurl');
% apply a 27-day moving mean to smooth
chl = movmean(chl,27,'omitnan');
Wcurl = movmean(Wcurl,27,'omitnan');

% plot the time series
figure('position',[1 1 1200 500]);
% chl
subplot(2,1,1);
plot(dt,chl,'r-','linew',2);
hold on; box on; grid on; ylabel('\langleSSChl\rangle (mg/m^3)');
set(gca,'fontsize',20,'linew',2); ylim([0 10]);
% Wcurl
subplot(2,1,2);
plot(dt,Wcurl*86400,'k-','linew',2); % *86400 to convert m/s to m/day
hold on; box on; grid on; ylabel('\langlew\rangle_c_u_r_l (m/day)');
set(gca,'fontsize',20,'linew',2); ylim([0 3.5]);






