% Compare Ug, Umit, Umit-Ug

clear
clc
folderVg = ['~dmenemen/llc_2160/regions/global/Vg_daily/'];
folderUg = ['~dmenemen/llc_2160/regions/global/Ug_daily/'];

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

n = 1;
strNames = {};

currFldrVg = ['~dmenemen/llc_2160/regions/global/Vg_daily'];
eval(['mkdir ' currFldrVg]);
currFldrUg = ['~dmenemen/llc_2160/regions/global/Ug_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);
    
    fnamVm = ['~dmenemen/llc_2160/regions/global/V_daily/V' '_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
    fnamUm = ['~dmenemen/llc_2160/regions/global/U_daily/U' '_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];

    fnamUg = [folderUg 'Ug_Eta_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
    fnamVg = [folderVg 'Vg_Eta_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
     
    if exist(fnamUg) && exist(fnamVg) && exist(fnamUm) && exist(fnamVm)
        n = n + 1
        strNames{n} = dy;
        
        Ug = readbin(fnamUg,siz);
        Vg = readbin(fnamVg,siz);
        Vg= Vg.*landV;
        Ug = Ug.*landU;     % land masks
        
        [Um, Vm] = quikreadpcolor_uv_llc(fnamUm,fnamVm,nx);
       
        
        
        break;
        
        
    end
end
su = stdfilt(Ug);
sv = stdfilt(Vg);
Ugf = Ug;
Vgf = Vg;
Ugf(su>.5) = 0;
Vgf(sv>.5) = 0;

Ugt = (Ugf(1:end-1,:)+Ugf(2:end,:))/2;
Ugt = [Ugt;(Ugf(end,:)+Ugf(1,:))/2];

Vgt = (Vgf(:,1:end-1)+Vgf(:,2:end))/2;
s = size(Vg);
Vgt = [Vgt,nan*ones(s(1),1)];

Mg = sqrt(Vgt.^2 + Ugt.^2);


Umt = (Um(1:end-1,:)+Um(2:end,:))/2;
Umt = [Umt;(Um(end,:)+Um(1,:))/2];

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

Mm = sqrt(Vmt.^2 + Umt.^2);

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

axes(hh(1))
quikpcolor(Mm');
%c = colorbar('southoutside');
caxis([0 .5])
title('Magnitude Total (MITgcm) Velocity 1 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');
c = colorbar('southoutside');
caxis([0 .5])
title('Magnitude Geostrophic Velocity 1 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(3))
quikpcolor(abs(Mm-Mg)');
%c = colorbar('southoutside');
caxis([0 .5])
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('-djpeg','-r3600','MAG_compare_2160.jpeg')