
addpath ~/AARONS_TOOLS/PLOTTING/
addpath ~/USGS_toolbox/ecomviz-sea-mat


nxx=1250;    nyy=800;
Maxis=[-82.2937  -81.7934   26.3304   27.0197];
pax=[177.85 178.1 -38.81 -38.56];

coupled = 1;  % to get wave height from history file

% getting todays date
today = date;

%mdir = '/data/users-tmp/abever/MODELS/WAIPAOA2/SVN30/';
%mdir = '/data/users-tmp/abever/MODELS/WAIPAOA2/OUTPUT/FEB2009_B/'
%mdir = '../OUTstable/'
mdir = '../OUTWET8/'
mgrid = '../FRC2008/swf_grid8.nc';
wind_file ='../FRC2008/swf_wind2008.nc';
river_file = '../FRC2008/swf_river2008.nc';
wave_file = '';

%bkgrn = get_logo('/data/users-tmp/abever/MODELS/WAIPAOA2/','North Island.png');
%bkgrn = get_logo('/data/users-tmp/abever/MODELS/WAIPAOA2/','new_zealand1.png');
%bkgrn = get_logo('/data/users-tmp/abever/MODELS/WAIPAOA2/','NI East coast.png');
% bkgrn = get_logo('/data/users-tmp/abever/MODELS/WAIPAOA2/','Poverty_bay3D.png');
% nsflogo=get_logo('/data/users-tmp/abever/MODELS/WAIPAOA2/','MARGINSsmall.jpg');
vlogo=get_logo('./','FGCU logo horizontal (bk).jpg');;
% onrlogo = get_logo('/data/users-tmp/abever/MODELS/WAIPAOA2/','onr_logo_small.jpg');;
%     onrlogo = onrlogo(1:2:end,1:2:end,:);
%  
%load([mdir,'sed_trans_paths_wet.mat']);  
load(['sed_trans_paths_wet.mat']);  %%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% loading in the already generated profile .mat file to show sediment
% transport in or out of the mouth of the bay. Also calculating the
% sediment transport based on if the profile runs along I or J lines.
% load([mdir,'profile.mat']);
% if(rect(3) == rect(4));
%     ST = SED.*U;
% elseif(rect(1) == rect(2));
%     ST = SED.*V;
% end
% ST(ST>1e10) = nan;
% 
% for i=1:20
%     D2(i,:) = D;
% end
% 
% D(end+1) = D(end);  D(end+1) = D(1);   % for filling in seabed.
% ZS = Z;
% H(end+1) = -100;    H(end+1) = -100;   % for filling in seabed.
% bed_elev=H; clear H
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% setting the axes positions and sizes
ax8=[0 0 1 1];
ax1=[0 .27 .25 .65];
ax2=[.25 .27 .25 .65];
ax3=[.5 .27 .25 .65];
ax4=[.75 .27 .25 .65];
ax5=[.5 .025 .23 .23];
ax6=[.045 .91 .5 .075];
ax7=[.78 .13 .175 .175];
ax9=[.275 .1 .23 .23];
ax10=[.4 .125 .23 .23];
ax11=[.3583 .155 .23 .23];
ax12 = [.05 .13 .22 .22];
nsfl_ax=[.15 0 .16 .07];
vl_ax=[.01 0 .16  .07];
onrl_ax = [.85 0 .16  .07];


theta = 0:.01:2*pi;   % for plotting a circle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% getting the names of the history files and the  grid spesifics
MncdfF = dir([mdir,'calus2008_his*.nc']) % all of them
%MncdfF = MncdfF(1);  %subset for test

