clear all
close all
addpath(genpath('/Users/Mau/Documents/MATLABhelp_new'));

%% Purpose: Construct seasonality of MLD from CTD data at M1
% This code reads a bunch of files (~1/year); and constrcuts a climatology 

bid='M1';
% path_data=['/Users/Mau/Documents/MBARIprojects/M1/'];
% path_data=['http://dods.mbari.org:80/opendap/data/ssdsdata/deployments/m1/202008/OS_M1_20200825tenMinute_CMSTV.nc'];
% NC=ncinfo(([path_data]));

path_data=['/Users/Mau/Documents/MBARIprojects/M1/'];
files=dir([path_data,'*.nc4']);
filename=files(1).name;
NC=ncinfo(([path_data,filename]));
nf=length(files);

path_plots='/Users/Mau/Documents/MBARIprojects/M1/plots/';

%% DATE=cell(nf,1);
Sal=cell(nf,1);
Temp=cell(nf,1);
Depth=cell(nf,1);

LON=cell(nf,1);
LAT=cell(nf,1);
DATE=cell(nf,1);
MONTH=cell(nf,1);
DOY=cell(nf,1);

for f=1:nf
    
    filename=files(f).name;
    
    T=squeeze(ncread([path_data,filename],'SEA_WATER_SALINITY_TM')); % readtable(filename);
    S=squeeze(ncread([path_data,filename],'SEA_WATER_TEMPERATURE_TM'));
    
    if f==1
        D=squeeze(ncread([path_data,filename],'DEPTH1'));
        lon=ncread([path_data,filename],'LONGITUDE1');
    else
        D=squeeze(ncread([path_data,filename],'DEPTH'));
        lon=ncread([path_data,filename],'LONGITUDE');
    end
    
    
    lat=ncread([path_data,filename],'LATITUDE'); % 'GPS_LATITUDE_TM');
    
    
    % Load time and convert to matlab time
    time_ctd=ncread([path_data,filename],'tm_time_ts'); % gregorian, units: "seconds since 1970-01-01 00:00:00";
    to=datenum([1970 1 1 0 0 0]);
    time=time_ctd./86400; %in fracion of day
    mtime=time+to; % adjust start time
    mDATE=datevec(mtime); % OK, data every ~10'
    month=mDATE(:,2);
    Date=datetime(mDATE);
    doy=day(Date,'dayofyear');
    
    missing_value = -1.0E34;
    T(T==missing_value)=NaN;
    S(S==missing_value)=NaN;
    
    Temp{f}=T;
    Sal{f}=S;
    Depth{f}=D;
    
    LON{f}=lon;
    LAT{f}=lat;
    DATE{f}=mDATE;
    MONTH{f}=month;
    DOY{f}=doy;
    
    % Read ADCP data
    % S=ncread([path_data,filename],'eastward_sea_water_velocity_TM');
    
    % Save in cell arrays to process all deployments
    % [nt,nc]=size(D.data);
    % D.data(D.data==999)=NaN;
    % DATE{f}=zeros(nt,6);
    % DATE{f}(:,1:5)=D.data(:,1:5);
    % WSPD{f}=D.data(:,7); % in m/s
    % WDIR{f}=D.data(:,12);% in deg, met convention (360,0 wind coming from the N or southerly wind)
    
end

%% Compute N2, density and MLD

Pres=cell(nf,1);
SA=cell(nf,1);
CT=cell(nf,1);
RHO=cell(nf,1);

N2=cell(nf,1);
p_mid=cell(nf,1);

MLD_0125=cell(nf,1);
MLD_003=cell(nf,1);
MLD_0005=cell(nf,1);
mN2=cell(nf,1);
indz_mld=cell(nf,1);

nz=nan(nf);
np=nan(nf);
%%
for f=1:nf
    [nz(f),np(f)]=size(Temp{f});
    
    for p=1:np(f)

%         % Compute CT, SA and density
%         Pres{f}= gsw_p_from_z(-1*Depth{f},LAT{f});
%         [SA{f}(:,p), in_ocean] = gsw_SA_from_SP(Sal{f}(:,p),Pres{f},LON{f},LAT{f});        % Absolute Salinity from Pactical salinity
%         CT{f}(:,p)=gsw_CT_from_t(SA{f}(:,p),Temp{f}(:,p),Pres{f});                         % Conservative Temperature
%         RHO{f}(:,p)= gsw_rho(SA{f}(:,p),CT{f}(:,p),0);   % Pres{f}(:,p)                         % In-situ Density, use p_ref=0 to calculate potential density

        % Compute N2
        [N2{f}(:,p), p_mid{f}] = gsw_Nsquared(SA{f}(:,p),CT{f}(:,p),Pres{f},LAT{f});
         mN2{f}(p,1) = nanmedian(N2{f}(:,p));

