% pom4cyc
%       read in POM's (T,U,V) binary data, and plot them.
%   
%-------
% input
%-------
%       T.YYMMDDHH - temperature volume 
%       U.yymmddhh
%       V.yymmddhh
%
%-------
% output
%-------
%       *.bmp
%
%--------
% history
%--------
%   7/8/2008   Hyun-Sook Kim
%--------------------------------------------------------------------------
clear all; close all; %pack

%load ClomClam;     % coastline
load('eastcoastmask_new.mat') %coastline

%--- Run information:
version='H050';
Hname='irene09';
ID='09';
yyyy='2011';
mm='08';
dd='26';
hh='00';
ds0=[yyyy,mm,dd,hh];

cycid=0;       % cycle

dn=datenum(ds0,'yyyymmddHH')+cycid/24;
ds=datestr(dn,'ddHH');

Ddir    = '/Volumes/boardwalk/seroka/data/HWRF/POM/';

%irene09l.2011082600.T.11082700 

%------ Working directory: compile data into a big data file
%Ddir=['Z:\',version,'\',Hname,'\',ds0,'\cyc',num2str(cycid,'%03g'),'-',ds,'\ocean\'];
Ddir    = '/Volumes/boardwalk/seroka/data/HWRF/POM/';

fprefix='irene09l';


ttl{1}=['Hurricane',' ',Hname,' ',  'POM Simulation: Initial time: ',datestr(dn,'yyyymmddHH'),'Z '];
ttl{2}=['Surface T:     Cycle',num2str(cycid,'%03g'),' run'];
%tit1=strcat('Temperature and Currents at 75m Depth:- Phase4:-: ',fcsthr,'h:  run');

%----- Track:
% HWRF or GFS storm track file
% frc='HWRF';
% switch lower(frc)
%     case {'hwrf'}
%         afile=[Ddir,fprefix,'.',datestr(dn,'yyyymmddHH'),'.trak.hwrf.atcfunix'];
%         afile=[Ddir,'trak.hwrf.atcfunix.',datestr(dn,'yyyymmddHH')];
%         dlook=dir(afile);
%         tfile=dlook.name;
%         hsk=caseread([Ddir,tfile]);            clear dlook tfile
%         hhr=str2num(hsk(:,31:33));
%         knt=str2num(hsk(:,65:66));      % reference knot
%         % find index for reference wind speed=34kt
%         ix=find(knt==34);
%         hsk=hsk(ix,:);
%         hx=-1*str2num(hsk(:,43:45))/10;
%         hy=str2num(hsk(:,36:38))/10; 
%         hhr=hhr(ix);
%         hdn=datenum(hsk(:,9:18),'yyyymmddHH')+hhr/24;
%         vkt=str2num(hsk(:,49:51));              % max. wind speed (kts)
%         pmbar=str2num(hsk(:,53:57));            % minimum Pressure (mbar)
%         Rads=str2num(hsk(:,74:95)) *1.8;        % size [nm] to [km]
%         Rsize=[];
%         for n=1:length(ix)
%             tmp=Rads(n,:);
%             mn=mean(tmp(tmp~=0));
%             Rsize(n)=mn;
%             clear tmp mn
%         end
%         clear Rads hhr B ix jx hsk
%     case {'gfs'}
%         nhc=load('G:\PROJECT07\OFS\GFStrack\Rita.best.track.gfs.txt');
%         hx=nhc(:,3); hy=nhc(:,2);
%         ymdh=num2str(nhc(:,1));
%         hdn=datenum(ymdh,'yyyymmddHH');
% end

% Rita Property: [year month day lat lon pres wind (kt) size(km)]
%
% 1  2  3   4   5        6       7      8    9
%-------------------
% dv=datevec(hdn);
% Hurricane=[dv(:,1:4) hy hx pmbar vkt Rsize(:)];
% 
% dL=fnddist(hx(2:end),hy(2:end),hx(1:end-1),hy(1:end-1));    % [km]
% Ut=dL./diff(hdn)*(1000/24/3600);      % Storm Translation Speed [m/s]
% clear dL
% 
% %--- best HWRF hurricane file
% nhc=caseread(['G:\PROJECT07\TPCtrack\y',yyyy,'\bal',ID,yyyy,'.dat']);
% indx=strmatch('  34,',nhc(:,63:67),'exact');        % wind speed threshold (kt)
% nhc=nhc(indx,:);    clear indx
% nhy=str2num(nhc(:,35:38))/10;
% nhx=(-1)*str2num(nhc(:,41:45))/10;
% npmin=str2num(nhc(:,53:57));
% nvmax=str2num(nhc(:,48:51));
% ndn=datenum(nhc(:,8:18),'yyyymmddHH');
% 
% best=[ndn nhx nhy npmin nvmax];
% clear nhc nhx nhy npmin nvmax ndn

%--------------------------------------------------------------------------
% prep:
%       model parameters
%--------------------------------------------------------------------------
load /Users/seroka/Research/data/HWRF/pom_config.mat

%--------------------------------------------------------------------------
% main:
%--------------------------------------------------------------------------
% model temperature field (3D)
Tfiles=dir([Ddir,fprefix,'.*.T.*']);

% model (east,north) velocity field (3D)
Ufiles=dir([Ddir,fprefix,'.*.U.*']);
Vfiles=dir([Ddir,fprefix,'.*.V.*']);

Nrcd=length(Tfiles);
%%
% work space
% Tdat=repmat(nan,[Nrcd,klv,jm,im]);
Tdat=repmat(nan,[Nrcd,1,jm,im]);
pdn=[];       % date number for POM output
for K=1:Nrcd
    
    ftmp=[Ddir,Tfiles(K).name];     % temp
    
    %<---- build model datanumber
    ymdh=ftmp(end-7:end);
    dn=datenum(ymdh,'yymmddHH');
    pdn(K)=dn;
    clear dn ymdh
    %-----> 
    
    for L=1:klv
        fid=fopen(ftmp,'r','b');  
        n=fread(fid,1,'int32');

        skip=(L-1)*jm*im*4;
        
        fseek(fid,skip,'cof');
        dtmp=fread(fid,[im,jm],'float');
        dtmp=dtmp';
        dtmp=dtmp+Toff;
        %if L==1
            Tdat(K,L,:,:)=dtmp;
        %end
        clear n skip dtmp
        
        fclose(fid);
    end
end
fclose all;

lnd=find(Tdat<-80);     % masking for NaNs including land
Tdat(lnd)=nan;     % land masking

% velocity x-component
Udat=repmat(nan,[Nrcd,1,jm,im]);
for K=1:Nrcd
    
    ftmp=[Ddir,Ufiles(K).name];     % temp
    for L=1:klv
        fid=fopen(ftmp,'r','b');  
        n=fread(fid,1,'int32');
        skip=(L-1)*jm*im*4;
        
        fseek(fid,skip,'cof');
        dtmp=fread(fid,[im,jm],'float');
        dtmp=dtmp';
        %if L==1
            Udat(K,L,:,:)=dtmp;
        %end
        clear n skip dtmp
        fclose(fid);
    end
end
fclose all;
Udat(lnd)=nan;

% velocity y-component
Vdat=repmat(nan,[Nrcd,1,jm,im]);
for K=1:Nrcd
    
    ftmp=[Ddir,Vfiles(K).name];     % temp
    for L=1:klv
        fid=fopen(ftmp,'r','b');  
        n=fread(fid,1,'int32');
        skip=(L-1)*jm*im*4;
        
        fseek(fid,skip,'cof');
        dtmp=fread(fid,[im,jm],'float');
        dtmp=dtmp';
        %if L==1;        % sst
            Vdat(K,L,:,:)=dtmp;
        %end
        clear n skip dtmp
        fclose(fid);
    end
end
fclose all;
Vdat(lnd)=nan;

%--- final stage ----------------------------------------------------------
disp('                 .........  rd_pomout.m done!')

sst=squeeze(Tdat(:,1:klv:end,:,:));     % model sea surface temperature

%--------------------------------------------------------------------------
% diagnostic plot
%-----------------
ref=4;      % reference velocity: 5 m/s
Udat(:,:,175:177,50:52)=sqrt(ref);      % reference velocity x-component
Vdat(:,:,175:177,50:52)=sqrt(ref);

gskip=3;    % subsampling velocity vectors just for plotting only
js=1:gskip:length(ys);
is=1:gskip:length(xs);
%%    
bc=[0.5 0.2 0.8];     % color for best track data
for K=6:6%Nrcd
    
    % a point on HWRF-POM track
    %hxi=interp1(hdn,hx,pdn(K),'linear');
    %hyi=interp1(hdn,hy,pdn(K),'linear');

    % a point on the best track
    %nhxi=interp1(best(:,1),best(:,2),pdn(K),'linear');
    %nhyi=interp1(best(:,1),best(:,3),pdn(K),'linear');
    
    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
    
    load('HWRF_POMdata_irene.mat')
    
    [xsgrid,ysgrid] = meshgrid(xs,ys);
    
    depthnew    = interp2(loni,lati,-depthi,romslon,romslat);
    
    dkm   = nan(225,254,82);
    dkm_min = nan(82,1);
    crosssectionlat = dkm_min;
    crosssectionlon = dkm_min;
    index_y = dkm_min;
    index_x = dkm_min;
    
    Tdatnew = nan(23,82);
    for m = 1:length(romslat)
        for i = 1:length(xs)
            for j = 1:length(ys)
                dkm(j,i,m) = lldistkm([romslat(m) romslon(m)],[ysgrid(j,i) xsgrid(j,i)]);
                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) = ysgrid(r,c);
        crosssectionlon(m) = xsgrid(r,c);
        index_y(m) = r;
        index_x(m) = c;
        Tdatnew(:,m) = Tdat(K,:,index_y(m),index_x(m));
        m = m+1;
    end
        
    %figure('position',[34   500   716   420]);
    close all
    figure;
    set(gcf,'Paperposition',[0 0 6 6]);
    %colorhsk(xs,ys,squeeze(Tdat(K,1,:,:))); shading flat; caxis([15 35]);
