%plot HYCOM angular cross section from NJ coast to shelf break
%to compare to POM cross section3

clear all; close all; %pack

%load('eastcoastmask_new.mat') %coastline

load('Users/seroka/Research/data/HWRF/ncep12_04_hur_beast_grid.mat');
%load('HWRF_HYCOMdata_irene_aug26_00Z.mat') 

hycom_ncfile = ('/Volumes/boardwalk/seroka/data/HWRF/HYCOM/hwrf_IRENE.09l.ofs_hwrf04.t00z.F120.3d.nc');
layth = nc_varget(hycom_ncfile,'LAYTH_120_HYBY');
hycomtemp     = nc_varget(hycom_ncfile,'WTMPC_120_HYBY');
hycomsst      = nc_varget(hycom_ncfile,'WTMPC_120_240');  %??
hycomlon      = blon;
hycomlat      = blat;
hycombathy    = bdh;

hycomdepth = nan(26,683,403);
for i = 1:size(hycomlon,1)
    for j = 1:size(hycomlon,2)
        for k = 1:size(hycomtemp,1)
            if k == 1
                hycomdepth(k,j,i) = layth(k,j,i);
            else
                hycomdepth(k,j,i) = layth(k,j,i) + hycomdepth(k-1,j,i);
            end
        end
    end
end

load('roms_jview_temp_greg.mat') 
romslat = squeeze(lat_i(1,:));
romslon = squeeze(lon_i(1,:));

load bathymetry_USeastcoast
baddepth = find(depthi == 999999);
if isempty(baddepth) ~= 1
    depthi(baddepth) = NaN;
end

depthnew    = interp2(loni,lati,-depthi,romslon,romslat);

%[xsgrid,ysgrid] = meshgrid(xs,ys);  %don't need for HYCOM- hycomlon and hycomlat are already gridded

dkm   = nan(683,403,82);
dkm_min = nan(82,1);
crosssectionlat = dkm_min;
crosssectionlon = dkm_min;
index_y = dkm_min;
index_x = dkm_min;

Tdatnew = nan(26,82);
hycomlon_new = nan(82,1);
hycombathynew = nan(82,1);
hycomdepth_new = nan(26,1);
for m = 1:length(romslat)
    for i = 1:size(hycomlon,1)
        for j = 1:size(hycomlon,2)
            dkm(j,i,m) = lldistkm([romslat(m) romslon(m)],[hycomlat(i,j) hycomlon(i,j)]);
            j = j+1;
        end
        i = i+1;
    end
    [r,c]   = find(dkm(:,:,m)==nanmin(nanmin(dkm(:,:,m))));   
    dkm_min(m) = nanmin(nanmin(dkm(:,:,m)));
    crosssectionlat(m) = hycomlat(r,c);
    crosssectionlon(m) = hycomlon(r,c);
    index_y(m) = r;
    index_x(m) = c;
    Tdatnew(:,m) = hycomtemp(:,index_y(m),index_x(m));
    hycomlon_new(m,1) = hycomlon(index_x(m),index_y(m));
    hycombathy_new(m,1) = hycombathy(index_x(m),index_y(m));
    hycomdepth_new(:,m) = hycomdepth(:,index_y(m),index_x(m));
    m = m+1;
end

hycomlon_new2 = repmat(hycomlon_new,1,26);

%%
close all
figure;
set(gcf,'Paperposition',[0 0 6 6]);
%pcolor(xs(125:140),zs,squeeze(Tdat(K,:,175,125:140))); shading flat; caxis([9 24]); %xslice
pcolor(hycomlon_new2',hycomdepth_new,Tdatnew); shading flat; caxis([9 30]); %jslice (from ROMS)
set(gca,'YDir','reverse')

hold on

plot(romslon,depthnew,'k','Linewidth',2)
plot(romslon,hycombathy_new,'k--','Linewidth',2)

xlim([-74.75 -72])
ylim([0 200])
%xlim([-74.75 -71])
%ylim([0 1000])
%ylim([0 4000])
%xlim([-74.75 -71.5])
colormap(jet); hc=colorbar; set(hc,'tickdir','out');
grid

set(gca,'tickdir','out','fontsize',12);

ttl{1}=['Hurricane irene09 HYCOM Simulation: Initial time: 2011082600Z'];
ttl{2}=['T:     Cycle000 run'];

title(ttl,'fontsize',12);

text(0.01,0.05,'120-hr (2011/08/31 00Z)',...
    'fontsize',12,'fontweight','bold','color','k','units','normalized');
%text(0.01,0.1,[num2str((K-1)*24),'-hr (',datestr(pdn(K),'yyyy/mm/dd HHZ'),')'],...
%    'fontsize',12,'fontweight','bold','color','k','units','normalized');
%text(0.01,0.05,['Xslice at lat=',num2str(ys(175))],...
%    'fontsize',12,'fontweight','bold','color','k','units','normalized');

wysiwyg
eval('print -dpng /Users/seroka/Research/data/HWRF/HYCOM_T_xslice_cyc000_11083100_romsslice_30C.png');
