% Making a script to compare different ROMS model runs.

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

diry1 = '/data/users-tmp/abever/MODELS/WAIPAOA2/OUTPUT/GRID300_4C/';
diry2 = '/data/users-tmp/abever/MODELS/WAIPAOA2/OUTPUT/GRID380/';
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 = '/data/users-tmp/abever/MODELS/WAIPAOA2/FORCING/INP_3AUG/modern_test_grid.nc';
    llm_all
    gname=fname;
    
    for i=1:length(StI)
        [mI(i),mJ(i)] = near2d(rlon,rlat,StI(i),StJ(i));
        mI(i)=mI(i)+3; mJ(i)=mJ(i)+3;
        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,diry1);
    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
    for i=1:length(mI);
        RSU(i,:,:) = get_ptvariable('u',fname,[mI(i) mJ(i)],4,diry1);
        RSV(i,:,:) = get_ptvariable('v',fname,[mI(i) mJ(i)],4,diry1);
        RST(i,:,:) = get_ptvariable('temp',fname,[mI(i) mJ(i)],4,diry1);
        RSS(i,:,:) = get_ptvariable('salt',fname,[mI(i) mJ(i)],4,diry1);
        %RSH(i,:) = get_ptvariable('Hwave',fname,[mI(i) mJ(i)],3);
        RSZ(i,:) = get_ptvariable('zeta',fname,[mI(i) mJ(i)],3,diry1);
    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('RU','var'))
    run2 ='GRID380';
    disp(run2)
    fname = '*_his*.nc';
    t2 = get_variable('ocean_time',fname,diry2);
    t2 = t2 ./ (3600*24) + 2440000;
    RU = ones(length(mI),length(t2),20) * nan;
    RV = ones(length(mI),length(t2),20) * nan;
    RT = ones(length(mI),length(t2),20) * nan;
    RS = ones(length(mI),length(t2),20) * nan;
    RH = ones(length(mI),length(t2)) * nan;
    RZ = ones(length(mI),length(t2)) * nan;
    % not adding sediment stuff yet
    for i=1:length(mI);
        RU(i,:,:) = get_ptvariable('u',fname,[mI(i) mJ(i)],4,diry2);
        RV(i,:,:) = get_ptvariable('v',fname,[mI(i) mJ(i)],4,diry2);
        RT(i,:,:) = get_ptvariable('temp',fname,[mI(i) mJ(i)],4,diry2);
        RS(i,:,:) = get_ptvariable('salt',fname,[mI(i) mJ(i)],4,diry2);                        
  %      RH(i,:) = get_ptvariable('Hwave',[mI(i) mJ(i)],3);
        RZ(i,:) = get_ptvariable('zeta',fname,[mI(i) mJ(i)],3,diry2);
    end
    disp('Still going, still')
    Atmp = ones(length(A),length(t2),length(RS(1,1,:)));
    for i=1:length(t2)
        for j=1:length(RSS(1,1,:))
            Atmp(:,i,j) = A;
        end
    end
    W = RU + sqrt(-1)*RV;
    W = W .* exp(sqrt(-1).*Atmp);
    RU = real(W); RV = imag(W);
end

if(~exist('S4N','var'))
    disp('Observations')
%    S4N = 'G:\WAIPAOA\PROCESSED_DATA\northS4.nc';
%    S4M = 'G:\WAIPAOA\PROCESSED_DATA\mouthS4.nc';
%    S4S = 'G:\WAIPAOA\PROCESSED_DATA\southS4.nc';
    S4N = '/data/users-tmp/abever/WAIPAOA/FIELD_DATA/OCEAN_DATA/S4/northS4.nc';
    S4M = '/data/users-tmp/abever/WAIPAOA/FIELD_DATA/OCEAN_DATA/S4/mouthS4.nc';
    S4S = '/data/users-tmp/abever/WAIPAOA/FIELD_DATA/OCEAN_DATA/S4/southS4.nc';


    spd=nc_varget(S4N,'speed_ba');
    dir=nc_varget(S4N,'dir_ba');
    [un,vn]=pol2cart(dir,spd);
    un = un./100; vn = vn./100;
    timen = nc_varget(S4N,'date');
    tn = nc_varget(S4N,'temp_ba');
    sn = nc_varget(S4N,'sal_ba');
    hn = nc_varget(S4N,'Hmo');
    dn = nc_varget(S4N,'depth_ba');

    spd=nc_varget(S4M,'speed_ba');
    dir=nc_varget(S4M,'dir_ba');
    [um,vm]=pol2cart(dir,spd);
    um = um./100; vm = vm./100;
    timem = nc_varget(S4M,'date');
    tm = nc_varget(S4M,'temp_ba');
    sm = nc_varget(S4M,'sal_ba');
    hm = nc_varget(S4M,'Hmo');
    dm = nc_varget(S4M,'depth_ba');

    spd=nc_varget(S4S,'speed_ba');
    dir=nc_varget(S4S,'dir_ba');
    [us,vs]=pol2cart(dir,spd);
    us = us./100; vs = vs./100;
    times = nc_varget(S4S,'date');
    ts = nc_varget(S4S,'temp_ba');
    ss = nc_varget(S4S,'sal_ba');
    hs = nc_varget(S4S,'Hmo');
    ds = nc_varget(S4S,'depth_ba');

    p = isnan(dn);
    dn = dn - mean(dn(~p));
    p = isnan(dm);
    dm = dm - mean(dm(~p));
    p = isnan(ds);
    ds = ds - mean(ds(~p));

