% ekman spiral
clear

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);
dyG = readbin(['~dmenemen/llc_2160/regions/' region_name 'grid/DXG_' sizstr],siz1);
dxG = readbin(['~dmenemen/llc_2160/regions/' region_name 'grid/DYG_' sizstr],siz1);
z = readbin('~dmenemen/llc_2160/grid/RC.data',[27]);

dxG = circshift(dxG,[0 -1]);
dxG(:,end) = nan;
fc = 2*7.2921*10^-5*sind(lats);

g = 9.806;

pin = ['/nobackupp9/dmenemen/llc_2160/regions/' region_name '/averages_winter/'];

eta = readbin([pin 'Eta_winterAVG'],siz1);

if strcmp(region_name,'DopplerScat/southern/')
    Umodel = readbin([pin 'U_winterAVG'],siz27);
    Vmodel = readbin([pin 'V_winterAVG'],siz27);
    dyG = readbin(['~dmenemen/llc_2160/regions/' region_name 'grid/DYG_' sizstr],siz1);
    dxG = readbin(['~dmenemen/llc_2160/regions/' region_name 'grid/DXG_' sizstr],siz1);
else
    Umodel = readbin([pin 'V_winterAVG'],siz27);
    Vmodel = -circshift(readbin([pin 'U_winterAVG'],siz27),[0 -1]);
    Vmodel(:,end) = nan;
end

Ty = readbin([pin 'oceTAUY_winterAVG'],siz1);
Tx= readbin([pin 'oceTAUX_winterAVG'],siz1);
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);

%% plot


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

axes(hh(1))
quikpcolor(Vmit')
colorbar
caxis([-.2 .15])
title('Vg')
axes(hh(2))
quikpcolor(Umit');
colorbar
caxis([-.2 .15])
title('Ug')
axes(hh(3))
quikpcolor(Vmodel(:,:,22)');
colorbar
caxis([-.2 .15])
title('Vmodel @ 100m')
axes(hh(4))
quikpcolor(Umodel(:,:,22)');
colorbar
caxis([-.2 .15])
title('Umodel @ 100m')
%% Compute Ekman Spiral

for ii = 1 : 27
    Udiff(:,:,ii) = Umodel(:,:,ii) - Umit;
    Vdiff(:,:,ii) = Vmodel(:,:,ii) - Vmit;
    uu = Udiff(:,:,ii);
    vv = Vdiff(:,:,ii);
    
    
    Uek(ii) = mean(mean(uu(~isnan(uu)&uu~=0)));
    Vek(ii) = mean(mean(vv(~isnan(vv)&vv~=0)));
   
end
Txm = mean(mean(Tx(~isnan(Tx)&Tx~=0)));
Tym = mean(mean(Ty(~isnan(Ty)&Ty~=0)));

angTau = (atan2(Tym,Txm))*180/pi;
angEk = (atan2(Vek,Uek))*180/pi;
UVek = sqrt(Uek.^2 + Vek.^2);

rotAngEk  =  angEk + angTau;
UVek_n = UVek./max(UVek);

figure
plot(UVek_n);

figure
quiver3(zeros(22,22),zeros(22,22),repmat(z(1:22),1,22),diag(Uek(1:22))./max(Uek(1:22)),diag(Vek(1:22))./max(Vek(1:22)),zeros(22,22)')

save(['/nobackupp9/dmenemen/llc_2160/regions/' region_name '/averages_winter/southern_winter_ang_mag.mat'],'angTau','angEk','rotAngEk','UVek','UVek_n','Uek','Vek');