% %     pcolor(xs,ys,squeeze(Tdat(K,1,:,:))); shading flat;caxis([14 30]);
    hold on
    %%slicexpom=[xs(125), xs(140)];
    %%sliceypom=[ys(175), ys(175)];
    %%plot(slicexpom,sliceypom,'k')
    %%%slicexroms=[romslon(1), romslon(end)];
    %%%sliceyroms=[romslat(1), romslat(end)];
    %%%plot(slicexroms,sliceyroms,'b')
    %pcolor(xs(125:end),zs,squeeze(Tdat(K,:,175,125:end))); shading flat; caxis([9 24]); %xslice
    pcolor(xs(index_x),zs,Tdatnew); shading flat; caxis([9 30]); %jslice (from ROMS)
    set(gca,'YDir','reverse')
    hold on
    
    plot(romslon,depthnew,'k','Linewidth',2)

    %plot(loni(968,797:end),-depthi(968,797:end),'k','Linewidth',2)
        
    %colormap(jet(40)); hc=colorbar; set(hc,'tickdir','out');
% %     axis([-77 -68 34 43])
    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');
    %quiver(xs(is),ys(js),squeeze(Udat(K,1,js,is)),squeeze(Vdat(K,1,js,is)),5,'k');
    grid
    %[cs,h]=contour(xs,ys,squeeze(Tdat(K,1,:,:)),[26 26],'r');
    %set(h,'linewidth',1.3); clabel(cs,h);
    % HWRF-POM track
    %plot(hx,hy,'-m','linewidth',1.5);
    %plot(hxi,hyi,'oy','linewidth',1.5);
