
% Surface Geostrophic Currents
% Alex Wineteer
% Email: wineteer@jpl.nasa.gov, awinetee@calpoly.edu

%{
    This code will calculate the surface geostrophic currents based on the 
    input ocean surface heights. 

    INPUTS: surface height folder locations (Eta), 
            size of grid [x,y,z] -- include z, even if it is just 1.
            grid folder containing lats, dxC, dyC, dxG, dyG

    OUTPUTS: NONE, though it will write the currents to disk under Ug/Ug
    directories. Also note that these output binary files will be written
    in a pcolor-like format -- not the default MITgcm format. THis means
    that reading them only requires a call to readbin.m
    Grid Output Locations: (NOTE: actual grid output is flipped up-down, and
                             transposed relative to this grid -- at least, in MATLAB)
    
    |         |         |    
    |____^____|_________|___
    |   V21   |         |
    |         |         |
    |         |         |
    |____^____|_________|___
    |   V11   |         |
    |     U11 ->    U12 ->
    |         |         |
    |_________|_________|___
    *

%}


clear
close all
%folder = ['~dmenemen/llc_1080/regions/global/Eta_daily'];
%siz = [4320 56160 1];
  
%fnames = dir([folder '/*.2*']);
%numFiles = length(fnames);

%currFldrVg = ['~dmenemen/llc_1080/regions/global/Vg_daily'];
currFldrVg = ['/nobackupp8/awinetee/global_1080_test/Vg_noIce_hamming_daily'];
eval(['mkdir ' currFldrVg]);
%currFldrUg = ['~dmenemen/llc_1080/regions/global/Ug_daily'];
currFldrUg = ['/nobackupp8/awinetee/global_1080_test/Ug_noIce_hamming_daily'];
eval(['mkdir ' currFldrUg]);
nx = 1080;
lats = quikreadpcolor_llc('~dmenemen/llc_1080/grid/YC.data',nx);
%lons = quikreadpcolor_llc('~dmenemen/llc_1080/grid/XC.data',nx);
landV = quikreadpcolor_llc('~dmenemen/llc_1080/grid/hFacS.data',nx);
landU = quikreadpcolor_llc('~dmenemen/llc_1080/grid/hFacW.data',nx);


[dx,dy] = quikreadpcolor_dxdy_llc('~dmenemen/llc_1080/grid/DXC.data','~dmenemen/llc_1080/grid/DYC.data',nx);
[dxG,dyG] = quikreadpcolor_dxGdyG_llc('~dmenemen/llc_1080/grid/DXG.data','~dmenemen/llc_1080/grid/DYG.data',nx);

fc = 2*7.2921*10^-5*sind(lats);

g = 9.806;
rho = 1027.5;
SEAICE_rhoIce = 910;
SEAICE_rhoSnow = 330;
pin='~dmenemen/llc_1080/MITgcm/run_day732_896/';
%fnames = dir('~dmenemen/llc_1080/regions/global/Eta_daily/E*');
fnames = dir('/nobackupp8/awinetee/global_1080_test/Eta_seaIceRemoved_30hLPF_hamming_daily/E*');

%numFiles = length(fnames);
strNames = {};
n = 0;
for ii=700800:(40*24*5):700800+(40*24*5)*30 
        
      
    dyy=ts2dte(ii,90,2010,1,1)
      
      
    
    fnam = ['/nobackupp8/awinetee/global_1080_test/Eta_seaIceRemoved_30hLPF_hamming_daily/' 'Eta_seaIceRemoved_' int2str((nx)) 'x' int2str((nx*13)) '.' dyy];
    fnamIce = [pin 'SIheff.' myint2str(ii,10) '.data'];
    %fnamSnow = [pin 'SIhsnow.' myint2str(ii,10) '.data'];
    
    if exist(fnam)% && exist(fnamIce) && exist(fnamSnow)
        n = n + 1;
        strNames{n} = dyy;
        
       SIheff = quikreadpcolor_llc([pin 'SIheff.' myint2str(ii,10) '.data'],nx);
       % SIhsnow = quikreadpcolor_llc([pin 'SIhsnow.' myint2str(ii,10) '.data'],nx);
      
        
        
        
        %% FIX THIS PORTION TO DO MEDIAN FILTERING FOR ICE BUT MAAKE IT DO IT FAST

        eta = quikreadpcolor_llc(fnam,nx);
        