mrlon=nc_varget(mgrid,'lon_rho',[0 0],[-1 -1]);
mrlat=nc_varget(mgrid,'lat_rho',[0 0],[-1 -1]);
mulon=nc_varget(mgrid,'lon_u',[0 0],[-1 -1]);
mulat=nc_varget(mgrid,'lat_u',[0 0],[-1 -1]);
mvlon=nc_varget(mgrid,'lon_v',[0 0],[-1 -1]);
mvlat=nc_varget(mgrid,'lat_v',[0 0],[-1 -1]);
mmask=nc_varget(mgrid,'mask_rho',[0 0],[-1 -1]);
mangle=nc_varget(mgrid,'angle',[0 0],[-1 -1]);
mh=nc_varget(mgrid,'h',[0 0],[-1 -1]);
pm=1./nc_varget(mgrid,'pm',[119 0],[1 50]);  pm=cumsum(pm);
pm(end+1)=max(pm);  pm(end+1)=min(pm);    % for shading purposes
pm = pm./1000;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% wind file, getting the time indicators
wtime=nc_varget(wind_file,'wind_time',[0],[-1])+2440000;
%%%%%%opening the river file and getting the time indicators
rtime=nc_varget(river_file,'river_time');
rtime=rtime+2440000;
RD_F=nc_varget(river_file,'river_transport');
RD_F=sum(RD_F,2);
rn=datenum(gregorian(rtime));
%%%%%%%%%%%%%%%%%%%%%%%%opening wave file
if(coupled==0)
    wavenc=netcdf(wave_file);
    wavet=wavenc{'wave_time'}+2440000;
else
    wavet = rtime; % acts as a place holder
end

load swfl_coast.mat;
coast=[lon lat];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is because griddata is so much slower than interp2, 
% method copied from stuff Rich Signell had written
[nx,ny] = size(mulon);
[iiu,jju] = meshgrid(1:ny,1:nx);
ibu = griddata(mulon,mulat,iiu,mrlon,mrlat);
jbu = griddata(mulon,mulat,jju,mrlon,mrlat);
keepiu = find(isnan(ibu) == 1);  ibu(keepiu) = 1;
keepju = find(isnan(jbu) == 1);  jbu(keepju) = 1;

[nx,ny] = size(mvlon);
[iiv,jjv] = meshgrid(1:ny,1:nx);
ibv = griddata(mvlon,mvlat,iiv,mrlon,mrlat);
jbv = griddata(mvlon,mvlat,jjv,mrlon,mrlat);
keepiv = find(isnan(ibv) == 1);  ibv(keepiv) = 1;
keepjv = find(isnan(jbv) == 1);  jbv(keepjv) = 1;

clear nx ny x y


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T=-1;
ztmp = [1:1:12];
nfiles = length(MncdfF);
sethelv14BW
clear i
begin=1;
for jj = 1:nfiles
    
    tmphis = [mdir,MncdfF(jj).name];
    t=nc_varget(tmphis,'ocean_time');
    t=t./(3600*24)+2440000;

    for ii=1:length(t)
        T=T+1;
        
        k=ii+(jj-1)*40
  %%%%%%%%%%%%%%%%%getting wind,wave, river, and sediment transport index
        w_ind=tind(wtime,gregorian(t(ii)));
        r_ind=tind(rtime,gregorian(t(ii)));
        wave_ind=tind(wavet,gregorian(t(ii)));
  ti=tind(transport_time,gregorian(t(ii)));
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% depths and grid angle of the stations
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for vertical velocity plots
  z=depths2(tmphis,1,0,ii);
        
        
        z=permute(z,[3,2,1]);
        

        Z1=z(:,92,10);   
        XX=ones(size(Z1)).*0;  YY=ones(size(Z1)).*0;
        AS1=mangle(92,10);
        
        u=nc_varget(tmphis,'u',[ii-1 0 91 9],[1 20 1 1]);
        v=nc_varget(tmphis,'v',[ii-1 0 91 9],[1 20 1 1]);
        W=nc_varget(tmphis,'w',[ii-1 0 91 9],[1 20 1 1]);
        wsta=u+sqrt(-1).*v;    wsta=wsta.*exp(sqrt(-1).*AS1);
        u=real(wsta);   v=imag(wsta);
        
        Z2=z(:,110,20);   
        XX2=ones(size(Z1)).*0;  YY2=ones(size(Z1)).*0;
        AS2=mangle(110,20);
        
        u2=nc_varget(tmphis,'u',[ii-1 0 109 19],[1 20 1 1]);
        v2=nc_varget(tmphis,'v',[ii-1 0 109 19],[1 20 1 1]);
        W2=nc_varget(tmphis,'w',[ii-1 0 109 19],[1 20 1 1]);
        wsta=u2+sqrt(-1).*v2;    wsta=wsta.*exp(sqrt(-1).*AS2);
        u2=real(wsta);   v2=imag(wsta);

        Z3=z(:,115,10);   
        XX3=ones(size(Z1)).*0;  YY3=ones(size(Z1)).*0;
        AS3=mangle(115,10);
        
        u3=nc_varget(tmphis,'u',[ii-1 0 114 9],[1 20 1 1]);
        v3=nc_varget(tmphis,'v',[ii-1 0 114 9],[1 20 1 1]);
        W3=nc_varget(tmphis,'w',[ii-1 0 114 9],[1 20 1 1]);
        wsta=u3+sqrt(-1).*v3;    wsta=wsta.*exp(sqrt(-1).*AS3);
        u3=real(wsta);   v3=imag(wsta);      
  %%%%%%%%%%%%%salinity
        Ms=nc_varget(tmphis,'salt',[ii-1 19 0 0],[1 1 200 250]);
        Ms(~mmask)=nan;
  %%%%%%%%%%%%%%%%%currents
        mubar=nc_varget(tmphis,'ubar',[ii-1 0 0],[1 200 249]);
        mvbar=nc_varget(tmphis,'vbar',[ii-1 0 0],[1 199 250]);
        mu=nc_varget(tmphis,'u',[ii-1 19  0 0],[1 1 200 249]);
        mv=nc_varget(tmphis,'v',[ii-1 19  0 0],[1 1 199 250]);
        mub=nc_varget(tmphis,'u',[ii-1 0  0 0],[1 1 200 249]);
        mvb=nc_varget(tmphis,'v',[ii-1 0  0 0],[1 1 199 250]);
        