%         % Calculate MLD
%         MLD_0125{f}(p,1)=mld_m(RHO{f}(:,p),Depth{f},0.125); % 0.125 density threshold(W_ctd{f}(1:100,:),0,1)';
%         MLD_003{f}(p,1)=mld_m(RHO{f}(:,p),Depth{f},0.03); % 0.125 density threshold(W_ctd{f}(1:100,:),0,1)';
%         MLD_0005{f}(p,1)=mld_m(RHO{f}(:,p),Depth{f},0.005); % 0.125 density threshold(W_ctd{f}(1:100,:),0,1)';
    end
end


%% 
save('M1_data_2016-2020.mat');
%% Plot PDF by month

pdffit2(cell2mat(MLD_0125),1,'k');


%%  Seasonality for MLD

y=cell2mat(MLD_0125);
month=cell2mat(MONTH);
doy=cell2mat(DOY);

bin=doy;
nb=length(unique(bin));

Nm=nan(nb,1);
my=Nm;
my_error=Nm;

for m=1:nb
    indm=find(bin==m);
    Nm(m)=nansum(~isnan(y(indm)));
    
    my(m)=nanmean(y(indm));
    my_error(m)=2*nanstd(y(indm),0,1); %./sqrt(Nm(m));
    
%     my(m)=nanmedian(y(indm));
%     my_error(m)=1.253.*nanstd(y(indm),0,1)./sqrt(Nm(m)); % error for the median
end

%% Plot MLD seasonal cycle and std

path_plots='/Users/Mau/Documents/MBARIprojects/M1/';

close all
errorbar([1:nb],my,my_error,'k'); hold on
xlim([0 nb])
ylim([0 80])

% xticks([1:12])
% xlabel(['month']);
datetick

ylabel(['MLD [m]']);
title('M1')

print_setup
maximize;
print([path_plots,'Seasonality_MLD0125_doy.eps'],'-depsc2')    

%% Max MLD

y=cell2mat(MLD_0125);
month=cell2mat(MONTH);
doy=cell2mat(DOY);

bin=doy;
nb=length(unique(bin));

Nm=nan(nb,1);
my=Nm;
my_max=Nm;
my_error=Nm;

for m=1:nb
    indm=find(bin==m);
    Nm(m)=nansum(~isnan(y(indm)));
    
    my(m)=nanmax(y(indm));
    my_max(m)=nanmax(y(indm));
%    my_error(m)=2*nanstd(nanmax(y(indm)),0,1)./sqrt(Nm(m));
    
    my(m)=nanmedian(y(indm));
    my_error(m)=1.253.*nanstd(y(indm),0,1)./sqrt(Nm(m)); % error for the median
end


close all
plot([1:nb],my,'k'); hold on
plot([1:nb],my_max,'r'); hold on
errorbar([1:nb],my,my_error,'k'); hold on
xlim([0 nb])
ylim([0 80])

% xticks([1:12])
% xlabel(['month']);
datetick

ylabel(['MLD [m]']);
title('M1')

print_setup(40,30)
maximize;
print([path_plots,'Seasonality_MLD0125_median_doy.eps'],'-depsc2')    


%% N2

y=cell2mat(mN2);
month=cell2mat(MONTH);
doy=cell2mat(DOY);

bin=doy;
nb=length(unique(bin));

Nm=nan(nb,1);
my=Nm;
my_min=Nm;
my_max=Nm;
my_error=Nm;

for m=1:nb
    indm=find(bin==m);
    Nm(m)=nansum(~isnan(y(indm)));
    
    my_min(m)=nanmin(y(indm));
    my_max(m)=nanmax(y(indm));
    
    my(m)=nanmedian(y(indm));
    my_error(m)=1.253.*nanstd(y(indm),0,1)./sqrt(Nm(m)); % error for the median
end


close all
plot([1:nb],my,'k'); hold on
% errorbar([1:nb],my,my_error,'k'); hold on

plot([1:nb],my_min,'b'); hold on
plot([1:nb],my_max,'r'); hold on

xlim([0 nb])
% ylim([0 80])

% xticks([1:12])
% xlabel(['month']);
datetick

ylabel(['N^2 [1/s]']);
title('M1')

legend('median','min','max','location','best')

print_setup(40,30)
maximize;
print([path_plots,'Seasonality_N2_median_doy.eps'],'-depsc2')    