%         etaSTD = stdfilt(eta,ones(25));
%         
% 
%         etaMedFilt = medfilt2(eta,[11,11]);
% 
% 
%         eta(etaSTD > .2&(lats>65|lats<-65)&SIheff>.01) = etaMedFilt(etaSTD > .2&(lats>65|lats<-65)&SIheff>.01);
        clear etaSTD etaMedFilt SIheff
      
%         figure
%         quikpcolor(etaMedFilt');
%         colorbar
%         caxis([-2 2]);
        %%
        
       % SIheff(isnan(SIheff)) = 0;
       % SIhsnow(isnan(SIhsnow)) = 0;
       % iceLoad = SIheff*SEAICE_rhoIce + SIhsnow*SEAICE_rhoSnow;
        
        %ssh = eta;% + iceLoad/rho;
        s = size(eta);
        ss= s;
        %eta = ssh;
        deltaEtaX = [diff(eta,1,1);eta(1,:)-eta(end,:)]; %[eta(2:end,:) - eta(1:end-1,:);eta(1,:)-eta(end,:)];  % diff between cell enters -- therefore centered on boundary
        deltaEtaY = [eta(:,2:end) - eta(:,1:end-1),nan*ones(s(1),1)];
        
        clear eta
        
        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;
        
        clear ffY ffX deltaEtaY deltaEtaX
        Us =  Us.*landV;
        Vs = Vs.*landU;     % switched b/c edges are switched right now
        
        % Vs(abs(lats)<15) = nan;
        % Us(abs(lats)<15) = nan;
        
      
        
        s= size(Us);
        Us1 = (Us(:,1:end-1) + Us(:,2:end))/2;
        Us1(:,end+1) = nan*ones(s(1),1);
        
        clear Us;
        
        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,:));
        
        clear Us1;
        
        Vs1 = (Vs(1:end-1,:) + Vs(2:end,:))/2;
        s = size(Vs1);
        Vs1(end+1,:) = (Vs(1,:)+Vs(end,:))/2;
        
        clear Vs;
        
        Vmit = Vs1(:,1:end-1) + (.5.*dyG(:,1:end-1)./dy(:,2:end)).*(Vs1(:,2:end)-Vs1(:,1:end-1));
        
        clear Vs1;
        
        s = size(Vmit);
        Vmit(:,end+1) = zeros(s(1),1);
        
        Vmit(isnan(Vmit)) = 0;
        Umit(isnan(Umit)) = 0;
        
        Vmit = Vmit.*landV;
        Umit = Umit.*landU;     % land masks
        s = ss;
        
%         
%         figure
%         quikpcolor(Umit');
%         caxis([-.2 .5])
%         figure
%         quikpcolor(Vmit');
%         caxis([-.2 .5])
%         pause
%         
        fileName = [currFldrUg '/Ug_noIce_hamming' int2str(s(1)) 'x' int2str(s(2)) '.' dyy];
        writebin(fileName,Umit);
        
        fileName = [currFldrVg '/Vg_noIce_hamming' int2str(s(1)) 'x' int2str(s(2)) '.' dyy];
        writebin(fileName,Vmit);
    end
end

%%
% figure(3)
% quikpcolor(sqrt(Uout.^2 + Vout.^2)');
% colorbar
% caxis([0 1])
% daspect([1 1 1])
% title('Umag 27 Day Avg')
% print('-dpng','-r800','MAGmit_27day_global.png')
% figure(4)
% quikpcolor(Vout');
% colorbar
% caxis([-.3 1])
% title('Vmit 27 Day Avg')
% daspect([1 1 1])
% print('-dpng','-r800','Vmit_27day_global.png')
