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

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

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


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




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);

siz = size(lats);

for ii=10368+50*144*24:(144*24):485568
    
    
    dy=ts2dte(ii,25,2011,9,10,29)
    
    fnamTy = ['/nobackupp8/awinetee/global_4320/oceTAUY_30hLPF_hamming_daily/oceTAUY' '_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
    fnamTx = ['/nobackupp8/awinetee/global_4320/oceTAUX_30hLPF_hamming_daily/oceTAUX' '_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
    
    if exist(fnamTy) && exist(fnamTx)
        [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_' num2str(siz(1)) 'x' num2str(siz(2)) '_' dy];
        writebin(fileName,Uek);
        
        fileName = [currFldrVek '/Vek_30hLPF_hamming_' num2str(siz(1)) 'x' num2str(siz(2)) '_' dy];
        writebin(fileName,Vek);
        
        clear Vek Uek Uc1 Uc2 Vc1 Vc2 Uc Vc Txc Tyc angleUV V0
    end
end