clear all
close all

addpath(genpath('/Users/Mau/Documents/MATLABhelp_new')); 
% print_setup.m for saving plots
% cbdate, maximize
% spectra.m

%% Purpose: Compare w_ctd from SBE41cp on float (2-m binnned vs 1Hz data
% from .eng files) and dp/dt from 64Hz pressure sensor on MicroRider

% Directory for CTD data:
path_ctd='/Users/Mau/Documents/MBARIprojects/CPF/CPF003/Data/matfiles/ParagonTests/MICRORIDER/2021-10-8T16-19-7/'; % profile and w_ctd data in .mat files
% Directory for MR data:
path_mr='/Users/Mau/Documents/MBARIprojects/MR/TurbulenceData/ParagonTests/2021-10-08/';

path_ctd='/Users/Mau/Documents/MBARIprojects/CPF/CPF003/Data/matfiles/ParagonTests/MICRORIDER/2021-10-13/';
path_mr='/Users/Mau/Documents/MBARIprojects/MR/TurbulenceData/ParagonTests/2021-10-13/'; % only 1 prof returned data in .msg files

% path_ctd='/Users/Mau/Documents/MBARIprojects/CPF/CPF003/Data/matfiles/ParagonTests/MICRORIDER/2021-10-15T16-10-6/';
% path_mr='/Users/Mau/Documents/MBARIprojects/MR/TurbulenceData/ParagonTests/2021-10-15/';

mission_date='2021-10-15T16-10-6';
% mission_date='2021-10-19T16-3-9';
path_ctd=['/Users/Mau/Documents/MBARIprojects/CPF/CPF003/Data/matfiles/ParagonTests/MICRORIDER/',mission_date,'/'];
path_mr=['/Users/Mau/Documents/MBARIprojects/MR/TurbulenceData/ParagonTests/',mission_date(1:10),'/'];

mission_date='2021-11-16'; % mission_date='2021-10-16T?'; % Vernight mission, 18 profiles
mission_date='2021-11-16T19-16-49'; % preserving folder name, % Overnight mission, 18 profiles
path_ctd=['/Users/Mau/Documents/MBARIprojects/CPF/CPF003/Data/matfiles_mr/ParagonTests/MICRORIDER/',mission_date,'/'];
path_mr=['/Users/Mau/Documents/MBARIprojects/MR/TurbulenceData/ParagonTests/',mission_date(1:10),'/matfiles/'];

% Define Directories to save merged data from CPF+MR
path_out=['/Users/Mau/Documents/MBARIprojects/CPF-MR/ParagonTests/',mission_date(1:10),'/matfiles/'];
% Define Directories to save plots (with MR data)
path_plots=[path_out(1:end-9),'plots','/'];

