% this script is addapted from compare_RS.m to pull out variables needed to
% make the comparison plots and saves it to a .mat file.
if(~matlabpool('size')>0);
    mpclose = 1;
    matlabpool open
end

%addpath /home/abever/matlab/AARONS_TOOLS/NETCDF_TOOLS/
%addpath /home/abever/matlab/AARONS_TOOLS/PLOTTING/

zl = 4;   % z level to get
lp=0;     % low-pass filter or not

cv2ll = 1;% convert progressive vectors to lat/lon from utm/xy
Sc = 6;    % scale factor to shorten the progressive lat/lon

%diry = '/data/users-tmp/abever/MODELS/WAIPAOA2/OUTPUT/GRID300_4C/';%
diry = '/home/beach/guest/bevera/POVERTY_BAY/OUTPUT/GRID300ub4C/';
offset = 0;

if(~exist('mI','var'))
    disp('Calculating I/J locations')
    StI = [177.95131 177.95373 177.96478 177.96231];
    StJ = [-38.730443 -38.720727 -38.696380 -38.721163];
    
    fname = '/home/beach/guest/bevera/POVERTY_BAY/INPUT/modern_test_grid.nc';
    llm_all
    gname=fname;
    
    parfor i=1:length(StI)
        [mI(i),mJ(i)] = near2d(rlon,rlat,StI(i),StJ(i));
        mI(i)=mI(i)+offset; mJ(i)=mJ(i)+offset;
        A(i) = nc_varget(gname,'angle',[mI(i)-1 mJ(i)-1],[1 1]);
    end
end

if(length(mI) < 5);
    %add a point near the eastern open boundary
    mI(end+1) = 100;  mJ(end+1) = 113;
    for i=1:length(mI)
        A(i) = nc_varget(gname,'angle',[mI(i)-1 mJ(i)-1],[1 1]);
    end
end


if(~exist('RSU','var'))
    run1 = 'GRID300';
    disp(run1)
    fname = '*_his*.nc';
    t = get_variable('ocean_time',fname,diry);
    t = t ./ (3600*24) + 2440000;
    RSU = ones(length(mI),length(t),20) * nan;
    RSV = ones(length(mI),length(t),20) * nan;
    RST = ones(length(mI),length(t),20) * nan;
    RSS = ones(length(mI),length(t),20) * nan;
    RSH = ones(length(mI),length(t)) * nan;
    RSZ = ones(length(mI),length(t)) * nan;
    % not adding sediment stuff yet
    parfor i=1:length(mI);
        RSU(i,:,:) = get_ptvariable('u',fname,[mI(i) mJ(i)],4,diry);
        RSV(i,:,:) = get_ptvariable('v',fname,[mI(i) mJ(i)],4,diry);
        RST(i,:,:) = get_ptvariable('temp',fname,[mI(i) mJ(i)],4,diry);
        RSS(i,:,:) = get_ptvariable('salt',fname,[mI(i) mJ(i)],4,diry);
        %RSH(i,:) = get_ptvariable('Hwave',fname,[mI(i) mJ(i)],3);
        RSZ(i,:) = get_ptvariable('zeta',fname,[mI(i) mJ(i)],3,diry);
    end 
    disp('still going')
    Atmp = ones(length(A),length(t),length(RSS(1,1,:)));
    for i=1:length(t)
        for j=1:length(RSS(1,1,:))
            Atmp(:,i,j) = A;
        end
    end
    W = RSU + sqrt(-1)*RSV;
    W = W .* exp(sqrt(-1).*Atmp);
    RSU = real(W); RSV = imag(W);
end

if(~exist('RSUP','var'));
    disp('Generating progressive vectors')
%    fname = 'G:\WAIPAOA\FORCING\INP_3AUG\modern_test_grid.nc';
    fname = '/home/beach/guest/bevera/POVERTY_BAY/INPUT/modern_test_grid.nc';
    X = nc_varget(fname,'x_rho');
    Y = nc_varget(fname,'y_rho');
    h = nc_varget(fname,'h');
    h(~mask)=nan;

    RSUP = cumsum((RSU*3600*3),2);   RSVP = cumsum((RSV*3600*3),2);

    if(~cv2ll)
        xtmp = X(mI(i),mJ(j(i)));
        ytmp = Y(mI(i),mJ(j(i)));
        parfor i=1:length(mI);
            xtmp = X(mI(i),mJ(j(i)));
            ytmp = Y(mI(i),mJ(j(i)));
            RSUP(i,:,:) = RSUP(i,:,:) + xtmp; %RSUP(i,:,:) = RSUP(i,:,:) - max(max(max(RSUP(i,:,:))));
            RSVP(i,:,:) = RSVP(i,:,:) + ytmp; %RSVP(i,:,:) = RSVP(i,:,:) - max(max(max(RSVP(i,:,:))));
        end
    else
        N = [5712706.5 5713782.5 5716474 5713726 5694767];
        E = [582688.0 582910.1 583898.94 583655.1 607429];
        m_proj('UTM','ellipsoid','wgs84','zone',60,'hem',1)
        RSUPL = ones(size(RSUP))*nan;    RSVPL = ones(size(RSVP))*nan;
        
        for i=1:length(N);
            [RSUPL(i,:,:),RSVPL(i,:,:)] = m_xy2ll(RSUP(i,:,:)./Sc+E(i),RSVP(i,:,:)./Sc+N(i));
        end
    end
end

save comparison_variables.mat RSU RSV RST RSS RSZ mI mJ run1 zl offset ...
    RSUPL RSVPL Sc t

if(exist('mpclose','var'))
    matlabpool close
end