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

%{
    This code will calculate the ekman currents based on the input ocean
    surface shear. 

    INPUTS: surface shear folder locations, both Tx and Ty
            size of grid [x,y,z] -- include z, even if it is just 1.

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

    Grid Output Locations: (NOTE: actual grid output is flipped up down and
    transposed relative 
    to this grid. This output matches the MITgcm U/V locations, but is
    missing the first and last rows. NaN has been input in the 
    first and last rows (latitude bands)
    
    |         |         |    
    |____^____|_________|___
    |   V21   |         |
    |         |         |
    |         |         |
    |____^____|_________|___
    |   V11   |         |
    |      U11->     U12->
    |         |         |
    |_________|_________|___
    *

%}
clear
region_name='global';
nx=2160;

folderTx = ['/nobackupp8/awinetee/global/oceTAUX_30hLPF_hamming_daily'];
folderTy = ['/nobackupp8/awinetee/global/oceTAUY_30hLPF_hamming_daily'];

% folderTxo = ['~dmenemen/llc_2160/regions/global/oceTAUX_daily'];
% folderTyo = ['~dmenemen/llc_2160/regions/global/oceTAUY_daily'];


currFldrVek = ['/nobackupp8/awinetee/global/Vek_30hLPF_hamming_daily'];
eval(['mkdir ' currFldrVek]);
currFldrUek = ['/nobackupp8/awinetee/global/Uek_30hLPF_hamming_daily'];
eval(['mkdir ' currFldrUek]);

%siz = [4320 56160 1];

fnamesTx = dir([folderTx '/*.2*']);
fnamesTy = dir([folderTy '/*.2*']);
% fnamesTxo = dir([folderTxo '/*.2*']);
% fnamesTyo = dir([folderTyo '/*.2*']);

lats = quikreadpcolor_llc('~dmenemen/llc_2160/grid/YC.data',nx);
lons = quikreadpcolor_llc('~dmenemen/llc_2160/grid/XC.data',nx);
[dx,dy] = quikreadpcolor_dxdy_llc('~dmenemen/llc_2160/grid/DXC.data','~dmenemen/llc_2160/grid/DYC.data',nx);
[dxG,dyG] = quikreadpcolor_dxdy_llc('~dmenemen/llc_2160/grid/DXG.data','~dmenemen/llc_2160/grid/DYG.data',nx);

numFiles = length(fnamesTy);
%started at numFiles = 526/750
for ii = 400:numFiles
    fnamTx = [folderTx '/' fnamesTx(ii).name];
    fnamTy = [folderTy '/' fnamesTy(ii).name];
%   fnamTxo = [folderTxo '/' fnamesTxo(ii).name];
%     fnamTyo = [folderTyo '/' fnamesTyo(ii).name];

    
    %txname = strsplit([fnamesTx(ii).name],'.');
    % might want to implent a check here to make sure Tx and Ty are the
    % same day. With the data I'm using, there doesn't look likee there's
    % a problem.
    
    [Tx, Ty]= quikreadpcolor_uv_llc(fnamTx,fnamTy,nx);
    %[Txo, Tyo]= quikreadpcolor_uv_llc(fnamTxo,fnamTyo,nx);
    Az = .0104;
    rho = 1027;
    f = 2*7.2921*10^-5*sind(lats);
    
    Txc = [(Tx(1:end-1,:) + Tx(2:end,:))./2; (Tx(end,:)+Tx(1,:))/2];
    s = size(Ty);
    Tyc = [(Ty(:,1:end-1)+Ty(:,2:end))/2,nan*ones(s(1),1)];
    normT = sqrt(Txc.^2 + Tyc.^2);
    
    angleUV = atan2(Tyc,Txc);
    V0 = normT./(sqrt(abs(f).*rho.^2.*Az));
    
    Uc1 = V0.*cos(angleUV - pi/4);
    Vc1 = V0.*sin(angleUV - pi/4);
    Uc2 = V0.*cos(angleUV + pi/4);
    Vc2 = V0.*sin(angleUV + pi/4);
    
    Uc1(lats<=0) = 0;
    Vc1(lats<=0) = 0;
    Uc2(lats>0) = 0;
    Vc2(lats>0) = 0;
    
    Uc = real(Uc1 + Uc2);
    Vc = real(Vc1 + Vc2);
    
    Uek = [(Uc(end,:) + Uc(1,:))/2;(Uc(1:end-1,:) + Uc(2:end,:))/2];
    Vek = [nan*ones(s(1),1),(Vc(:,1:end-1) + Vc(:,2:end))/2;];
     
    fileName = [currFldrUek '/Uek_30hLPF_hamming_' fnamesTx(ii).name];
    writebin(fileName,Uek);
    
    fileName = [currFldrVek '/Vek_30hLPF_hamming_' fnamesTy(ii).name];
    writebin(fileName,Vek);
end