%         mubar=griddata(ulong,ulatg,mubar,rlong,rlatg);
%         mvbar=griddata(vlong,vlatg,mvbar,rlong,rlatg);
        
%         mu=griddata(ulong,ulatg,mu,rlong,rlatg);
%         mv=griddata(vlong,vlatg,mv,rlong,rlatg);
        
%         mub=griddata(ulong,ulatg,mub,rlong,rlatg);
%         mvb=griddata(vlong,vlatg,mvb,rlong,rlatg);
        mubar = interp2(iiu,jju,mubar,ibu,jbu);
        mvbar = interp2(iiv,jjv,mvbar,ibv,jbv);

        mu = interp2(iiu,jju,mu,ibu,jbu);
        mv = interp2(iiv,jjv,mv,ibv,jbv);
        
        mub = interp2(iiu,jju,mub,ibu,jbu);
        mvb = interp2(iiv,jjv,mvb,ibv,jbv);
             
        mw=mu+sqrt(-1).*mv;   mw=mw.*exp(sqrt(-1).*mangle);
        mwb=mub+sqrt(-1).*mvb;   mwb=mwb.*exp(sqrt(-1).*mangle);
        mwbar=mubar+sqrt(-1).*mvbar;   mwbar=mwbar.*exp(sqrt(-1).*mangle);
        
        mw(~mmask)=nan;      mw(find(mw>50))=nan;  
        mwb(~mmask)=nan;     mwb(find(mwb>50))=nan;
        mwbar(~mmask)=nan;   mwbar(find(mwbar>50))=nan;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Sediment Transport
%%%%%%%%%%%%%%%%%%% Just making the arrows didn't work well, some were HUGE
%%%%%%%%%%%%%%%%%%% and others non-existent.....making them log10 arrows...
        tmpu = squeeze(Ut(ti,:,:));  tmpu2 = tmpu./abs(tmpu);
        tmpv = squeeze(Vt(ti,:,:));  tmpv2 = tmpv./abs(tmpv);
        
        [tmpTH,tmpR] = cart2pol(tmpu,tmpv);
        %tmpR = log10(tmpR+.000001);
        %tmpR(tmpR<.05) = nan;
        tmpR(tmpR>20) = nan;
        tmpR(tmpR>1) = 1;
        tmpR(~mmask) = nan;
        p=isnan(tmpR);
        if(min(min(p))==1)
            tmpR(:,:) = .00000001;
        end
        [tmpu,tmpv] = pol2cart(tmpTH,tmpR);
        
