`% Seasonal maps of angle between Ekman surface velocity and surface wind
% stress for four wind stress regimes

clear
clc

nx = 2160;
siz = [8640 7560 1];



lats = quikreadpcolor_llc('~dmenemen/llc_2160/grid/YC.data',nx);
lons = quikreadpcolor_llc('~dmenemen/llc_2160/grid/XC.data',nx);

n = 0;
Txs = 0;
Tys = 0;
Uds = 0;
Vds = 0;
strNames = {};
TXY = 0;
angUV = 0;
angT = 0;
angT1t = 0;
angT2t = 0;
angT3t = 0;
angT4t = 0;

angUV1t = 0;
angUV2t = 0;
angUV3t = 0;
angUV4t = 0;

Utek = 0;
Vtek = 0;
UekT = 0;
VekT = 0;
TxT = 0;
TyT = 0;
n14 = 0;
TyT1t = 0;
TyT2t = 0;
TyT3t = 0;
TyT4t = 0;

TxT1t = 0;
TxT2t = 0;
TxT3t = 0;
TxT4t =0;

UekT1t = 0;
UekT2t =0;
UekT3t =0;
UekT4t = 0;

VekT1t = 0;
VekT2t =0;
VekT3t = 0;
VekT4t =  0;

for ts = 92160:(24*80):92160+90*80*24%:1586400;
  
    dy=ts2dte(ts,45,2011,1,17,30);
    fld=0; nn=0;
    
   
    fnamTy = ['~dmenemen/llc_2160/regions/global/oceTAUY_daily/oceTAUY' '_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
    fnamTx = ['~dmenemen/llc_2160/regions/global/oceTAUX_daily/oceTAUX' '_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
    
    fnamUg = ['~dmenemen/llc_2160/regions/global/Ug_daily/Ug_' int2str((siz(1))) 'x' int2str((siz(2))) '.' dy];
    fnamVg = ['~dmenemen/llc_2160/regions/global/Vg_daily/Vg_' int2str((siz(1))) 'x' int2str((siz(2))) '.' dy];
     
    fnamU = ['~dmenemen/llc_2160/regions/global/U_daily/U_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
    fnamV = ['~dmenemen/llc_2160/regions/global/V_daily/V_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
 
    if exist(fnamUg) && exist(fnamVg) && exist(fnamTy) && exist(fnamTx)
        
        Ug = readbin(fnamUg,siz);
        Vg = readbin(fnamVg,siz);
        [Tx, Ty] = quikreadpcolor_uv_llc(fnamTx,fnamTy,nx);
        [U, V] = quikreadpcolor_uv_llc(fnamU,fnamV,nx);
        
        if sum(sum(Tx)) ~=0 && sum(sum(Ty)) ~=0 && sum(sum(Ug))~=0 && sum(sum(U)) ~=0 && sum(sum(V)) ~=0 && sum(sum(Vg)) ~=0
            n = n + 1;
            strNames{n} = dy;
            
            Vek = V-Vg;
            Uek = U-Ug;
            %TXY = TXY + sqrt(Tx.^2 + Ty.^2);
            
            UekT = UekT + Uek;
            VekT = VekT + Vek;
            TxT = TxT + Tx;
            TyT = TyT + Ty;

            if mod(n,14) == 0
                n14 = n14 + 1
                %TXY = TXY./14;
                TxT = TxT./14;
                TyT = TyT./14;
                UekT = UekT./14;
                VekT = VekT./14;
                TXY = sqrt(TxT.^2 + TyT.^2);
                
                weights = ones(111)./(111*111);
                VekT = conv2(VekT,weights);
                UekT= conv2(UekT,weights);
                TxT = conv2(TxT,weights);
                TyT = conv2(TyT,weights);
                TXY = conv2(TXY,weights);
                
                TyT1 = TyT; TyT2 = TyT; TyT3 = TyT; TyT4 = TyT;
                TxT1 = TxT; TxT2 = TxT; TxT3 = TxT; TxT4 = TxT;
                UekT1 = UekT; UekT2 = UekT; UekT3 = UekT; UekT4 = UekT;
                VekT1 = VekT; VekT2 = VekT; VekT3 = VekT; VekT4 = VekT;
                
                TyT1(TXY<.0099|TXY>.0397) = 1i;
                TyT2(TXY<.0397|TXY>.0892) = 1i;
                TyT3(TXY<.0892|TXY>.1586) = 1i;
                TyT4(TXY<.1586|TXY>.2478) = 1i;
                
                TxT1(TXY<.0099|TXY>.0397) = 1i;
                TxT2(TXY<.0397|TXY>.0892) = 1i;
                TxT3(TXY<.0892|TXY>.1586) = 1i;
                TxT4(TXY<.1586|TXY>.2478) = 1i;
               
                UekT1(TXY<.0099|TXY>.0397) = 1i;
                UekT2(TXY<.0397|TXY>.0892) = 1i;
                UekT3(TXY<.0892|TXY>.1586) = 1i;
                UekT4(TXY<.1586|TXY>.2478) = 1i;
                
                VekT1(TXY<.0099|TXY>.0397) = 1i;
                VekT2(TXY<.0397|TXY>.0892) = 1i;
                VekT3(TXY<.0892|TXY>.1586) = 1i;
                VekT4(TXY<.1586|TXY>.2478) = 1i;
                
                TyT1t = TyT1t + TyT1;
                TyT2t = TyT2t + TyT2;
                TyT3t = TyT3t + TyT3;
                TyT4t = TyT4t + TyT4;
                
                TxT1t = TxT1t + TxT1;
                TxT2t = TxT2t + TxT2;
                TxT3t = TxT3t + TxT3;
                TxT4t = TxT4t + TxT4;
               
                UekT1t = UekT1t + UekT1;
                UekT2t = UekT2t + UekT2;
                UekT3t = UekT3t + UekT3;
                UekT4t = UekT4t + UekT4;
                
                VekT1t = VekT1t + VekT1;
                VekT2t = VekT2t + VekT2;
                VekT3t = VekT3t + VekT3;
                VekT4t = VekT4t + VekT4;
                
 %{               
%                 TxT = colfilt(TxT,[111,111],'sliding',@mean);
%                 TXY = colfilt(TXY,[111,111],'sliding',@mean);
%                 TyT = colfilt(TyT,[111,111],'sliding',@mean);
%                 UekT = colfilt(UekT,[111,111],'sliding',@mean);
%                 VekT = colfilt(VekT,[111,111],'sliding',@mean);

%                 for ii = 1 : siz(1)
%                     for jj = 1 : siz(2)
                        
%                         TxT(ii,jj) =  mean(mean(TxT(~isnan(TxT)&TxT~=0&abs(TxT)~=inf&abs(lats-lats(ii,jj))<2.2&abs(lons-lons(ii,jj))<2.2)));
%                         TXY(ii,jj) = mean(mean(TXY(~isnan(TXY)&TXY~=0&abs(TXY)~=inf&abs(lats-lats(ii,jj))<2.2&abs(lons-lons(ii,jj))<2.2)));
%                         TyT(ii,jj) =  mean(mean(TyT(~isnan(TyT)&TyT~=0&abs(TyT)~=inf&abs(lats-lats(ii,jj))<2.2&abs(lons-lons(ii,jj))<2.2)));
%                         UekT(ii,jj) = mean(mean(UekT(~isnan(UekT)&UekT~=0&abs(UekT)~=inf&abs(lats-lats(ii,jj))<2.2&abs(lons-lons(ii,jj))<2.2)));
%                         VekT(ii,jj) = mean(mean(VekT(~isnan(VekT)&VekT~=0&abs(VekT)~=inf&abs(lats-lats(ii,jj))<2.2&abs(lons-lons(ii,jj))<2.2)));
% 
%                     end
%                 end
%                 
 %}               
                
%                 angT = atan2(TyT,TxT)/pi*180;
%                 angUV = atan2(VekT,UekT)/pi*180;
%                 
%                 angT1 = angT;
%                 angT2 = angT;
%                 angT3 = angT;
%                 angT4 = angT;
%                 
%                 angUV1 = angUV;
%                 angUV2 = angUV;
%                 angUV3 = angUV;
%                 angUV4 = angUV;
%                   
%                 angT1(TXY<.0099|TXY>.0397) = 1i;
%                 angT2(TXY<.0397|TXY>.0892) = 1i;
%                 angT3(TXY<.0892|TXY>.1586) = 1i;
%                 angT4(TXY<.1586|TXY>.2478) = 1i;
%                 
%                 angUV1(TXY<.0099|TXY>.0397) = 1i;
%                 angUV2(TXY<.0397|TXY>.0892) = 1i;
%                 angUV3(TXY<.0892|TXY>.1586) = 1i;
%                 angUV4(TXY<.1586|TXY>.2478) = 1i;
% 
%                 angT1t = angT1t + angT1;
%                 angT2t = angT2t + angT2;
%                 angT3t = angT3t + angT3;
%                 angT4t = angT4t + angT4;
%                 
%                 angUV1t = angUV1t + angUV1;
%                 angUV2t = angUV2t + angUV2;
%                 angUV3t = angUV3t + angUV3;
%                 angUV4t = angUV4t + angUV4;
               
                TXY = 0;
                angT = 0;
                angUV = 0;
                TxT = 0;
                TyT = 0;
                VekT = 0;
                UekT = 0;
                
            end
%             Utek = Utek + Uek;
%             Vtek = Vtek + Vek;
%             Txt = Txt + Tx;
%             Tyt = Tyt + Ty;
%           
        end
    end
end

clear Uek Vek TXY U Ug V Vg angT1 angT UVd angT2 angT3 angT4 angUV1 angUV2 angUV3 angUV4 angUV


nn1 = n14.*ones(8640+100,7560+100) - imag(TxT1t);
nn2 = n14.*ones(8640+100,7560+100) - imag(TxT2t);
nn3 = n14.*ones(8640+100,7560+100) - imag(TxT3t);
nn4 = n14.*ones(8640+100,7560+100) - imag(TxT4t);

angT1t = atan2(real(TyT1t)./nn1,real(TxT1t)./nn1)*180/pi;
angT2t = atan2(real(TyT2t)./nn2,real(TxT2t)./nn2)*180/pi;
angT3t = atan2(real(TyT3t)./nn3,real(TxT3t)./nn3)*180/pi;
angT4t = atan2(real(TyT4t)./nn4,real(TxT4t)./nn4)*180/pi;

angUV1t = atan2(real(Vek1t)./nn1,real(Uek1t)./nn1)*180/pi;
angUV2t = atan2(real(Vek2t)./nn2,real(Uek2t)./nn2)*180/pi;
angUV3t = atan2(real(Vek3t)./nn3,real(Uek3t)./nn3)*180/pi;
angUV4t = atan2(real(Vek4t)./nn4,real(Uek4t)./nn4)*180/pi;

angT1t = angT1t(51:8640+50,51:7560+50);
angT2t = angT2t(51:8640+50,51:7560+50);
angT3t = angT3t(51:8640+50,51:7560+50);
angT4t = angT4t(51:8640+50,51:7560+50);
angUV1t = angUV1t(51:8640+50,51:7560+50);
angUV2t = angUV2t(51:8640+50,51:7560+50);
angUV3t = angUV3t(51:8640+50,51:7560+50);
angUV4t = angUV4t(51:8640+50,51:7560+50);

angd1 = nan*ones(8640,7560);
angd2 = angd1; angd3 = angd1; angd4 = angd1;

angd1 = (angT1t-angUV1t);
angd2 = (angT2t-angUV2t);
angd3 = (angT3t-angUV3t);
angd4 = (angT4t-angUV4t);


% angd1(angT1t>angUV1t) = ((angT1t(real(angT1t)>real(angUV1t)) - angUV1t(real(angT1t)>real(angUV1t)))./nn1(angT1t>angUV1t));
% angd2(angT2t>angUV2t) = ((angT2t(real(angT2t)>real(angUV2t)) - angUV2t(real(angT2t)>real(angUV2t)))./nn2(angT2t>angUV2t));
% angd3(angT3t>angUV3t) = ((angT3t(real(angT3t)>real(angUV3t)) - angUV3t(real(angT3t)>real(angUV3t)))./nn3(angT3t>angUV3t));
% angd4(angT4t>angUV4t) = ((angT4t(real(angT4t)>real(angUV4t)) - angUV4t(real(angT4t)>real(angUV4t)))./nn4(angT4t>angUV4t));
% 
% angd1(angT1t<angUV1t) = -real((-angT1t(real(angT1t)<real(angUV1t)) + angUV1t(real(angT1t)<real(angUV1t)))./nn1(angT1t<angUV1t));
% angd2(angT2t<angUV2t) = -real((-angT2t(real(angT2t)<real(angUV2t)) + angUV2t(real(angT2t)<real(angUV2t)))./nn2(angT2t<angUV2t));
% angd3(angT3t<angUV3t) = -real((-angT3t(real(angT3t)<real(angUV3t)) + angUV3t(real(angT3t)<real(angUV3t)))./nn3(angT3t<angUV3t));
% angd4(angT4t<angUV4t) = -real((-angT4t(real(angT4t)<real(angUV4t)) + angUV4t(real(angT4t)<real(angUV4t)))./nn4(angT4t<angUV4t));

angd1(angd1==0) = nan;
angd2(angd2==0) = nan;
angd3(angd3==0) = nan;
angd4(angd4==0) = nan;

angd1(abs(lats)<2) = nan;
angd2(abs(lats)<2) = nan;
angd3(abs(lats)<2) = nan;
angd4(abs(lats)<2) = nan;

        
%% plot



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

axes(hh(1))
quikpcolor(angd1')
colorbar
caxis([-90 90])
title('2.5-5 m/s wind')
daspect([1 1 1])
%set(get(gca,'title'),'fontsize',16);
set(gca,'xtick',[])
set(gca,'ytick',[])
axis([0 8640 1500 6500])
axes(hh(2))
quikpcolor(angd2');
colorbar
caxis([-90 90])
title('5-7.5 m/s wind')
daspect([1 1 1])
%set(get(gca,'title'),'fontsize',16);
set(gca,'xtick',[])
set(gca,'ytick',[])
axis([0 8640 1500 6500])
axes(hh(3))
quikpcolor(angd3');
colorbar
caxis([-90 90])
title('7.5-10 m/s wind')
daspect([1 1 1])
%set(get(gca,'title'),'fontsize',16);
set(gca,'xtick',[])
set(gca,'ytick',[])
axis([0 8640 1500 6500])
axes(hh(4))
quikpcolor(angd4');
colorbar
caxis([-90 90])
title('10-12.5 m/s wind')
daspect([1 1 1])
%set(get(gca,'title'),'fontsize',16);
set(gca,'xtick',[])
set(gca,'ytick',[])
axis([0 8640 1500 6500])
colormap([0 0 0 ; redblue(128)]);