%
% This file can be used to generate sediment transport vectors for calculating
% sediment fluxes and for using in figures.
%
%       Aaron Bever, 2009,  DCF 2010

%diry='C:\Documents and Settings\abever\Desktop\Transfer\ROMS_SWAN\';
%diry='/data/users-tmp/abever/MODELS/ROMS/S_BF2_12/'

mdir = '../OUTAUGSed/'
ncdfF = dir([mdir,'calus2008_his*.nc']);
%grid = '/home/beach/guest/bevera/SWF/FORCING/calus_grid2.nc'
grid = '../FRC2008/swf_grid8.nc';
nfiles=length(ncdfF);

fname = [mdir ncdfF(1).name];
llm_all;


if(~exist('iiu','var'))
        % this is because griddata is so much slower than interp2, 
        % method copied from stuff Rich Signell had written
        [nx,ny] = size(ulon);
        [iiu,jju] = meshgrid(1:ny,1:nx);
        % x = [ ulong(:,1); ulong(end,:).'; ulong([end:-1:1],end); ulong(1,[end:-1:1]).'];
        % y = [ ulatg(:,1); ulatg(end,:).'; ulatg([end:-1:1],end); ulatg(1,[end:-1:1]).'];
        % tmp = insider(rlong(:),rlatg(:),x,y);
        % in = find(tmp == 1);
        ibu = griddata(ulon,ulat,iiu,rlon,rlat);
        jbu = griddata(ulon,ulat,jju,rlon,rlat);
        keepiu = find(isnan(ibu) == 1);  ibu(keepiu) = 1;
        keepju = find(isnan(jbu) == 1);  jbu(keepju) = 1;

        [nx,ny] = size(vlon);
        [iiv,jjv] = meshgrid(1:ny,1:nx);
        % x = [ vlong(:,1); vlong(end,:).'; vlong([end:-1:1],end); vlong(1,[end:-1:1]).'];
        % y = [ vlatg(:,1); vlatg(end,:).'; vlatg([end:-1:1],end); vlatg(1,[end:-1:1]).'];
        % tmp = insider(rlong(:),rlatg(:),x,y);
        % in = find(tmp == 1);
        ibv = griddata(vlon,vlat,iiv,rlon,rlat);
        jbv = griddata(vlon,vlat,jjv,rlon,rlat);
        keepiv = find(isnan(ibv) == 1);  ibv(keepiv) = 1;
        keepjv = find(isnan(jbv) == 1);  jbv(keepjv) = 1;
end

Ut = [];
Vt = [];

for j=1:nfiles;
    clear U V Z dz tmp
    fname=[mdir ncdfF(j).name];
    disp(fname);
     
    u=nc_varget(fname,'u');         u(:,:,~umask) = nan;
    v=nc_varget(fname,'v');         v(:,:,~vmask) = nan;
    m1=nc_varget(fname,'mud_01');
    m2=nc_varget(fname,'mud_02');
    m3=nc_varget(fname,'mud_03');
    
      
    M = m1+m2+m3;                      M(:,:,~mask) = nan;

    disp('gridding')
    for k=1:length(u(:,1,1,1))
        for K=1:length(u(1,:,1,1))
            U(k,K,:,:) = interp2(iiu,jju,squeeze(u(k,K,:,:)),ibu,jbu);
            V(k,K,:,:) = interp2(iiv,jjv,squeeze(v(k,K,:,:)),ibv,jbv);
        end
    end

    for i=1:length(u(:,1,1,1))
        Z(i,:,:,:) = depths(fname,grid,5,1,i);
    end
    dz = diff(Z,1,4);
    dz = permute(dz,[1 4 2 3]);
    
    tmp = squeeze(sum(U.*M.*dz,2));
    Ut = [Ut;tmp];
    tmp = squeeze(sum(V.*M.*dz,2));
    Vt = [Vt;tmp];

end

Utm = squeeze(mean(Ut));
Vtm = squeeze(mean(Vt));

Uts = squeeze(sum(Ut));
Vts = squeeze(sum(Vt));

README = {'This is the U and V directional sediment transport.'
    'The units of Ut and Vt are kg/ms. So they are depth integrated'
    'throughout the water column.'}

tmpDIR = '../OUTAUGSed';
tmpDIR = [tmpDIR,'/'];
t = get_variable('ocean_time','calus2008_*.nc',tmpDIR);
t = t./(3600*24)+2440000;
transport_time = t;

save sed_trans_paths_AUGSed.mat Ut Vt README transport_time
    
caca

    
