% plot time series surface ageostrophic V and wind stress

clear
clc

nx = 2160;
% subtropic size = [72 81 27]
% subpolar size = [72 109 27]
% southern size = [72 119 27];
region_name = 'DopplerScat/southern/';
siz1 = [72 119 1];
sizstr = '72x119';
sizstr27 = '72x119x27';
siz27 = [72 119 27];
lats = (readbin(['~dmenemen/llc_2160/regions/' region_name 'grid/YC_' sizstr],siz1));
lons = readbin(['~dmenemen/llc_2160/regions/' region_name 'grid/XC_' sizstr],siz1);
landV = readbin(['~dmenemen/llc_2160/regions/' region_name 'grid/hFacS_' sizstr27],siz1);
landU = readbin(['~dmenemen/llc_2160/regions/' region_name 'grid/hFacW_' sizstr27],siz1);
dy = readbin(['~dmenemen/llc_2160/regions/' region_name 'grid/DYC_' sizstr],siz1);
dx = readbin(['~dmenemen/llc_2160/regions/' region_name 'grid/DXC_' sizstr],siz1);
z = readbin('~dmenemen/llc_2160/grid/RC.data',[27]);

fc = 2*7.2921*10^-5*sind(lats);
g = 9.806;
n = 0;
for ts = 92160:(24*80):1586400;%92160+90*80*24
    
    dyy=ts2dte(ts,45,2011,1,17,30);
    
    fnamTy = ['~dmenemen/llc_2160/regions/' region_name 'oceTAUY/oceTAUY' '_' sizstr '.' dyy];
    fnamTx = ['~dmenemen/llc_2160/regions/' region_name 'oceTAUX/oceTAUX' '_' sizstr '.' dyy];
    
    fnamEta = ['~dmenemen/llc_2160/regions/' region_name 'Eta/Eta' '_' sizstr 'x1.' dyy];
    
    fnamU = ['~dmenemen/llc_2160/regions/' region_name 'U/U_' sizstr 'x1' '.' dyy];
    fnamV = ['~dmenemen/llc_2160/regions/' region_name 'V/V_' sizstr 'x1' '.' dyy];
    
    if exist(fnamU) && exist(fnamV) && exist(fnamTy) && exist(fnamTx)
        
        eta = readbin(fnamEta,siz1);
        
        if strcmp(region_name,'DopplerScat/southern/')
            Umodel = readbin(fnamU,siz1);
            Vmodel = readbin(fnamV,siz1);
            dyG = readbin(['~dmenemen/llc_2160/regions/' region_name 'grid/DYG_' sizstr],siz1);
            dxG = readbin(['~dmenemen/llc_2160/regions/' region_name 'grid/DXG_' sizstr],siz1);
            Ty = readbin(fnamTy,siz1);
            Tx = readbin(fnamTx,siz1);
        else
            Umodel = readbin(fnamV,siz1);
            Vmodel = -circshift(readbin(fnamU,siz1),[0 -1]);
            Vmodel(:,end) = nan;
            Ty = -circshift(readbin(fnamTx,siz1),[0 -1]);
            Tx= readbin(fnamTy,siz1);
            dyG = readbin(['~dmenemen/llc_2160/regions/' region_name 'grid/DXG_' sizstr],siz1);
            dxG = readbin(['~dmenemen/llc_2160/regions/' region_name 'grid/DYG_' sizstr],siz1);
            dxG = circshift(dxG,[0 -1]);
            dxG(:,end) = nan;
        end
        
        if sum(sum(Tx)) ~=0 && sum(sum(Ty)) ~=0 && sum(sum(Umodel)) ~=0 && sum(sum(Vmodel)) ~=0
            n = n + 1;
            
            %pin = ['/nobackupp9/dmenemen/llc_2160/regions/' region_name '/averages_winter/'];
            %%%%%% NEED A FOR LOOP HERE TO LOOP THROUGH FILES
            
            s = size(eta);
            
            deltaEtaX = [diff(eta,1,1);eta(1,:)-eta(end,:)];
            deltaEtaY = [eta(:,2:end) - eta(:,1:end-1),nan*ones(s(1),1)];
            
            ffY = [((fc(:,1:end-1) + fc(:,2:end))/2),nan*ones(s(1),1)];
            ffX = [(fc(1:end-1,:) + fc(2:end,:)) / 2;(fc(end,:)+fc(1,:))/2];
            
            Us = -(g ./ ffY .* deltaEtaY ./ dy);
            Vs = g ./ ffX .* deltaEtaX ./ dx;
            
            Us =  Us.*landV;
            Vs = Vs.*landU;     % switched b/c edges are switched right now
            
            %% interpolate to Umit grid points
            s= size(Us);
            Us1 = (Us(:,1:end-1) + Us(:,2:end))/2;
            Us1(:,end+1) = nan*ones(s(1),1);
            
            Umit = Us1(1:end-1,:) + (.5.*dxG(1:end-1,:)./dx(2:end,:)).*(Us1(2:end,:)-Us1(1:end-1,:));
            s = size(Umit);
            Umit(end+1,:) = Us1(end,:) + (.5.*dxG(end,:)./dx(1,:)).*(Us1(1,:)-Us1(end,:));
            
            Vs1 = (Vs(1:end-1,:) + Vs(2:end,:))/2;
            s = size(Vs1);
            Vs1(end+1,:) = (Vs(1,:)+Vs(end,:))/2;
            
            Vmit = Vs1(:,1:end-1) + (.5.*dyG(:,1:end-1)./dy(:,2:end)).*(Vs1(:,2:end)-Vs1(:,1:end-1));
            s = size(Vmit);
            Vmit(:,end+1) = zeros(s(1),1);
            
            %% Compute ekman
            
            Uek = Umodel-Umit;
            Vek = Vmodel-Vmit;
            
            UVek = sqrt(Uek.^2 + Vek.^2);
            UVekm(n) = mean(mean(UVek(~isnan(UVek)&abs(UVek)~=inf)));
            
            TXY = sqrt(sqrt(Tx.^2 + Ty.^2));
            TXYm(n) = mean(mean(TXY(~isnan(TXY)&abs(TXY)~=inf)));
            
        end
    end