end

if(~exist('RUP','var'));
    disp('Generating progressive vectors')
%    fname = 'G:\WAIPAOA\FORCING\INP_3AUG\modern_test_grid.nc';
    fname = '/data/users-tmp/abever/MODELS/WAIPAOA2/FORCING/INP_3AUG/modern_test_grid.nc';
    X = nc_varget(fname,'x_rho');
    Y = nc_varget(fname,'y_rho');
    h = nc_varget(fname,'h');
    h(~mask)=nan;

    unx = cumsum((un*3600));   vnx = cumsum((vn*3600));
    umx = cumsum((um*3600));   vmx = cumsum((vm*3600));
    usx = cumsum((us*3600));   vsx = cumsum((vs*3600));

    unx = unx + X(mI(3),mJ(3));     vnx = vnx + Y(mI(3),mJ(3));
    umx = umx + X(mI(2),mJ(2));     vmx = vmx + Y(mI(2),mJ(2));
    usx = usx + X(mI(1),mJ(1));     vsx = vsx + Y(mI(1),mJ(1));
    
    RUP = cumsum((RU*3600*3),2);   RVP = cumsum((RV*3600*3),2);
    RSUP = cumsum((RSU*3600*3),2);   RSVP = cumsum((RSV*3600*3),2);

    if(~cv2ll)
        for i=1:length(mI);
            RUP(i,:,:) = RUP(i,:,:) + X(mI(i),mJ(i)); %RUP(i,:,:) = RUP(i,:,:) - max(max(max(RUP(i,:,:))));
            RVP(i,:,:) = RVP(i,:,:) + Y(mI(i),mJ(i)); %RVP(i,:,:) = RVP(i,:,:) - max(max(max(RVP(i,:,:))));
            RSUP(i,:,:) = RSUP(i,:,:) + X(mI(i),mJ(i)); %RSUP(i,:,:) = RSUP(i,:,:) - max(max(max(RSUP(i,:,:))));
            RSVP(i,:,:) = RSVP(i,:,:) + Y(mI(i),mJ(i)); %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)
        RUPL = ones(size(RUP))*nan;    RVPL = ones(size(RVP))*nan;
        RSUPL = ones(size(RSUP))*nan;    RSVPL = ones(size(RSVP))*nan;
        
        for i=1:length(N);
            [RUPL(i,:,:),RVPL(i,:,:)] = m_xy2ll(RUP(i,:,:)./Sc+E(i),RVP(i,:,:)./Sc+N(i));
            [RSUPL(i,:,:),RSVPL(i,:,:)] = m_xy2ll(RSUP(i,:,:)./Sc+E(i),RSVP(i,:,:)./Sc+N(i));
        end
        load /data/users-tmp/abever/MODELS/WAIPAOA2/nz_coast.mat
%     [coast(:,1) coast(:,2)] = m_ll2xy(coast(:,1),coast(:,2));
    end
end


