% 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
folder = ['~dmenemen/llc_4320/regions/global/Eta_daily'];
siz = [17280 15120 1];
  
fnames = dir([folder '/*.2*']);
numFiles = length(fnames);

currFldrVg = ['~dmenemen/llc_4320/regions/global/Vg_daily'];
eval(['mkdir ' currFldrVg]);
currFldrUg = ['~dmenemen/llc_4320/regions/global/Ug_daily'];
eval(['mkdir ' currFldrUg]);
nx = 4320;
lats = quikreadpcolor_llc('~dmenemen/llc_4320/grid/YC.data',nx);
lons = quikreadpcolor_llc('~dmenemen/llc_4320/grid/XC.data',nx);
landV = quikreadpcolor_llc('~dmenemen/llc_4320/grid/hFacS.data',nx);
landU = quikreadpcolor_llc('~dmenemen/llc_4320/grid/hFacW.data',nx);
[dx,dy] = quikreadpcolor_dxdy_llc('~dmenemen/llc_4320/grid/DXC.data','~dmenemen/llc_4320/grid/DYC.data',nx);
[dxG,dyG] = quikreadpcolor_dxGdyG_llc('~dmenemen/llc_4320/grid/DXG.data','~dmenemen/llc_4320/grid/DYG.data',nx);

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

g = 9.806;
fnames = dir('~dmenemen/llc_4320/regions/global/Eta_daily/E*');
numFiles = length(fnames);
for ii = 1 : 5
    ii
    fnam = ['~dmenemen/llc_4320/regions/global/Eta_daily/' fnames(ii).name];%'~dmenemen/llc_4320/regions/global/Eta_daily/Eta_4320x56160.2011-10-24';
    
    eta = quikreadpcolor_llc(fnam,nx);
    s = size(eta);
    
    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)];
    
    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;
    
    % Vs(abs(lats)<15) = nan;
    % Us(abs(lats)<15) = nan;
    
    clear eta ffY ffX deltaEtaX deltaEtaY
    
    s= size(Us);
    Us1 = (Us(:,1:end-1) + Us(:,2:end))/2;
    Us1(:,end+1) = nan*ones(s(1),1);
    
    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,:));
    
    Vs1 = (Vs(1:end-1,:) + Vs(2:end,:))/2;
    s = size(Vs1);
    Vs1(end+1,:) = (Vs(1,:)+Vs(end,:))/2;
    
    Vmit = Vs1(:,1:end-1) + (.5.*dyG(:,1:end-1)./dy(:,2:end)).*(Vs1(:,2:end)-Vs1(:,1:end-1));
    s = size(Vmit);
    Vmit(:,end+1) = zeros(s(1),1);
    
    Vmit = Vmit.*hFacS;
    Umit = Umit.*hFacW;
    
    Vmit(isnan(Vmit)) = 0;
    Umit(isnan(Umit)) = 0;
    
    fileName = [currFldrUg '/Ug_' fnames(ii).name];
    writebin(fileName,Umit);
    
    fileName = [currFldrVg '/Vg_' fnames(ii).name];
    writebin(fileName,Vmit);
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')