end

%% Plot


figure('color','white','units','normalized','outerposition',[0 0 1 1],'name','1');
hh = tight_subplot(1,2);

axes(hh(1))
plot(UVekm);
title('Ageostrophic Velocity')
%daspect([1 1 1])
set(get(gca,'title'),'fontsize',16);

axes(hh(2))
plot(TXYm);
title('Surface Shear')
%daspect([1 1 1])
set(get(gca,'title'),'fontsize',16);

hf = hamming(7);
UVf = filter(hf/7,1,UVekm);
Tf = filter(hf/7,1,TXYm);

figure('color','white','units','normalized','outerposition',[0 0 1 1],'name','1');
[ax,h1,h2] = plotyy(10:375,Tf(10:375),10:375,UVf(10:375));
set(h1,'marker','.','color','b','linewidth',2,'markersize',2);
set(h2,'marker','.','linewidth',2,'markersize',2);
legend('Sqrt Surface Shear','Ageostrophic Speed (m/s)')
%title('Subpolar (45N, 180W) Ageostrophic Speed and Surface Shear March 2011-March 2012, 7d LPF')
title('Southern (50S, 90E) Ageostrophic Speed and Surface Shear March 2011-March 2012, 7d LPF')
%title('Subtropic (15N, 130W) Ageostrophic Speed and Surface Shear March 2011-March 2012, 7d LPF')
set(get(gca,'title'),'fontsize',16);


 set(ax,'XTick',[0:365/12:365])
 set(ax,'XTickLabel', ['Mar';'Apr';'May';'Jun';'Jul';'Aug';'Sep';'Oct';'Nov';'Dec';'Jan';'Feb';'Mar']);
 set(ax,'FontSize',14);
 set(ax(2),'YLim',[-.1 .25]);
 set(ax(1),'YLim',[.15 .5]);
 set(get(ax(1),'xlabel'),'FontSize',14);
 set(get(ax(1),'ylabel'),'fontsize',16,'String','Sqrt Surface Shear');
 set(get(ax(2),'ylabel'),'fontsize',16,'String','Ageostrophic Speed (m/s)');
 set(get(ax(2),'title'),'fontsize',16);
 
print('-djpeg','-r300','southern_ageo_sqrtstress_year.jpeg')
 R = corrcoef([UVf', Tf'])