if(lp)
    if(~exist('LP','var'));
        LP=[];
        ctoff = 1/33;
        for i=1:length(RU(:,1,1))
            for j=1:length(RU(1,1,:));
                RU(i,:,j) = lpfilt(squeeze(RU(i,:,j))',3,ctoff);
                RV(i,:,j) = lpfilt(squeeze(RV(i,:,j))',3,ctoff);
            end
                RZ(i,:) = lpfilt(squeeze(RZ(i,:))',3,ctoff);
        end
        for i=1:length(RSU(:,1,1))
            for j=1:length(RSU(1,1,:));
                RSU(i,:,j) = lpfilt(squeeze(RSU(i,:,j))',3,ctoff);
                RSV(i,:,j) = lpfilt(squeeze(RSV(i,:,j))',3,ctoff);
            end
            RSZ(i,:) = lpfilt(squeeze(RSZ(i,:))',3,ctoff);
        end
    p=isnan(us);
    us(~p) = lpfilt(us(~p),1,ctoff);
    p=isnan(um);
    um(~p) = lpfilt(um(~p),1,ctoff);
    p=isnan(un);
    un(~p) = lpfilt(un(~p),1,ctoff);
    p=isnan(vs);
    vs(~p) = lpfilt(vs(~p),1,ctoff);
    p=isnan(vm);
    vm(~p) = lpfilt(vm(~p),1,ctoff);
    p=isnan(vn);
    vn(~p) = lpfilt(vn(~p),1,ctoff);
    p=isnan(ds);
    ds(~p) = lpfilt(ds(~p),1,ctoff);
    p=isnan(dm);
    dm(~p) = lpfilt(dm(~p),1,ctoff);
    p=isnan(dn);
    dn(~p) = lpfilt(dn(~p),1,ctoff);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Plotting
sethelvB
n = datenum(gregorian(t));
% n2 = ones(length(n),20) * nan;
% for i=1:20
%     n2(:,i) = n;
%     
% end
n2 = datenum(gregorian(t2));
%n2 = ones(length(n),20) * nan;


% order of I,J's = south,mouth,north,ADCP

%  currents and waves figure2
for i = 1:length(mI)
figure

y0 = zeros(size(n));
subplot(5,1,1,'align');
P = plot(n2,squeeze(RU(i,:,zl)),'k',n,squeeze(RSU(i,:,zl)),'b');
set(P,'linewidth',2)
axis tight
hold on
if(i==1)
    P = plot(times,us,'r');
elseif(i==2);
    P = plot(timem,um,'r');
elseif(i==3);
    P = plot(timen,un,'r');
end
if(i==1 || i==2 || 1==3)
set(P,'linewidth',2)
end
plot(n,y0,'k:')
axis tight
set(gca,'Xlim',[min(n) max(n)])
datetick('x',6,'keeplimits')
% tmp_title = [run1,' and ',run2,sprintf(' Comparison, PT = %d',i);
% title(tmp_title)
ylabel('U Speed, m s^-^1')

subplot(5,1,2,'align');
P = plot(n2,squeeze(RV(i,:,zl)),'k',n,squeeze(RSV(i,:,zl)),'b');
set(P,'linewidth',2)
axis tight
hold on
if(i==1)
    P = plot(times,vs,'r');
elseif(i==2);
    P = plot(timem,vm,'r');
elseif(i==3);
    P = plot(timen,vn,'r');
end
if(i==1 || i==2 || 1==3)
set(P,'linewidth',2)
end
plot(n,y0,'k:')
axis tight
set(gca,'Xlim',[min(n) max(n)])
datetick('x',6,'keeplimits')
ylabel('V Speed, m s^-^1')

subplot(5,1,3,'align');
P = plot(n2,squeeze(RT(i,:,zl)),'k',n,squeeze(RST(i,:,zl)),'b');
set(P,'linewidth',2)
axis tight
hold on
if(i==1)
    P = plot(times,ts,'r');
elseif(i==2);
    P = plot(timem,tm,'r');
elseif(i==3);
    P = plot(timen,tn,'r');
end
if(i==1 || i==2 || 1==3)
set(P,'linewidth',2)
end
plot(n,y0,'k:')
axis tight
set(gca,'Xlim',[min(n) max(n)])
datetick('x',6,'keeplimits')
ylabel('Temperature')

subplot(5,1,4,'align');
P = plot(n2,squeeze(RS(i,:,zl)),'k',n,squeeze(RSS(i,:,zl)),'b');
set(P,'linewidth',2)
hold on
if(i==1)
    P = plot(times,ss,'r');
elseif(i==2);
    P = plot(timem,sm,'r');
elseif(i==3);
    P = plot(timen,sn,'r');
end
if(i==1 || i==2 || 1==3)
set(P,'linewidth',2)
end
plot(n,y0,'k:')
axis tight
set(gca,'Xlim',[min(n) max(n)])
datetick('x',6,'keeplimits')
ylabel('Salinity')
legend(run2,run1,3)

subplot(5,1,5,'align');
P = plot(n2,squeeze(RH(i,:)),'k',n,squeeze(RSH(i,:)),'b');
set(P,'linewidth',2)
axis tight
hold on
if(i==1)
    P = plot(times,hs,'r');
elseif(i==2);
    P = plot(timem,hm,'r');
elseif(i==3);
    P = plot(timen,hn,'r');
end
if(i==1 || i==2 || 1==3)
set(P,'linewidth',2)
end
axis tight
set(gca,'Xlim',[min(n) max(n)])
datetick('x',6,'keeplimits')
xlabel('Date')
ylabel('Wave Height')

set(gcf,'pos',[17   146   950   800])
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% progressive vector figure
tmp=sprintf('Scale Factor = %d',Sc);
figure
%fillseg(coast)
hold on
if(~cv2ll)
    for i=1:length(mI)
        P = plot(squeeze(RUP(i,:,zl)),RVP(i,:,zl),'k');
        set(P,'linewidth',2)
        plot(squeeze(RUP(i,1,zl)),RVP(i,1,zl),'r<','markersize',15);

        P = plot(squeeze(RSUP(i,:,zl)),RSVP(i,:,zl),'b');
        set(P,'linewidth',2)
        plot(squeeze(RSUP(i,1,zl)),RSVP(i,1,zl),'c>','markersize',15);
    end
    title('Progressive Vectors')
else
    fillseg(coast)
    for i=1:length(mI)
        P = plot(squeeze(RUPL(i,:,zl)),RVPL(i,:,zl),'k');
        set(P,'linewidth',2)
        plot(squeeze(RUPL(i,1,zl)),RVPL(i,1,zl),'r<','markersize',15);

        P = plot(squeeze(RSUPL(i,:,zl)),RSVPL(i,:,zl),'b');
        set(P,'linewidth',2)
        plot(squeeze(RSUPL(i,1,zl)),RSVPL(i,1,zl),'c>','markersize',15);
    end
    title({'Progressive Vectors',tmp})
end    
legend(run2,run2,run1,run1)

P=plot(rlon(end,:),rlat(end,:));
set(P,'linewidth',2,'color','k');
P=plot(rlon(:,end),rlat(:,end));
set(P,'linewidth',2,'color','k');
dasp(-38)
axis equal
grid on
set(gcf,'pos',[100 100 1100 800])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sea surface height figure
figure
set(gcf,'pos',[100 100 1100 800])
subplot(4,1,1,'align');
P = plot(n2,squeeze(RZ(1,:)),'k',n,squeeze(RSZ(1,:)),'b',times,ds,'r');
set(P,'linewidth',2);
axis tight
%set(gca,'Xlim',[datenum([2006 4 23 0 0 0]) datenum([2006 5 21 0 0 0])]...
%    ,'Ylim',[-1 1]);
datetick('x',6,'keeplimits')
title('Sea Surface Height (zeta), South')
legend(run2,run1,'Observations');
grid on

subplot(4,1,2,'align');
P = plot(n2,squeeze(RZ(2,:)),'k',n,squeeze(RSZ(2,:)),'b',timem,dm,'r');
set(P,'linewidth',2);
axis tight
%set(gca,'Xlim',[datenum([2006 4 23 0 0 0]) datenum([2006 5 21 0 0 0])]...
%    ,'Ylim',[-1 1]);
datetick('x',6,'keeplimits')
title('Mouth')
ylabel('Sea Surface Height, m')
grid on

subplot(4,1,3,'align');
P = plot(n2,squeeze(RZ(3,:)),'k',n,squeeze(RSZ(3,:)),'b',timen,dn,'r');
set(P,'linewidth',2);
axis tight
%set(gca,'Xlim',[datenum([2006 4 23 0 0 0]) datenum([2006 5 21 0 0 0])]...
%    ,'Ylim',[-1 1]);
datetick('x',6,'keeplimits')
title('North')
grid on

subplot(4,1,4,'align');
P = plot(n2,squeeze(RZ(5,:)),'k',n,squeeze(RSZ(5,:)),'b');
set(P,'linewidth',2);
axis tight
%set(gca,'Xlim',[datenum([2006 4 23 0 0 0]) datenum([2006 5 21 0 0 0])]...
%    ,'Ylim',[-1 1]);
datetick('x',6,'keeplimits')
title('Near Eastern OB.')
xlabel('Date')
grid on

