% 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;

Utek_l = 0;
Vtek_l =0;
Txt_l =0;
Tyt_l =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); % around 4.5 deg box
                VekT = conv2(VekT,weights);
                UekT= conv2(UekT,weights);
                TxT = conv2(TxT,weights);
                TyT = conv2(TyT,weights);
                TXY = conv2(TXY,weights);

                angT = atan2(TyT,TxT)/pi*180;
                angUV = atan2(VekT,UekT)/pi*180;
                
                angUV(angUV<0) = angUV(angUV<0)  + 360;
                angT(angT<0) = angT(angT<0) + 360;

                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_l = Utek_l + Uek;
            Vtek_l = Vtek_l + Vek;
            Txt_l = Txt_l + Tx;
            Tyt_l = Tyt_l + Ty;
%           
        end
    else
        disp('Does not exist');
    end
end

clear Uek Vek TXY U Ug V Vg angT1 angT UVd angT2 angT3 angT4 angUV1 angUV2 angUV3 angUV4 angUV
               
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);

nn1 = n14.*ones(8640,7560) - imag(angUV1t);
nn2 = n14.*ones(8640,7560) - imag(angUV2t);
nn3 = n14.*ones(8640,7560) - imag(angUV3t);
nn4 = n14.*ones(8640,7560) - imag(angUV4t);
angd1 = nan*ones(8640,7560);
angd2 = angd1; angd3 = angd1; angd4 = angd1;
% 
angUV1t = real((angUV1t)./nn1);
angUV2t = real((angUV2t)./nn2);
angUV3t = real((angUV3t)./nn3);
angUV4t = real((angUV4t)./nn4);

angT1t = real((angT1t)./nn1);
angT2t = real((angT2t)./nn2);
angT3t = real((angT3t)./nn3);
angT4t = real((angT4t)./nn4);

angUV1t(angUV1t<0) = angUV1t(angUV1t<0)  + 360;
angUV2t(angUV2t<0) = angUV2t(angUV2t<0)  + 360;
angUV3t(angUV3t<0) = angUV3t(angUV3t<0)  + 360;
angUV4t(angUV4t<0) = angUV4t(angUV4t<0)  + 360;

angT1t(angT1t<0) = angT1t(angT1t<0) + 360;
angT2t(angT2t<0) = angT2t(angT2t<0) + 360;
angT3t(angT3t<0) = angT3t(angT3t<0) + 360;
angT4t(angT4t<0) = angT4t(angT4t<0) + 360;

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

angd1(angd1>180) = angd1(angd1>180) -360;
angd2(angd2>180) = angd2(angd2>180) -360;
angd3(angd3>180) = angd3(angd3>180) -360;
angd4(angd4>180) = angd4(angd4>180) -360;


angd1(angd1<-180) = angd1(angd1<-180) +360;
angd2(angd2<-180) = angd2(angd2<-180) +360;
angd3(angd3<-180) = angd3(angd3<-180) +360;
angd4(angd4<-180) = angd4(angd4<-180) +360;


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;

  
%% Lat data

Utek_l = Utek_l./n;
Vtek_l = Vtek_l./n;
Txt_l = Txt_l./n;
Tyt_l = Tyt_l./n;
T_l = sqrt(Txt_l.^2 + Tyt_l.^2);
UV_l = sqrt(Utek_l.^2 + Vtek_l.^2);
ratio_l = UV_l./T_l;

angUV_l = atan2(Vtek_l,Utek_l)*180/pi;
anngT_l = atan2(Tyt_l,Txt_l)*180/pi;
angUV_l(angUV_l<0) = angUV_l(angUV_l<0) + 360;
angT_l(angT_l<0) = angT_l(angT_l<0) + 360;

angd_l = angUV_l - angT_l;
angd_l(angd_l>180) = angd_l(angd_l>180) -360;


for lts = 1 : 180
    
    angLats(ii) = mean(mean(angd_l(abs(lats-ii-90)<2&~isnan(angd_l)&angd_l~=0&abs(ratio_l)~=inf)));
    ratioLats(ii) = mean(mean(ratio_l(abs(lats-ii-90)<2&~isnan(ratio_l)&ratio_l~=0&abs(ratio_l)~=inf)));
end
   
    
figure('color','white','units','normalized','outerposition',[0 0 1 1],'name','1');
plot(-90:90,angLats,'r.')
figure('color','white','units','normalized','outerposition',[0 0 1 1],'name','1');
plot(-90:90,ratioLats,'r.')
%% 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)]);

print('-dpng','-r500','ekmanAngle_all2160.png')