%         tmpu(tmpu<-3) = nan;    tmpv(tmpv<-3) = nan;
%         tmpu(tmpu>10) = nan;    tmpv(tmpv>10) = nan;
        
        WT = tmpu+sqrt(-1).*tmpv;   WT=WT.*exp(sqrt(-1).*mangle);
        tmpmag = sqrt(real(WT).^2 + imag(WT).^2);
        WT(tmpmag<.005) = nan;  WT(1:10,1:10) = .1+sqrt(-1)*.1;
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%WAVES
        if(coupled==1)
            hsig = nc_varget(tmphis,'Hwave',[ii-1 0 0],[1 200 250]);
        else
            hsig=nc_varget(wave_file,'Hwave',[wave_ind-1 0 0],[1 200 250]);
        end
        hsig(~mmask)=nan;
   %%%%%%%%%%%%%%%%%%%%%%sediment
        MmDI=nc_varget(tmphis,'mud_01',[ii-1 0 0 0],[1 20 200 250]) + ...
            nc_varget(tmphis,'mud_02',[ii-1 0 0 0],[1 20 200 250]) + ...
            nc_varget(tmphis,'mud_03',[ii-1 0 0 0],[1 20 200 250]);
        
        Mm=squeeze(MmDI(20,:,:));
        
        mZ=depths([mdir,MncdfF(jj).name],mgrid,5,ii);
        mZ=permute(mZ,[ 3 1 2]);    
        mZ=diff(mZ);

        MmDI=MmDI.*mZ;      MmDI=squeeze(sum(MmDI));
        MmDI(~mmask)=nan;
    %%%%%%%%%%%%%%deposition
    if(exist('begin','var'))
        Mdep1=nc_varget(tmphis,'mudmass_01',[ii-1 0 0 0],[1 8 200 250]) +...
            nc_varget(tmphis,'mudmass_02',[ii-1 0 0 0],[1 8 200 250]);
            nc_varget(tmphis,'mudmass_03',[ii-1 0 0 0],[1 8 200 250]);
        Mdep1=squeeze(sum(Mdep1));
        Mdep1(~mmask)=nan;
    end
        Mdep=nc_varget(tmphis,'mudmass_01',[ii-1 0 0 0],[1 8 200 250]) +...
            nc_varget(tmphis,'mudmass_02',[ii-1 0 0 0],[1 8 200 250]);
            nc_varget(tmphis,'mudmass_03',[ii-1 0 0 0],[1 8 200 250]);
        
        Mdep=squeeze(sum(Mdep));
        Mdep(~mmask)=nan;
        
        dep = Mdep - Mdep1;
    %%%%%%%%%%%%%Wind
        u_wind=nc_varget(wind_file,'Uwind',[w_ind-1 0 0],[1 1 1]);
        v_wind=nc_varget(wind_file,'Vwind',[w_ind-1 0 0],[1 1 1]);
        
        w=(u_wind+sqrt(-1).*v_wind).*exp(sqrt(-1).*mangle(1,1));
    %%%%%%%%%%%River Discharge
        
        RD=nc_varget(river_file,'river_transport',[r_ind-1 0],[1 2]);
        RD=sum(RD);
%%%%%%%%%%%%%%%%%%%%%%%% sea surface elevation
        zeta = nc_varget(tmphis,'zeta',[ii-1 0 0],[1 200 250]);
        zeta(~mmask) = nan;
        
%%%%%%%%%%%%%%%%%%%%PLOTTING
        clf
        set(gcf,'pos',[20 100 nxx nyy],'color','k');
%%%%%%%%%%%%%%%%%%%%%%%%Adding the background, VIMS, and NSF logos
         a8=axes('pos',ax8);
         text(.05,.9,'June-July 2008','color','w');
         text(.05,.1,{'David Fugate',today},'color','w');
         text(0.2,.1,'City of Sanibel and Lee Co','color','w')
         text(0.05,0.25,datestr(gregorian(t(ii))))
         axis off
%         imshow(bkgrn);
%         axis off
%         
        axvl=axes('pos',vl_ax);
        imshow(vlogo);
        axis off
