% compare geostrophic kinetic energy flux to total KE flux

% Alex Wineteer
% wineteer@jpl.nasa.gov

clear
clc
folderVgp = ['~dmenemen/llc_2160/regions/global/Vg_power_daily2/'];
folderUgp = ['~dmenemen/llc_2160/regions/global/Ug_power_daily2/'];

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

lats = quikreadpcolor_llc('~dmenemen/llc_2160/grid/YC.data',nx);
[aV aU] = quikreadpcolor_RASRAW_llc('~dmenemen/llc_2160/grid/RAS.data','~dmenemen/llc_2160/grid/RAW.data',nx);
landV = quikreadpcolor_llc('~dmenemen/llc_2160/grid/hFacS.data',nx);
landU = quikreadpcolor_llc('~dmenemen/llc_2160/grid/hFacW.data',nx);

p_tot = 0;
n = 1;
strNames = {};

currFldrVg = ['~dmenemen/llc_2160/regions/global/Vg_power_daily'];
eval(['mkdir ' currFldrVg]);
currFldrUg = ['~dmenemen/llc_2160/regions/global/Ug_power_daily'];
eval(['mkdir ' currFldrUg]);
n = 0;
Vdiff = 0;
Udiff = 0;
Vgpt = 0;
Ugpt  = 0;
Upmt = 0;
Vpmt = 0;

for ii = 92160:(80*24):1586400;
                
    dy=ts2dte(ii,45,2011,1,17,30);
    
    fnamVp = ['~dmenemen/llc_2160/regions/global/V_oceTAUY_daily/V_oceTAUY' '_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
    fnamUp = ['~dmenemen/llc_2160/regions/global/U_oceTAUX_daily/U_oceTAUX' '_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];

    fnamUgp = [folderUgp 'Ug_power_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
    fnamVgp = [folderVgp 'Vg_power_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
     
    if exist(fnamUgp) && exist(fnamVgp) && exist(fnamUp) && exist(fnamVp)
        n = n + 1
        strNames{n} = dy;
        
        Ugp = readbin(fnamUgp,siz);
        Vgp = readbin(fnamVgp,siz);
        
        [Upm, Vpm] = quikreadpcolor_upvp_llc(fnamUp,fnamVp,nx);
        
        Ugp(Ugp == 0 | abs(Ugp) == inf) = nan;
        Vgp(Vgp == 0 | abs(Vgp) == inf) = nan;
        Upm(Upm == 0 | abs(Upm) == inf) = nan;
        Vpm(Vpm == 0 | abs(Vpm) == inf) = nan;
       
        
        Vpmt = Vpmt + Vpm;
        Upmt = Upmt + Upm;
        Ugpt = Ugpt + Ugp;
        Vgpt = Vgpt + Vgp;
        
        Vdiff = Vdiff + (Vpm-Vgp);
        Udiff = Udiff + (Upm-Ugp);
        
       
    end
end
%%
Vgpt = Vgpt.*landV;
Ugpt = Ugpt.*landU;

Ugp = Ugp.*landU;
Vgp = Vgp.*landV;

v = Vgpt;
u = Ugpt;

Vgpt(abs(lats)>65) = nan;
Ugpt(abs(lats)>65) = nan;

%% interpolate to tracers for mag


Ugp_tracer = (Ugpt(1:end-1,:)+Ugpt(2:end,:))/2;
Ugp_tracer = [Ugp_tracer;(Ugpt(end,:)+Ugpt(1,:))/2];

Vgp_tracer = (Vgpt(:,1:end-1)+Vgpt(:,2:end))/2;
s = size(Vgpt);
Vgp_tracer = [Vgp_tracer,nan*ones(s(1),1)];

Mg = Vgp_tracer + Ugp_tracer; %sqrt(Vgp_tracer.^2 + Ugp_tracer.^2);


Ump_tracer = (Upmt(1:end-1,:)+Upmt(2:end,:))/2;
Ump_tracer = [Ump_tracer;(Upmt(end,:)+Upmt(1,:))/2];

Vmp_tracer = (Vpmt(:,1:end-1)+Vpmt(:,2:end))/2;
s = size(Vpmt);
Vmp_tracer = [Vmp_tracer,nan*ones(s(1),1)];

Mm = Vmp_tracer + Ump_tracer; %sqrt(Vmp_tracer.^2 + Ump_tracer.^2);

mnFg = mean(mean(Mg(~isnan(Mg)&abs(Mg)~=inf&Mg~=0&abs(lats)<60&abs(lats)>3)./n));
mnFm = mean(mean(Mm(~isnan(Mm)&abs(Mm)~=inf&Mm~=0&abs(lats)<60&abs(lats)>3)./n));

figure
dp = abs(((((Upmt+Vpmt)./n)-((Ugpt+Vgpt)./n)))./((Upmt+Vpmt)./n));
quikpcolor(dp');
caxis([0 1.5])
figure('color','white','units','normalized','outerposition',[0 0 1 1],'name','1');
hh = tight_subplot(1,3);

axes(hh(1))
quikpcolor((Mm./n)');
%c = colorbar('southoutside');
caxis([-.02 .05])
title('Total Model (MITgcm) KE Flux 90 Day Avg')
%xlabel(c,'Velocity (m/s)')
set(get(gca,'title'),'fontsize',16);
%a = get(c,'xlabel');
%set(a,'fontsize',16);
set(gca,'xtick',[])
set(gca,'ytick',[])
axis([0 8640 1500 6500])
daspect([1 1 1])

axes(hh(2))
quikpcolor((Mg./n)');
c = colorbar('southoutside');
caxis([-.02 .05])
title('Total Geostrophic KE Flux 90 Day Avg')
xlabel(c,'Kinetic Energy Flux (W/m^2)')
set(get(gca,'title'),'fontsize',16);
a = get(c,'xlabel');
set(a,'fontsize',16);
set(gca,'xtick',[])
set(gca,'ytick',[])
axis([0 8640 1500 6500])
daspect([1 1 1])

axes(hh(3))
quikpcolor(abs((Mm-Mg)./n)');
%c = colorbar('southoutside');
caxis([-.02 .05])
title('Absolute Difference: Total - Geostrophic')
%xlabel(c,'Velocity (m/s)')
set(get(gca,'title'),'fontsize',16);
%a = get(c,'xlabel');
%set(a,'fontsize',16);
set(gca,'xtick',[])
set(gca,'ytick',[])
axis([0 8640 1500 6500])
daspect([1 1 1])
        
%print('-dpdf','-r1200','PWR_compare_2160_hi.pdf')