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

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

    INPUTS: surface height folder locations, 
            size of grid [x,y,z] -- include z, even if it is just 1.

    OUTPUTS: NONE, though it will write the currents to disk under
    the same folders used for input. (inputFolder/surfGeoCur)

    Grid Output Locations: (NOTE: actual grid output is transposed relative
    to this grid.)
    
    |         |         |    
    |____>____|_________|___
    |   U21   |         |
    |         |         |
    |         |         |
    |____>____|_________|___
    |   U11   |         |
    |      V11^      V12^
    |         |         |
    |_________|_________|___
    *

%}
function surfGeoCur(folder,siz)

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

currFldr = ['~dmenemen/llc_4320/regions/global/surfGeoCur'];
eval(['mkdir ' currFldr]);
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);
lats = quikread_llc('~dmenemen/llc_4320/grid/YC.data',4320,1,'real*4','~dmenemen/llc_4320/grid/',-90,90,0,360);
lons = quikread_llc('~dmenemen/llc_4320/grid/XC.data',4320,1,'real*4','~dmenemen/llc_4320/grid/',-90,90,0,360);

for ii = 1 : length(dy)
    f{ii} = 2*7.2921*10^-5*sind(lats{ii});
end
g = 9.806;

for ii = 1 : 1
    
    clear  deltaEtaX deltaEtaY;
    
    %fnam = [folder '/' fnames(ii).name];
    %fnam ='~dmenemen/llc_4320/regions/global/V_daily/V_4320x56160.2011-10-24';
    fnam ='~dmenemen/llc_4320/regions/global/Eta_daily/Eta_4320x56160.2011-10-24';
    %fld = rot90(readbin(fnam,siz),1);
    [fld,~,~,~] = quikread_llc(fnam,4320,1,'real*4','~dmenemen/llc_4320/grid/',-90,90,0,360);
    for fc = 1 : length(fld)
        if fc <= 4 && fc ~=3
            
            %deltaEtaX = [diff(fld{fc},1,1);nan*ones(1,12960)];  % diff between cell enters -- therefore centered on boundary
            %deltaEtaY = [diff(fld{fc},1,2),nan*ones(4320,1)];
            
            deltaEtaX = [fld{fc}(2:end,:) - fld{fc}(1:end-1,:);nan*ones(1,12960)];  % diff between cell enters -- therefore centered on boundary
            deltaEtaY = [fld{fc}(:,2:end) - fld{fc}(:,1:end-1),fld{fc+1}(:,1) - fld{fc}(:,end)];
             
            ffY = [((f{fc}(:,1:end-1) + f{fc}(:,2:end))/2),(f{fc+1}(:,1)+f{fc}(:,end))/2];
            ffX = [(f{fc}(1:end-1,:) + f{fc}(2:end,:)) / 2;nan*ones(1,12960)];  
      
        elseif fc > 4
            %deltaEtaX = [diff(fld{fc},1,1);nan*ones(1,12960)];
            %deltaEtaY = [diff(fld{fc},1,2),nan*ones(4320,1)];
            
            deltaEtaX = [fld{fc}(2:end,:) - fld{fc}(1:end-1,:);nan*ones(1,12960)];  % diff between cell enters -- therefore centered on boundary
            deltaEtaY = [fld{fc}(:,2:end) - fld{fc}(:,1:end-1),nan*ones(4320,1)];
            ffY = [((f{fc}(:,1:end-1) + f{fc}(:,2:end))/2),nan*ones(4320,1)];
            ffX = [(f{fc}(1:end-1,:) + f{fc}(2:end,:)) / 2;nan*ones(1,12960)];  
      
        else
            continue;
        end
        
        Us = -(g ./ ffY .* deltaEtaY ./ dy{fc});
        Vs = g ./ ffX .* deltaEtaX ./ dx{fc};
        
        Vs(isnan(Vs)) = 0;
        Vg{fc} = Vs;
        Us(isnan(Us)) = 0;
        Ug{fc} = Us;
        
        %fv{fc} = deltaEtaX;
        %Us{fc} = cat(1,nan*ones(1,siz(1)),Us);
        %Vs{fc} = cat(2,nan*ones(siz(2),1),Vs);
    end
%     Us = rot90(Us,3);   % Rotate things back to the orientation they came in on
%     Vs = rot90(Vs,3);
    
    fileName = [currFldr '/Us_' fnames(ii).name];
    writebin(fileName,Us);
     
    fileName = [currFldr '/Vs_' fnames(ii).name];
    writebin(fileName,Vs);
        
end