%--- BEST track
    %plot(best(:,2),best(:,3),'-x','linewidth',1.8,'color',bc);
    %plot(nhxi,nhyi,'oy','linewidth',1.8,'markersize',7);

    %plot coastline
    %plot(lon,lat,'k')

    set(gca,'tickdir','out','fontsize',12);
    %ylabel('Latitude');
    %xlabel('Longitude');
    
%     cbfreeze(hc)
%     freezeColors
%     MAB_no_contour_generator(1)
%     grid on;
%     set(gca,'layer','top')
    
    title(ttl,'fontsize',12);
    
    %text(xs(50),ys(170),[num2str(ref),' m/s'],'fontweight','bold','color','b');

    text(0.01,0.05,[num2str((K-1)*24),'-hr (',datestr(pdn(K),'yyyy/mm/dd HHZ'),')'],...
        'fontsize',12,'fontweight','bold','color','k','units','normalized');
% %     text(0.01,0.95,[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');
    
    pause(0.1);
    
% %     project_mercator
    wysiwyg
    %set(gcf,'paperpositionmode','auto');
    %eval(['print -dpng /Users/seroka/Research/data/HWRF/POM_TUVsrf_cyc',num2str(cycid,'%03g'),'_',datestr(pdn(K),'yymmddHH')]);
% %     eval(['print -dpng /Users/seroka/Research/data/HWRF/POM_Tsrf_cyc',num2str(cycid,'%03g'),'_',datestr(pdn(K),'yymmddHH'),'_3']);
    eval(['print -dpng /Users/seroka/Research/data/HWRF/POM_T_xslice_cyc',num2str(cycid,'%03g'),'_',datestr(pdn(K),'yymmddHH'),'_romsslice_30C_zoomoutfull']);

    close
end

%savefile    = strcat('/Users/seroka/Research/data/HWRF/HWRF_POMdata_irene.mat');
%save(savefile,'Tdat','xs','ys','zs')

%%
%--------------------------------------------------------------------------
%  save output
%--------------
ofname1=[Ddir,'\hwrfpom',ds0,'Tdat'];
eval(['save ',ofname1,' Hurricane Ut Tdat hdn im jm klv xs ys zs'])

ofname2=[Ddir,'\hwrfpom',ds0,'Udat'];
eval(['save ',ofname2,' Hurricane Ut Udat hdn im jm klv xs ys zs'])

ofname3=[Ddir,'\hwrfpom',ds0,'Vdat'];
eval(['save ',ofname3,' Hurricane Ut Vdat pdn hdn im jm klv xs ys zs'])



