%% Power calculations for MITgcm output. 

% LOAD DATA
clear
clc
folderVp = ['~dmenemen/llc_2160/regions/global/V_oceTAUY_daily/'];
folderUp = ['~dmenemen/llc_2160/regions/global/U_oceTAUX_daily/'];
nx = 2160;
fnamesVp = dir([folderVp '/*.2*']);
fnamesUp = dir([folderUp '/*.2*']);

numFiles = length(fnamesVp);
%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);
%aC = quikreadpcolor_llc('~dmenemen/llc_4320/grid/RAC.data',nx);
[aV aU] = quikreadpcolor_RASRAW_llc('~dmenemen/llc_2160/grid/RAS.data','~dmenemen/llc_2160/grid/RAW.data',nx);

p_tot = 0;
Up_tot = 0;
Vp_tot = 0;
n = 1;
strNames = {};
%for ii = 92160+(80*24):1586400;
for  ii=92160:(80*24):92160+1*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];
    
    %fnamVp = ['~dmenemen/llc_4320/regions/global/V_oceTAUY_daily/' fnamesVp(ii).name]
    %fnamUp = ['~dmenemen/llc_4320/regions/global/U_oceTAUX_daily/' fnamesUp(ii).name]
    
    if exist(fnamVp) && exist(fnamUp)
        strNames{n} = dy;
        [Up Vp] = quikreadpcolor_upvp_llc(fnamUp,fnamVp,nx);
        % Up = quikreadpcolor_llc(fnamUp,nx);
        % power per area
        
        Vp(Vp == 0) = nan;
        Up(Up == 0) = nan;
        
        ppV = Vp.*aV ;
        ppU =  Up.*aU;
        
        Up_tot(n) = sum(sum(ppU(~isnan(ppU))));
        Vp_tot(n) = sum(sum(ppV(~isnan(ppV))));
        
        p_tot(n) = Vp_tot(n) + Up_tot(n);
        n = n + 1
    end
    
            
end
% startDate = datenum('03-07-2011');
% endDate = datenum('11-06-2011');
% xData = linspace(startDate,endDate,9);
% 
% 
% figure('color','white')
% plot(p_tot,'r.-','markeredgecolor','r')
% hold on
% plot(Up_tot,'k.-','markeredgecolor','k')
% plot(Vp_tot,'b.-','markeredgecolor','b')
% 
% axis([0 334 10^11 5.5*10^12])
% text(50,2.25*10^12,['Mean Power: ' num2str(mean(p_tot)/10^12) ' TW'],'color','r')
% text(50,2*10^12,['Median Power: ' num2str(median(p_tot)/10^12) ' TW'],'color','r')
% text(150,2.25*10^12,['Mean Power: ' num2str(mean(Vp_tot)/10^12) ' TW'],'color','b')
% text(150,2*10^12,['Median Power: ' num2str(median(Vp_tot)/10^12) ' TW'],'color','b')
% text(250,2.25*10^12,['Mean Power: ' num2str(mean(Up_tot)/10^12) ' TW'],'color','k')
% text(250,2*10^12,['Median Power: ' num2str(median(Up_tot)/10^12) ' TW'],'color','k')
% 
% 
% %xlabel('~Days Since 2160 Sim Start')
% ylabel('Power (Watts)')
% title('Global Power Input for 1/24 degree model')
% ax = gca;
% set(ax,'XTick',[0:334/11:334])
% set(ax,'XTickLabel', ['Mar';'Apr';'May';'Jun';'Jul';'Aug';'Sep';'Oct';'Nov';'Dec';'Jan';'Feb']);
% set(ax,'FontSize',14);
% set(get(ax,'ylabel'),'fontsize',16);
% set(get(ax,'title'),'fontsize',16);
% legend('Total Power','U Power','V Power')
% print('-dpng','-r300','time_series_globalPwr_mar_feb_2011_2160.png')

fig = figure('color','white','units','normalized','outerposition',[0 0 1 1],'name','1');
quikpcolor((Vp + Up)');
c = colorbar;
caxis([-.02 .05])
title('Total Global 1/24 degree Wind Power Input for 8 day hourly snapshot avg (March 7 start)')
xlabel(c,'Power Per Area (W/m^2)')
set(get(gca,'title'),'fontsize',16);
a = get(c,'xlabel');
set(a,'fontsize',16);
set(gca,'xtick',[])
set(gca,'ytick',[])
daspect([1 1 1])
print('-dpng','-r300','TOTAL_globalPwr_8d_hourly_snapshot_2160.png')

fig = figure('color','white','units','normalized','outerposition',[0 0 1 1],'name','1');
quikpcolor(Vp');
c = colorbar;
caxis([-.02 .05])
title('V Portion of Global 1/24 degree Wind Power Input for March 7, 2011 (1d avg)')
xlabel(c,'Power Per Area (W/m^2)')
set(get(gca,'title'),'fontsize',16);
a = get(c,'xlabel');
set(a,'fontsize',16);
set(gca,'xtick',[])
set(gca,'ytick',[])
daspect([1 1 1])
%print('-dpng','-r300','V_globalPwr_mar_07_2011_2160.png')


fig = figure('color','white','units','normalized','outerposition',[0 0 1 1],'name','1');
quikpcolor(Up');
c = colorbar;
caxis([-.02 .05])
title('U Portion of Global 1/24 degree Wind Power Input for March 7, 2011 (1d avg)')
xlabel(c,'Power Per Area (W/m^2)')
set(get(gca,'title'),'fontsize',16);
a = get(c,'xlabel');
set(a,'fontsize',16);
set(gca,'xtick',[])
set(gca,'ytick',[])
daspect([1 1 1])
%print('-dpng','-r300','U_globalPwr_mar_07_2011_2160.png')