%% Create directory to save data & plots (if it doesn't exist already)
if isfolder(path_out)==0 & numel(dir(path_out))<=2 % will create directory only if it is empty!
    mkdir(path_out)
elseif isfolder(path_plots)==0 & numel(dir(path_plots))<=2 % will create directory only if it is empty!
    mkdir(path_plots)
end
%%
% NOTE: should consolidate folder naming conventions so data from the same
% mission matches. "done" for MR data.

%% Define colors for plots
cbar=cbrewer('seq','Blues',101);
% cbar=cbar(5:5:end,:);
blue=cbar(50,:);

Reds=cbrewer('seq','Reds',101);
red=Reds(50,:);

Greens=cbrewer('seq','Greens',101);
green1=Greens(40,:);
green2=Greens(60,:);

%% List files and loop through
files=dir([path_ctd,'*.mat']); % list matfiles with CPF data for a given mission
nfiles=length(files);
nfiles=nfiles-1; % Note, last CPF profile usually is garbage (if mission was recovered before writing data to file)

Lat=nan(nfiles,1);
Lon=nan(nfiles,1);
Time=nan(nfiles,1);

W_1hz=nan(800,nfiles);
P_1hz=nan(800,nfiles);

W_2mb=nan(200,nfiles);
D_2mb=nan(200,nfiles);
Nbin=nan(200,nfiles);

RHO=nan(200,nfiles);
CT=nan(200,nfiles);
SA=nan(200,nfiles);
P=nan(200,nfiles);

N2gsw=nan(200,nfiles);
Ngsw=nan(200,nfiles);
p_mid=nan(200,nfiles);
N2=nan(200,nfiles);
N=nan(200,nfiles);

mr_paxis=nan(300,nfiles);
W_mr=nan(300,nfiles);
P_mr=nan(300,nfiles);
T_mr=nan(300,nfiles);
epsilon=nan(300,nfiles,2);

% LR data
BETA_BACKSCATTERING700=nan(100,nfiles);
CHLA=nan(100,nfiles);
NITRATE=nan(100,nfiles);
PH=nan(100,nfiles);
DOXY=nan(100,nfiles);
PRES=nan(100,nfiles);
TEMP=nan(100,nfiles);
PSAL=nan(100,nfiles);

ppl=1; % do not plot profiles

for f=1:nfiles-1

    % NOte: should fix how file names when saved in .mat format to match
    % CPF matfile data ('.i.e. '0' extra)
    
    % Load data
    if f<10
        load([files(f).folder,'/','CPF003.00',num2str(f),'.mat']); % CPF data
        load([path_mr,'CPF003.00',num2str(f),'.MR.mat']); % MR data
    elseif (f>=10 && f<=99)
        load([files(f).folder,'/','CPF003.0',num2str(f),'.mat']); % CPF data
        load([path_mr,'CPF003.00',num2str(f),'.MR.mat']); % MR data
    end
    % Note: should fix how file names when saved in .mat format to match
    % CPF matfile data ('.i.e. '0' extra)

    maxp=100;
    
    % Compute W from SBE41cp 1Hz pressure signal (wihout smoothing, the RSI way)
    P1hz=ENG.P1hz;
    P1hz=P1hz(~isnan(P1hz));
    
    dP1Hz=diff(P1hz); dt=1;
    w1hz=-1*dP1Hz./dt;
    P1hz(w1hz==0)=[]; 
    w1hz(w1hz==0)=[]; %check w/Gene, is this because of the sensitivity of sensor or writing values twice in the .eng files?
    
    % NOTE: it makes sense to remove zeros (consecutive values with equal
    % pressure. Otherwise, a discontinuity will show up in the w profile at
    % ~80m)
    W_1hz(1:length(w1hz),f)=flipud(w1hz); %flip, such that surface is on top of matrix
    P_1hz(1:length(w1hz),f)=flipud(P1hz(1:end-1));
    
    % How does RSI calculate platform speed?
%     fc = 0.5;
%     fs_slow=63.99;
%     [b,a]=butter(4,fc/(fs_slow/2));
%     dP_dt = gradient(P_slow,1/fs_slow);
%     dP_dt_LP = filtfilt(b,a,dP_dt);
    
    % Save all profiles from agiven misison in a matrix
    W_2mb(1:length(HR.w_ctd),f)=HR.w_ctd;
    D_2mb(1:length(HR.w_ctd),f)=HR.d_axis;
    Nbin(1:length(HR.N),f)=HR.N;
    
    
    if ppl==1
    plot(w1hz,P1hz(1:end-1),'-o','color',green1); % computed doing the difference between samples at depth
    hold on
    plot(ENG.W,ENG.Plp,'-o','color',green2); % computed by Gene, SlopeLPFpressure
%     plot(ENG.Wlp,ENG.Plp,'-o','color','r'); % wrong units? these platVel
%     don't make much sense
%     plot(ENG.Wdiff,ENG.Plp,'-o','color','g');
    
    plot(MR.speed,MR.P,'-','color',red); %cbar(10,:)); % from MicroRider (LowPassed)
    plot(HR.w_ctd,HR.d_axis,'color',blue); %cbar(10,:)); %SBE41cp 2-m binned
    
    axis ij

    legend('SBE41cp 1Hz','SBE41cp 1Hz LP','MR P 64Hz LP0.5','SBE41cp 2m binned',...
        'location','W')

    xlim([range(HR.w_ctd(:))])
    ylim([0 maxp])
    xlabel('w_{dpdt} [m/s]');

    print_setup(20,40)
    print([path_plots,'w_2sensors_p_',num2str(f),'.jpg'],'-djpeg');

    close all
    
    end
    
    % Compute N2 (using gsw function)
    %     if (LAT{f}(p)<90 && LAT{f}(p)>-90) % quick fix to compute MLD where lat/lon is within range only
    g=9.81; % [m/s]
    [N2gsw(1:length(HR.P)-1,f), p_mid(1:length(HR.P)-1,f)] = gsw_Nsquared(HR.SA,HR.CT,HR.P,HR.LAT);
    
    nnan=find(~isnan(N2gsw(:,f)));
    Ngsw(nnan,f)=sqrt(N2gsw(nnan,f));
    
    % Compute N2 from density (i.e. to sort profile and avoid
    % negative N^2)
    drho=diff(sort(HR.RHO));
    dz=diff(HR.DEPTH);
    drhodz=drho./dz;
    rho_ref=nanmean(HR.RHO); % mean densitify for each profile
    N2(1:length(HR.P)-1,f)=drhodz.*[g./rho_ref];
    
    N(1:length(HR.P)-1,f)=sqrt(N2(1:length(HR.P)-1,f));
    % NOTE: already in HR! local version
    
    
    % Save epsilon and platform velocity from the MR
    ne=length(MR.e(2,:));
    epsilon(1:ne,f,1)=flipud(MR.e(1,:)');
    epsilon(1:ne,f,2)=flipud(MR.e(2,:)');
    P_mr(1:ne,f)=flipud(MR.P);
    T_mr(1:ne,f)=flipud(MR.T);
    W_mr(1:ne,f)=flipud(MR.speed);
    
    % Save Coordinates for plotting (from HR CPF structure)
    Lat(f)=HR.LAT;
    Lon(f)=HR.LON;
    % Time(f)=HR.TIME; % bug?; deployement 10/15 doesn't look good
    Time(f)=INFO.sdn;
    
    % SAVE other sensors variables:

    % HR CTD data
    RHO(1:length(HR.P),f)=HR.RHO;
    CT(1:length(HR.P),f)=HR.T;
    SA(1:length(HR.P),f)=HR.S;
    P(1:length(HR.P),f)=HR.P;
    
    % LR BGC data
    PRES(1:length(LR.PRES_QC),f)=flipud(LR.PRES_ADJUSTED(LR.PRES_ADJUSTED_QC<=1));
    TEMP(1:length(LR.PRES_QC),f)=flipud(LR.TEMP_ADJUSTED); % LR.TEMP_ADJUSTED_QC=NaN
    PSAL(1:length(LR.PRES_QC),f)=flipud(LR.PSAL_ADJUSTED); % LR.PSAL_ADJUSTED_QC

    NITRATE(1:length(LR.PRES_QC),f)=flipud(LR.NITRATE(LR.NITRATE_QC<=3));
    DOXY(1:length(LR.PRES_QC),f)=flipud(LR.DOXY(LR.DOXY_QC<=3));
    PH(1:length(LR.PRES_QC),f)=flipud(LR.PH_IN_SITU_TOTAL);
    CHLA(1:length(LR.PRES_QC),f)=flipud(LR.CHLA);
    BETA_BACKSCATTERING700(1:length(LR.PRES_QC),f)=flipud(LR.BETA_BACKSCATTERING700);
    
    % NOTE: Should apply QC flags here! and use ADJUSTED data, but many
    % variables were not populated. Need to revise and uncomment code in
    % calc_BGC.m

    save([path_out,files(f).name(1:end-3),'CPF-MR.mat'],'HR','LR','INFO','ENG','MR'); % save individual .mat files (for each profile)
    display(['Merging HR, LR, ENG and MR data... ',path_out])
end

%% SAVE all profiles in columns in one .mat file
% along with DERIVED VARIABLES
display(['Saving all profiles in one .mat file... ',path_out])
save([path_out,'CPF_MR_allprofiles.mat'],...
    'Lat','Lon','Time',...
    'W_1hz','W_2mb','W_mr','P_1hz','D_2mb','P_mr','P',...
    'epsilon','T_mr',...
    'RHO','CT','SA','N2','N','Nbin',...
    'NITRATE','DOXY','CHLA','BETA_BACKSCATTERING700','PH','PRES','TEMP','PSAL');

%% CODE for plotting data below, moved to plot_CPFMRdata_depthvstime.m


% %%  Load merged CPF-MR data
% load([path_out,'CPF_MR_allprofiles.mat']);
% %%
% % Calculate Local Time
% LocalTime=Time-8/24;
% 
% %%  Bin data in 2m bin before plotting
% [sz,np]=size(P);
% 
% % Define Variable to plot
% Z=NITRATE; Zname='Nitrate'; Zunits='[\mu mol/kg]';
% 
% % Define X axis (i.e. Time or Lon or Lat)
% [X,Y]=meshgrid(LocalTime,1:100);
% 
% %% Interpolate data to a common depth grid
% Zaxis=0:5:100';
% nz=length(Zaxis);
% 
% Pint=nan(nz,np);
% Nint=nan(nz,np);
% CHLAint=nan(nz,np);
% Bint=nan(nz,np);
% DOXYint=nan(nz,np);
% 
% Z2maxis=0:2:100';
% nz=length(Z2maxis);
% 
% CTint=nan(nz,np);
% SAint=nan(nz,np);
% RHOint=nan(nz,np);
% N2int=nan(nz,np);
% P2int=nan(nz,np);
% 
% for p=1:np
% 
%     % Interpolate LR data
%     indnani=find(~isnan(PRES(:,p)),1,'first');
%     indnanf=find(~isnan(PRES(:,p)),1,'last');
% 
%     Pint(:,p)=interp1(PRES(indnani:indnanf,p),PRES(indnani:indnanf,p),Zaxis,'nearest','extrap');
%     Nint(:,p)=interp1(PRES(indnani:indnanf,p),NITRATE(indnani:indnanf,p),Zaxis,'nearest','extrap');
%     CHLAint(:,p)=interp1(PRES(indnani:indnanf,p),CHLA(indnani:indnanf,p),Zaxis,'nearest','extrap');
%     Bint(:,p)=interp1(PRES(indnani:indnanf,p),BETA_BACKSCATTERING700(indnani:indnanf,p),Zaxis,'nearest','extrap');
%     DOXYint(:,p)=interp1(PRES(indnani:indnanf,p),DOXY(indnani:indnanf,p),Zaxis,'nearest','extrap');
% 
%     % Interpolate HR data ('RHO','CT','SA','N2','N', etc)
%     indnani=find(~isnan(D_2mb(:,p)),1,'first'); 
%     indnanf=find(~isnan(D_2mb(:,p)),1,'last');
% 
%     P2int(:,p)=interp1(D_2mb(indnani:indnanf,p),D_2mb(indnani:indnanf,p),Z2maxis,'nearest','extrap');
%     CTint(:,p)=interp1(D_2mb(indnani:indnanf,p),CT(indnani:indnanf,p),Z2maxis,'nearest','extrap');
%     SAint(:,p)=interp1(D_2mb(indnani:indnanf,p),SA(indnani:indnanf,p),Z2maxis,'nearest','extrap');
%     RHOint(:,p)=interp1(D_2mb(indnani:indnanf,p),RHO(indnani:indnanf,p),Z2maxis,'nearest','extrap');
%     N2int(:,p)=interp1(D_2mb(indnani:indnanf,p),N2(indnani:indnanf,p),Z2maxis,'nearest','extrap');
% end
% 
% % Calculate MLD
% for p=1:np
%     % Calculate MLD
%     MLD(p,3)=mld_m(RHOint(:,p),Z2maxis',0.125); % 
%     MLD(p,2)=mld_m(RHOint(:,p),Z2maxis',0.03); % 
%     MLD(p,1)=mld_m(RHOint(:,p),Z2maxis',0.005); % 
% end
% 
% %%  Example plot: binned data HR
% 
% Z=CTint; Zname='Temperature'; Zunits='[^{\circ}C]';
% Z=SAint; Zname='Salinity'; Zunits='[]';
% Z=RHOint-1000; Zname='Sigma'; Zunits='[mg m^{-3}]';
% Z=sqrt(N2int); Zname='N'; Zunits='[s^{-1}]';
% 
% [X,Y]=meshgrid(LocalTime,Z2maxis);
% 
% close all
% 
% figure(1)
% h=pcolor(X,Y,Z);
% h.EdgeColor = 'none';
% hold on
% axis ij
% 
% hc=colorbar; 
% set(get(hc,'Title'),'String',Zunits,'fontsize',20);
% cmin=floor(nanmin(Z(:))); cmax=floor(nanmax(Z(:)));
% cmin=nanmin(Z(:)); cmax=nanmax(Z(:));
% caxis([cmin cmax])
% if strcmp(Zname,'Nitrate')
%     cmap = cmocean('dense',25);
% elseif strcmp(Zname,'Oxygen')
%     cmap = cmocean('matter',25);
% elseif strcmp(Zname,'Oxygen')
%     cmap = cmocean('oxy',25);
% elseif strcmp(Zname,'Chla')
%     h=pcolor(X,Y,log10(Z));
%     h.EdgeColor = 'none';
%     % cmin=ceil(nanmin(log10(Z(:)))); cmax=ceil(nanmax(log10(Z(:))));
%     cmin=0; cmax=10;
%     cmap = cmocean('delta',25);
%     hc.Ticks = log10([cmin, 1, 3, 5, cmax]); % linspace(cmin, cmax, 1) ; 
%     hc.TickLabels = num2cell([cmin, 1, 3, 5, cmax]); %num2cell(10^cmin:10^cmax);
%     caxis(log10([cmin cmax]))
% elseif strcmp(Zname,'Sigma')
%     cmap = cmocean('deep',25);
%     caxis([24.8 25.8])
% elseif strcmp(Zname,'Temperature')
%     cmap=flipud(cbrewer('div','RdYlBu',25));
% elseif strcmp(Zname,'N')
%     cmap = m_colmap('jet','step',25);
% elseif strcmp(Zname,'w_ctd')
%     cmap = cmocean('speed',25);
% end
% colormap(cmap)
% 
% % % Overlay Density Contours
%     hold on
%     v=24:0.02:27;
%     v1=[floor(nanmin(RHOint(:)-1000)):0.05:26.7];
%     v2=[26.7:0.05:ceil(nanmax(RHOint(:)-1000))];
%     v=[v1 v2];
%     [C,h]=contour(X,Y,RHOint-1000,v); %contourf
%     % v=floor(nanmin(nanmin(Temp))):0.4:ceil(nanmax(nanmax(Temp)));
%     % [C,h]=contour(Time,Depth,Temp,v); %contourf
%     set(h,'color','k');
%     clabel(C,h,v,'FontSize',18,'LabelSpacing',1200); %,'FontSize',16,'LabelSpacing',1200); %,'manual);
%     
% %   % Overlay MLD (largest threshold)
%     hold on
%     plot(X(1,:),MLD(:,3),'--o','color',[0.9 0.9 0.9],'linewidth',2,'markersize',10,...
%         'markerfacecolor',[0.9 0.9 0.9],'markeredgecolor',[0.2 0.2 0.2]);
% 
% ylabel('Depth [m]');
% title(Zname)
% 
% xlim([nanmin(X(:)) nanmax(X(:))])
% datetick('keeplimits') % 'keepticks',
% day=datetime(datevec(LocalTime),'Format','MM/dd');
% xlabel([char(day(1)),' - ',char(day(end))]);
% % xlabel('Longitude [^{\circ} W]');
% set(gca,'xMinortick','on')
% % xticklabel_rotate([],45,[]);
% 
% print_setup(60,30)
% print([path_plots,Zname,'_CPF003_p',num2str(1),'-',num2str(np),'_HR.eps'],'-depsc2');
% 
% Xhr=X; Yhr=Y; % saving X,Y axis for contouring HR on top of LR BGC data below
% 
% %%  Example plot: binned data LR
% 
% Z=Nint; Zname='Nitrate'; Zunits='[\mumol/kg]';
% Z=DOXYint; Zname='Oxygen'; Zunits='[\mumol/kg]';
% % Z=Bint; Zname='B_{bp}'; Zunits='[m^{-1}]';
% Z=CHLAint; Zname='Chla'; Zunits='[mg/m^3]';
% 
% [X,Y]=meshgrid(LocalTime,Zaxis);
% 
% close all
% 
% figure(1)
% h=pcolor(X,Y,Z);
% h.EdgeColor = 'none';
% hold on
% axis ij
% 
% hc=colorbar; 
% set(get(hc,'Title'),'String',Zunits,'fontsize',20);
% cmin=floor(nanmin(Z(:))); cmax=floor(nanmax(Z(:)));
% caxis([cmin cmax])
% if strcmp(Zname,'Nitrate')
%     cmap = cmocean('dense',25);
% elseif strcmp(Zname,'Oxygen')
%     cmap = cmocean('matter',25);
% elseif strcmp(Zname,'Oxygen')
%     cmap = cmocean('oxy',25);
% elseif strcmp(Zname,'Chla')
%     h=pcolor(X,Y,log10(Z));
%     h.EdgeColor = 'none';
%     % cmin=ceil(nanmin(log10(Z(:)))); cmax=ceil(nanmax(log10(Z(:))));
%     cmin=0; cmax=10;
%     cmap = cmocean('delta',25);
%     hc.Ticks = log10([cmin, 1, 3, 5, cmax]); % linspace(cmin, cmax, 1) ; 
%     hc.TickLabels = num2cell([cmin, 1, 3, 5, cmax]); %num2cell(10^cmin:10^cmax);
%     caxis(log10([cmin cmax]))
% elseif strcmp(Zname,'w_ctd')
%     cmap = cmocean('speed',25);
% end
% colormap(cmap)
% 
% 
% %             % adjust colorbar
% %             c = cmocean('thermal',25);colormap(c)
% %             if strcmp(cind(j),'temp')==1
% %             c = cbrewer('div','RdYlBu',25);colormap(c)
% %             elseif strcmp(cind(j),'sal')==1
% %             c = cmocean('haline',25);colormap(c)
% %             elseif strcmp(cind(j),'sigt')==1
% %             c = cmocean('dense',25);colormap(c)
% %             elseif strcmp(cind(j),'oxy')==1
% %             c = cmocean('oxy',25);colormap(c)
% %             elseif strcmp(cind(j),'chla')==1
% %             c = cmocean('algae',25);colormap(c)
% %             elseif strcmp(cind(j),'poc')==1
% %             c = cmocean('delta',25);colormap(c)
% %             elseif strcmp(cind(j),'cdom')==1
% %             c = cmocean('matter',25);colormap(c)
% %             end
% 
% % % Overlay Density Contours
%     hold on
%     v=24:0.02:27;
%     v1=[floor(nanmin(RHOint(:)-1000)):0.05:26.7];
%     v2=[26.7:0.05:ceil(nanmax(RHOint(:)-1000))];
%     v=[v1 v2];
%     [C,h]=contour(Xhr,Yhr,RHOint-1000,v); %contourf
%     % v=floor(nanmin(nanmin(Temp))):0.4:ceil(nanmax(nanmax(Temp)));
%     % [C,h]=contour(Time,Depth,Temp,v); %contourf
%     set(h,'color','k');
%     clabel(C,h,v,'LabelSpacing',1200); %'FontSize',16,'LabelSpacing',1200); %,'manual);
%     
% %   % Overlay MLD (largest threshold)
%     hold on
%     plot(Xhr(1,:),MLD(:,3),'--o','color',[0.9 0.9 0.9],'linewidth',2,'markersize',10,...
%         'markerfacecolor',[0.9 0.9 0.9],'markeredgecolor',[0.2 0.2 0.2]);
% 
% ylabel('Depth [m]');
% title(Zname)
% 
% xlim([nanmin(X(:)) nanmax(X(:))])
% datetick('keeplimits') % 'keepticks',
% day=datetime(datevec(LocalTime),'Format','MM/dd');
% xlabel([char(day(1)),' - ',char(day(end))]);
% % xlabel('Longitude [^{\circ} W]');
% 
% print_setup(60,30)
% print([path_plots,Zname,'_CPF003_p',num2str(1),'-',num2str(np),'_LR.eps'],'-depsc2');
% 
% %% Define grid for BINNING data
% % NOTE: this 2-d binning will coarsen the resolution in th eX axis as well
% % % Define Variable to plot
% % Z=NITRATE; Zname='Nitrate'; Zunits='[\mu mol/kg]';
% % 
% % % Define X axis (i.e. Time or Lon or Lat)
% % [X,Y]=meshgrid(LocalTime,1:100)
% % yres=5;
% % northl=0;
% % southl=100; %
% % 
% % xres=0.05; % fraction of a day
% % westl=nanmin(X(:))-0.01;
% % eastl=nanmax(X(:))+0.01;
% % 
% % xin=X(:);
% % yin=PRES(:);
% % zin=Z(:);
% % 
% % [A,X,Y] = binning(xin,yin,zin,xres,yres,westl,eastl,southl,northl);
% % 
% % close all
% % 
% % figure(1)
% % h=pcolor(X,Y,A);
% % h.EdgeColor = 'none';
% % hold on
% % 
% % hc=colorbar; 
% % set(get(hc,'Title'),'String',Zunits,'fontsize',20);
% % 
% % if strcmp(Zname,'Nitrate')
% %     cmap = cmocean('dense',22);
% % elseif strcmp(Zname,'w_ctd')
% %     cmap = cmocean('speed',22);
% % end
% % colormap(cmap)
% % 
% % cmin=nanmin(Z(:)); cmax=nanmax(Z(:));
% % caxis([cmin cmax])
% % 
% % ylim([ceil(tsnanmode(Y(:))) 0]);
% % ylabel('Depth [m]');
% % title(Zname)
% % 
% % datetick('keeplimits')
% % xlabel('Time');
% % % xlabel('Longitude [^{\circ} W]');
% % 
% % print_setup(40,30)