%                 
%         axnsfl=axes('pos',nsfl_ax);
%         imshow(nsflogo);
%         axis off
% 
%         axonrl=axes('pos',onrl_ax);
%         imshow(onrlogo);
%         axis off

%%%%%%%%%%%%%%%%%%%%%%%Plotting the 4 main panels        
        a1=axes('pos',ax1);
        pslice2(mrlon,mrlat,-Ms,[-35 -10]);
                hold on;
        %fillseg(coast);
        plot(coast(:,1),coast(:,2),'k-','linewidth',1.5)
        axis(Maxis);
        set(gca,'XtickLabel',[],'YtickLabel',[]);
        title({'Salinity';'Surface Currents'});
        pp=psliceuv(mrlon,mrlat,mw,8,.04,'r');
        set(pp,'linewidth',1.3);
        pp=psliceuv(177.92,-38.74,0+i*.3,1,.04,'w');
        text(177.905,-38.75,'0.3');
        text(177.905,-38.76,'m s^-^1');
        set(pp,'linewidth',1.3);
%        T2=sprintf('Hour %d',T);
%        text(177.91,-38.645,T2);
        hpc=pclegen2([-35:.05:-10],[0.04 .76 .2 .07]);
        hlb=get(hpc,'xticklabel');
        set(hpc,'xticklabel',-1*str2num(hlb),'Xcolor','w');
        contour(mrlon,mrlat,mh,[5:5:60],'k');
        dasp(26) %dcf cant find
        axis off
        
        a2=axes('pos',ax2);
        pslice2(mrlon,mrlat,log10(MmDI+.000001),[-4 0]);
        hold on;
        %fillseg(coast);
        plot(coast(:,1),coast(:,2),'k-','linewidth',1.5)
        axis(Maxis);
        set(gca,'XtickLabel',[],'YtickLabel',[]);
        title({'Depth Int. Sediment, log_1_0 kg m^-^2';'DI Sediment Transport, kg m^-^1s^-^1'});
        hold on
        pp=psliceuv(mrlon,mrlat,WT,6,4,'r');
        set(pp,'linewidth',2);
        hpc=pclegen2([-4:.04:0],[-0.02 .76 .2 .07]);
        set(hpc,'Xcolor','w')
        text(177.918,-38.698,'log_1_0');
        pp=psliceuv(177.92,-38.74,0+i*.5,1,.04,'w');
        text(177.905,-38.75,'0.5');
        text(177.905,-38.76,'kg m^-^1s^-^1');
        set(pp,'linewidth',1.3);
        contour(mrlon,mrlat,mh,[5:5:60],'k');
        [C1,C2]=contour(mrlon,mrlat,log10(MmDI+.000001),[0 0]);
        set(C2,'linewidth',2,'color',[.3 .3 .3]);
        [C1,C2]=contour(mrlon,mrlat,log10(MmDI+.000001),[1 1]);
        set(C2,'linewidth',2,'color',[.3 .3 .3]);
        [C1,C2]=contour(mrlon,mrlat,log10(MmDI+.000001),[2 2]);
        set(C2,'linewidth',2,'color',[.3 .3 .3]);
        dasp(26)
        axis off
