% 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=4320;

folderTx = ['~dmenemen/llc_4320/regions/global/oceTAUX_daily'];
folderTy = ['~dmenemen/llc_4320/regions/global/oceTAUY_daily'];


currFldrVek = ['~dmenemen/llc_4320/regions/global/Vek_daily'];
eval(['mkdir ' currFldrVek]);
currFldrUek = ['~dmenemen/llc_4320/regions/global/Uek_daily'];
eval(['mkdir ' currFldrUek]);

siz = [4320 56160 1];

fnamesTx = dir([folderTx '/*.2*']);
fnamesTy = dir([folderTy '/*.2*']);
lats = quikreadpcolor_llc('~dmenemen/llc_4320/grid/YC.data',nx);
lons = quikreadpcolor_llc('~dmenemen/llc_4320/grid/XC.data',nx);
[dx,dy] = quikreadpcolor_dxdy_llc('~dmenemen/llc_4320/grid/DXC.data','~dmenemen/llc_4320/grid/DYC.data',nx);
[dxG,dyG] = quikreadpcolor_dxdy_llc('~dmenemen/llc_4320/grid/DXG.data','~dmenemen/llc_4320/grid/DYG.data',nx);

numFiles = length(fnamesTx);
for ii = 2:5%numFiles
    fnamTx = [folderTx '/' fnamesTx(ii).name];
    fnamTy = [folderTy '/' fnamesTy(ii).name];
    [Tx, Ty]= quikreadpcolor_uv_llc(fnamTx,fnamTy,nx);
    
   Az = .013;
    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_' fnamesTx(ii).name];
    writebin(fileName,Uek);
    
    fileName = [currFldrVek '/Vek_' fnamesTy(ii).name];
    writebin(fileName,Vek);
end