clear
clc
% folderVg = ['~dmenemen/llc_2160/regions/global_4320/Vg_daily/'];
% folderUg = ['~dmenemen/llc_2160/regions/global_4320/Ug_daily/'];
% folderTx = ['~dmenemen/llc_2160/regions/global_4320/oceTAUX_daily/'];
% folderTy = ['~dmenemen/llc_2160/regions/global_4320/oceTAUY_daily/'];

nx = 4320;
siz = [17280 15120 1];
% fnamesVg = dir([folderVg '/*E*']);
% fnamesUg = dir([folderUg '/*E*']);
% fnamesTx = dir([folderTx '/*T*']);
% fnamesTy = dir([folderTy '/*T*']);
% numFiles = length(fnamesVg);
%dx = quikread_llc('~dmenemen/llc_4320/grid/DXC.data',4320,1,'real*4','~dmenemen/llc_4320/grid/',-90,90,0,360);
%dy = quikread_llc('~dmenemen/llc_4320/grid/DYC.data',4320,1,'real*4','~dmenemen/llc_4320/grid/',-90,90,0,360);
% aU = quikreadpcolor_llc('~dmenemen/llc_4320/grid/RAW.data',nx);
% aV = quikreadpcolor_llc('~dmenemen/llc_4320/grid/RAS.data',nx);
lats = quikreadpcolor_llc('~dmenemen/llc_4320/grid/YC.data',nx);
landV = quikreadpcolor_llc('~dmenemen/llc_4320/grid/hFacS.data',nx);
landU = quikreadpcolor_llc('~dmenemen/llc_4320/grid/hFacW.data',nx);

[aV, aU] = quikreadpcolor_RASRAW_llc('~dmenemen/llc_4320/grid/RAS.data','~dmenemen/llc_4320/grid/RAW.data',nx);

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

pin='~dmenemen/llc_4320/MITgcm/run/'; % for ice read in - assume an hourly reading is good enough

currFldrVg = ['/nobackupp8/awinetee/global_4320/Vg_power_30hLFP_hamming_daily'];
eval(['mkdir ' currFldrVg]);
currFldrUg = ['/nobackupp8/awinetee/global_4320/Ug_power_30hLPF_hamming_daily'];
eval(['mkdir ' currFldrUg]);

l = 1-landU;
l = l + circshift(1-landU,3);
l = l+ circshift(1-landU,-3);
l = l+ circshift(1-landU,[0,-3]);
l = l+ circshift(1-landU,[0,3]);
l(l>0) = 1;
landUPad = 1-l;

l = 1-landV;
l = l + circshift(1-landV,3);
l = l+ circshift(1-landV,-3);
l = l+ circshift(1-landV,[0,-3]);
l = l+ circshift(1-landV,[0,3]);
l(l>0) = 1;
landVPad = 1-l;

clear landU landV l 
for ii=10368:(144*24):485568 
        
            
    dy=ts2dte(ii,25,2011,9,10,29)      
        
    fnamTy = ['/nobackupp8/awinetee/global_4320/oceTAUY_30hLPF_hamming_daily/oceTAUY' '_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
    fnamTx = ['/nobackupp8/awinetee/global_4320/oceTAUX_30hLPF_hamming_daily/oceTAUX' '_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
    
    fnamUg = ['/nobackupp8/awinetee/global_4320/Ug_noIce_hamming_daily' '/Ug_noIce_hamming'  int2str((17280)) 'x' int2str((15120)) '.' dy];
    fnamVg = ['/nobackupp8/awinetee/global_4320/Vg_noIce_hamming_daily' '/Vg_noIce_hamming' int2str((17280)) 'x' int2str((15120)) '.' dy];
    
    if exist(fnamUg) && exist(fnamVg) && exist(fnamTy) && exist(fnamTx)
        strNames{n} = dy;
        SIheff = quikreadpcolor_llc([pin 'SIheff.' myint2str(ii,10) '.data'],nx);
        
        Ug = readbin(fnamUg,siz);
        Vg = readbin(fnamVg,siz);
        [Tx, Ty] = quikreadpcolor_uv_llc(fnamTx,fnamTy,nx);
        
        %Remove land(ish)
%         aUT = aU;
%         aVT = aV;
%         aUT(Ug==0|isnan(Ug)) = nan;
%         aVT(Vg==0|isnan(Vg)) = nan;
       
        % remove +-2deg lat 
        Ug(abs(lats)<3) = nan;
        Vg(abs(lats)<3) = nan;
        
        
        
        % remove land
        Ug(landUPad==0) = nan;
        Vg(landVPad==0) = nan;
        
        %su = stdfilt(Ug);
        %sv = stdfilt(Vg);
        
        %Ug(su>.5) = 0;
        %Vg(sv>.5) = 0;
        
        %Ug(abs(lats)>75&su>.1) = 0;
        %Vg(abs(lats)>75&sv>.1) = 0;
        
        %Calculate power per area
        Up = Tx.*Ug;
        Vp = Ty.*Vg;
        
        Up(SIheff>.1) = nan;
        Vp(SIheff>.1) = nan;
        %save
        
        fileName = [currFldrUg '/Ug_power_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
        writebin(fileName,Up);
        
        fileName = [currFldrVg '/Vg_power_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
        writebin(fileName,Vp);
        
        
        
        %Find mean power per area to replace equator region
        %mppa = mean(mean(Up(~isnan(Up)&Up~=0&abs(Up)~=inf))) + mean(mean(Vp(~isnan(Vp)&Vp~=0&abs(Vp)~=inf)));
        mppaU = mean(mean(Up(~isnan(Up)&abs(Up)<1)));
        mppaV = mean(mean(Vp(~isnan(Vp)&abs(Vp)<1)));
        
        upwr(n) = sum(sum(mppaU.*aU(landUPad~=0&~isnan(aU)&SIheff<.1)))
        vpwr(n) =sum(sum(mppaV.*aV(landVPad~=0&~isnan(aV)&SIheff<.1)))
        
        pwr(n) = upwr(n) + vpwr(n)
        
        %sum(sum(mppaU.*aU(landUPad~=0&~isnan(aU)&SIheff<.001))) + sum(sum(mppaV.*aV(landVPad~=0&~isnan(aV)&SIheff<.001)))
        
        % Find area in equator band missiing
        %         missArea = sum(sum(aUT(abs(lats)<2)));
        %         % Calc total power
        %         Up(isnan(Up)|abs(Up)==inf) = 0;
        %         Vp(isnan(Vp)|abs(Vp)==inf) = 0;
        %
        %         aUT(isnan(aUT)) = 0;
        %         aVT(isnan(aVT)) = 0;
        %         pwr(n) = mppa*missArea + sum(sum(Up.*aUT)) + sum(sum(Vp.*aVT))
        
        n = n +1;

      
    end
end
save('geoPowertimeSeries_4320.mat','upwr','vpwr','pwr');
%