%Load in the Ocean T data
%Get info out of netcdf
file = 'thetao_Omon_GFDL-ESM2M_piControl_r1i1p1_049601-050012.nc';
info = ncinfo(file);

%Extract T, lat, lon, lev, lat_vertices, lon_vertices, data from netcdf
T = ncread(file,'thetao',[1 1 1 1], [inf inf inf inf]); %want all levels of data
LAT = ncread(file,'lat');
LON = ncread(file,'lon');
LON(LON<0) = LON(LON<0)+360;
LEV = ncread(file,'lev');
lat_vert = ncread(file,'lat_vertices');
lon_vert = ncread(file,'lon_vertices');

%use area of each grid box to calculate the weight it should have for
%taking an area weighted average
weights = areaquad(squeeze(lat_vert(1,:,:)),squeeze(lon_vert(1,:,:)),...
    squeeze(lat_vert(3,:,:)),squeeze(lon_vert(2,:,:)));
weights = weights/(sum(sum(weights)));

%calculate Nino3 (5S-5N, 90W-150W) signal over time
%keep weights only in Nino3 Region and normalize
weights_nino3 = weights;
weights_nino3(LAT>5) = NaN;
weights_nino3(LAT<-5) = NaN;
weights_nino3(LON<210) = NaN;
weights_nino3(LON>270) = NaN;
weights_nino3 = weights_nino3/nansum(nansum(weights_nino3));
%weighted average of SST in Nino3 region over time
Nino3 = squeeze(nansum(nansum(squeeze(T(:,:,1,:)).*repmat(weights_nino3,[1 1 size(T,4)]))));

%Remove Seasonal Cycle
for i=1:12
    Ninoclim(:,i) = nanmean(Nino3);
end
Nino3 = Nino3-repmat(Ninoclim', [length(Nino3)/12 1]);

%Interpolate depth of thermocline (13C) over Equatorial Pacific: 130E-280E, 15S-15N
%Keep only temperatures in the Equatorial Pacific Region
T(:,abs(LAT(1,:))>15,:,:) = [];
T(LON(:,1)<130,:,:,:) = [];
lon = LON(:,1); lon(lon<130) = [];
T(lon>280,:,:,:) = [];
lon(lon>280) = [];
lat = LAT(1,:); lat(abs(lat)>15) = [];
%itereate through each location at each time to interpolate the depth of
%the thermocline (Note, there is probably a much faster way to do this)
DepthThermo = nan(size(T,1), size(T,2), size(T,4));
for i = 1:size(T,1)
    for j = 1:size(T,2)
        for k=1:size(T,4)
            try
                level = LEV;
                Tprofile = squeeze(T(i,j,:,k));
                level(isnan(Tprofile)) = [];
                Tprofile(isnan(Tprofile)) = [];
                DepthThermo(i,j,k) = interp1(Tprofile,level,13+273.15);
            end
        end
    end
end
            
%Remove the seasonal cycle
for i=1:12
    DepthThermoClim(:,:,i) = nanmean(DepthThermo,3);
end
DepthThermo = DepthThermo-repmat(DepthThermoClim,[1 1 size(DepthThermo,3)/12]);

%Calculate regression of Nino3 on thermocline depth at each location
Reg = nan(size(DepthThermo,1), size(DepthThermo,2));
for i=1:size(DepthThermo,1)
    for j=1:size(DepthThermo,2)
        dt = squeeze(DepthThermo(i,j,:));
        Reg(i,j) = regress(dt(:),Nino3(:));
    end
end

Reg(Reg==0) = NaN;

%latitude and longitude for mapping data
mlat = lat;
mlon = lon;

%plot the Nino3 regression on the Thermocline depth over time (slider
%controls plot time index)
figure; load coast;
axesm('miller','MapLonLimit',[110 300], 'MapLatLimit',[-20 20]);
plotm(lat,long,'k'); hold on;
title('Nino3 Regression on Thermocline Depth');
tightmap;
hplot = pcolorm(double(mlat),double(mlon),Reg*Nino3(1)); colormap(b2r_noWhite(-15,15));
h = uicontrol('style','slider','value',1,'max',length(Nino3),'min',1,'sliderstep',[1 1]/(length(Nino3)-1));
set(h,'Callback', @(hObject,event) updateplot(hObject,event,hplot,Nino3,Reg));