%         if(rect(1)==rect(2))
%             plot(mrlon(rect(1),rect(3):rect(4)),mrlat(rect(1),rect(3):rect(4)),'k--','linewidth',2);
%         else
%             plot(mrlon(rect(1):rect(2),rect(3)),mrlat(rect(1):rect(2),rect(3)),'k--','linewidth',2);
%         end
        
        a3=axes('pos',ax3);
        pslice2(mrlon,mrlat,dep,[-.2 .2]);
        hold on;
        %fillseg(coast);
        plot(coast(:,1),coast(:,2),'k-','linewidth',1.5);
        axis(Maxis);
        set(gca,'XtickLabel',[],'YtickLabel',[]);
        title({'Deposition, kg m^-^2';'Bottom Currents'});
        hpc=pclegen2([-.2:.04:.2],[-0.02 .76 .2 .07]);
        set(hpc,'Xcolor','w');
        hold on
        contour(mrlon,mrlat,mh,[5:5:60],'k');
        [C1,C2]=contour(mrlon,mrlat,Mdep,[0 0]);
        set(C2,'linewidth',2,'color',[.3 .3 .3]);
        [C1,C2]=contour(mrlon,mrlat,Mdep,[.1 .1]);
        set(C2,'linewidth',2,'color',[.3 .3 .3]);
        [C1,C2]=contour(mrlon,mrlat,Mdep,[-.1 -.1]);
        set(C2,'linewidth',2,'color',[.3 .3 .3]);
        text(177.918,-38.698,'log_1_0');
        pp=psliceuv(mrlon,mrlat,mwb,8,.2,'r');
        set(pp,'linewidth',1.3);
        pp=psliceuv(177.92,-38.74,0+i*.1,1,.2,'w');
        text(177.905,-38.75,'0.1');
        text(177.905,-38.76,'m s^-^1');
        set(pp,'linewidth',1.3);    
        dasp(26)
        axis off
        
        a4=axes('pos',ax4);
        pslice2(mrlon,mrlat,zeta,[-.5 .5]);
        hold on
        %fillseg(coast);
        plot(coast(:,1),coast(:,2),'k-','linewidth',1.5)
        axis(Maxis);
        set(gca,'XtickLabel',[],'YtickLabel',[]);
        title('Sea Surface Height, m');
        hpc=pclegen2([-.5:.01:.5],[-0.02 .76 .2 .07]);
        set(hpc,'Xcolor','w')
        contour(mrlon,mrlat,mh,[5:5:60],'k');
        dasp(26)
        tmp = gregorian(t(ii));
        ts = datestr(tmp,0);
        axis off

%%%%%%%%%%%%%%%Plotting wind, and river forcing
        a5=axes('pos',ax5);
        p=plot(20*cos(theta),20*sin(theta));
        set(p,'linewidth',2,'color','c');
        axis off;
        dasp;
        hold on;
        plot(0,0,'bo','Markersize',12,'markerfacecolor','b');
        H=quiver(0,0,real(w),imag(w));
        set(H,'linewidth',3,'color','c');
        axis([-20 20 -20 20]);
        mag=sqrt(real(w).^2+imag(w).^2);
        mag=sprintf('%2.1f m s^-^1',mag);
        text(21.5,6,'Wind');
        text(21.5,2,'Spd/Dir');
        text(21.5,-2,mag);
        
%         a6=axes('pos',ax6);
%         P=plot(rn,RD_F,'c','linewidth',2);
%         axis tight;
%         %set(p,'linewidth'
%         hold on;
%         plot(rn(r_ind),RD,'r*','markersize',20);
%         xlabel('Date');
%         text(rn(1)+5,1000,'River Discharge, m^3 s^-^1');
%         datetick('x',12,'keeplimits');
%         set(gca,'Ytick',[400 800 1200],'Xcolor','w','Ycolor','w','color','none')
%         grid on
        
%%%%%%%%%%%%%%%Plotting SSH        
%         a7=axes('pos',ax7);
%         fill(pm,zeta,[.5 .2 .1]);
%         hold on
%         plot(pm(1:50),zeta(1:50),'b','linewidth',2);
%         axis([0 12 -1 1]);
%         ax=axis;
%         xtmp = [ax(1) ax(2)];
%         ytmp = [0 0];
%         plot(xtmp,ytmp,'w--','linewidth',2.5);
%         plot(ztmp,z2,'g-','linewidth',2);
%         xlabel({'Distance, km','Time, hrs'});
%         ylabel('Meters');
%         title('Sea Surface Height');
%         set(gca,'Xcolor','w','Ycolor','w','color','none');
%        axis off
%         a7=axes('pos',ax7);
%         fill(pm,zeta,[.5 .2 .1]);
%         hold on
%         plot(pm(1:50),zeta(1:50),'b','linewidth',2);
%         axis([0 12 -.020 .020]);
%         set(gca,'Ytick',[-.020 0 0.020],'YtickLabel',[-20 0 20],'YAxisLocation','right');
%         ax=axis;
%         xtmp = [ax(1) ax(2)];
%         ytmp = [0 0];
%         plot(xtmp,ytmp,'w--','linewidth',2.5);
% %         plot(ztmp,z2,'g-','linewidth',2);
%         xlabel({'Distance, km'});%,'Time, hrs'});
%         ylabel('mm');
%         title('Sea Surface Height');
%         set(gca,'Xcolor','w','Ycolor','w','color','none');
%%%%%%%%%%Plotting vertical profile of currents
%         Ltmp = length(Z1);
%         a9=axes('pos',ax9);
%         Q=quiver3(XX,YY,Z1/1,u,v,W(1:Ltmp));
%         set(Q,'linewidth',2,'color','c');
%         hold on
%         plot3(XX,YY,Z1/1,'w-');
%         view(10,40);
%         %axis equal
%         axis([-2 2 -2 2 min(Z1/1) 0]);
%         set(gca,'dataaspectratio',[1.2 1.2 3]);
%         axis off;
%         
%         a10=axes('pos',ax10);
%         Q=quiver3(XX2,YY2,Z2/1,u2,v2,W2(1:Ltmp));
%         set(Q,'linewidth',2,'color','c');
%         hold on
%         plot3(XX2,YY2,Z2/1,'w-');
%         view(10,40);
%         %axis equal
%         axis([-2 2 -2 2 min(Z2/1) 0]);
%         set(gca,'dataaspectratio',[1.2 1.2 3]);
%         axis off;
% 
%         a11=axes('pos',ax11);
%         Q=quiver3(XX3,YY3,Z3/1,u3,v3,W3(1:Ltmp));
%         set(Q,'linewidth',2,'color','c');
%         hold on
%         plot3(XX3,YY3,Z3/1,'w-');
%         view(10,40);
%         %axis equal
%         axis([-2 2 -2 2 min(Z3/1) 0]);
%         set(gca,'dataaspectratio',[1.2 1.2 3]);
%         axis off;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Plotting the sediment transport out of the bay
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%mouth
%         a12 = axes('pos',ax12);
%         tmpr = range(squeeze(ST(k,:,:)));
%         if(tmpr(3)>1)
%             tmpca = [-.5 .5];
%             axlabel = {'-.5','-.25','0','.25','.5'};
%         elseif(tmpr(3)>.1)
%             tmpca = [-.25 .25];
%             axlabel = {'-.25','-.125','0','.125','.25'};
%         elseif(tmpr(3)>.01)
%             tmpca = [-.01 .01];
%             axlabel = {'-.01','-.005','0','.005','.01'};
%         else
%             tmpca = [-.001 .001];
%             axlabel = {'-.001','-.0005','0','.0005','.001'};
%         end
%         %[P,PP] = contour(D2,squeeze(ZS(k,:,:)),squeeze(ST(k,:,:)),tmpca,'k');
%         pcolor(D2,squeeze(ZS(k,:,:)),squeeze(ST(k,:,:)));
%         caxis(tmpca);
%         shading flat
%         hold on
% %        set(PP,'linewidth',1.5');
%         P = colorbar;
%         set(P,'Yticklabel',axlabel);
%         [P,PP] = contour(D2,squeeze(ZS(k,:,:)),squeeze(ST(k,:,:)),[0 0],'k');
%         set(PP,'linewidth',2');
%         fill(D,bed_elev,[.5 .2 .1])
%         P = plot(D,bed_elev,'k','linewidth',3);
%         axis([min(D) max(D) -30 1]);
%         set(gca,'Xcolor','w','Ycolor','w','color','none','Xtick',[2000 4000 6000 8000 10000],...
%             'Xticklabel',[2 4 6 8 10]);
%         ylabel('Depth, m')
%         xlabel('Distance, km')
%         title('Sediment Transport (kg m^-^2 s^-^1)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Setting the figure size and saving
        %set(gcf,'pos',[20 100 nxx nyy],'color','k');
        outfile=sprintf('swf_wet.%4d.png',k);

f = getframe(gcf);
        [X, map] = frame2im(f);
        imwrite(X,outfile,'png');
        
        clear u v w Mdep MmDI Ms Mm begin
    end
end
% profiled on Poverty at 892s, not local = 1900


      
        
        
        
        
        
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
        
        
        
        
        
        
        
