% puffsb.m
%
% puff.m adapted for use during Sea Breeze cruises (SRV, May 24, 2004)
%
% Models the distribution of tracer from a source, 
% given a time series of advection and constant diffusivities.
% 
% SRV, 3/13/01; revised 10/31/01, 12/02/01, 12/04/01, 1/21/02
% 2/19/02 added ability to create c-t-series at particular point(s)
% 8/9/02 saved backup in preparation for adding ability to create c-t-series at moving points
% 8/26/02 saved backup in prep for adding Mojfeld-like current time series summed with range of means

clear;
close all;

mofit=0;  %If changed to 1, then uncomment mofmean for and end statements below 
%meanvres=1; for mofmean=[0:meanvres:12],
%for mofmean=[0:5],
if mofit==1,
  mofmean
else,
  mofmean=[]
end;


%% Indicate desired routines (plots and calculations) with 1's
puffplot=0;    % plots each frame on screen, rather than just calculating matrices and related statistics...
printframes=0; % prints each frame as a .eps file 
oblabel=0;     % checks to see if FlowMow observation was made during current frame; overlays label is so...
fluxcalc=1;    % gets flux time series and total flux through control volume walls (1 depth only)
fluxplot=0;    % Plot instantaneous net flux and fluxes through particular walls
vdthplot=0;    % Plot "bow tie" plot (observed v versus modeled deltaTheta) 
displace=0;    % Plot integrated displacement of a particle througout a run (the larval perspective)
endfluxplot=0; % Plot a summary of instantaneous net flux and fluxes through particular walls (including means)
pt_c=0;        % option of tracking concentration at a particular point throughout a run (specified below)

%% Optionally, indicate a point where you want to track concentration throughout a run
if pt_c==1,
  if 0,
  % Based on GPS coords given in Thomson et al May 16,2002 abstract; geoconverted 
  xpt_utm=493104.11;
  ypt_utm=5311711.82;
  pt_name=['Wasp2161']   % position of WASP current meter at 2161m in summer 2000
  end;
  if 0,
  xpt_utm=493031.3;
  ypt_utm=5310766.1;
  pt_name=['MZ_cm_95']   % Mixing Zephyr current meter in 1995 (NE corner MEF)
  end;
  if 1,
  xpt_utm=492995;
  ypt_utm=5311225;
  pt_name=['NoMEFStn']   % A point just north of MEF, near FlowMow station 28
  end;
end;

%% Initialize some variables... 
southest=[];
northest=[];
westest=[];
eastest=[];
near=[];

% gather times/locations of Flow Mow observations 
time_table
load /snap/flowmow/ABE/Metadata/abe_mjd_loc.mat;
% initialize counters (for main iterative loop)
ll=1; mm=1; nn=1; rmflag=0; 
presentabe=0;
abeobstart=abe_mjd_loc(:,7);
abeobend=abe_mjd_loc(:,8);
load /snap/flowmow/ctd/ctdstats.txt;
ctdstart=ctdstats(:,1);
ctdend=ctdstats(:,2);
presentstn=0;
ctdlist=[1:40];
%nearmeflist=[3 4 5 6 7 8 9 10 11 12 13 14 15 20 21 22 24 25 28 30 31 32 33 35 36 37 38 39 40];
% translate and rotate observation locations into the model space (020True)
%rot=20; % degrees, note overall trend is 20~030, but this is local to MEF
load /snap/flowmow/locations/fmbox.xy
load /snap/flowmow/locations/fmbox_center.xy  % this is taken as the puff origin
rot=90-(180/pi)*atan2(fmbox(2,2)-fmbox(3,2),fmbox(2,1)-fmbox(3,1));  % +17.35 degrees

% find ctd box centers (warning: need better algorithm to find centers of arbitrary parallelograms
bbcenterx=ctdstats(:,3)+(ctdstats(:,5)-ctdstats(:,3))/2;  % bounding box centerx
bbcentery=ctdstats(:,4)+(ctdstats(:,6)-ctdstats(:,4))/2;  % bounding box centery

%translate so each ctd bounding box is centered on itself
ctdloc(:,1)=ctdstats(:,3)-bbcenterx;  %sw x
ctdloc(:,2)=ctdstats(:,4)-bbcentery;  %sw y
ctdloc(:,3)=ctdstats(:,3)-bbcenterx;  %nw x
ctdloc(:,4)=ctdstats(:,6)-bbcentery;  %nw y
ctdloc(:,5)=ctdstats(:,5)-bbcenterx;  %ne x
ctdloc(:,6)=ctdstats(:,6)-bbcentery;  %ne y
ctdloc(:,7)=ctdstats(:,5)-bbcenterx;  %se x
ctdloc(:,8)=ctdstats(:,4)-bbcentery;  %se y

% code to rotate each point of each box
for pp=1:4,
  newangle=atan2(ctdloc(:,pp*2),ctdloc(:,(pp*2)-1))+(rot*pi/180);
  newmag=sqrt(ctdloc(:,(pp*2)-1).^2+ctdloc(:,pp*2).^2);
  ctdloc(:,(pp*2)-1)= cos(newangle).*newmag;
  ctdloc(:,pp*2)= sin(newangle).*newmag;
end;

%translate to puff model origin 
offsetx=bbcenterx-fmbox_center(1);
offsety=bbcentery-fmbox_center(2);
newoffsetdir=atan2(offsety,offsetx)+(rot*pi/180);
offmag=sqrt(offsetx.^2+offsety.^2);
newoffsetx=cos(newoffsetdir).*offmag;
newoffsety=sin(newoffsetdir).*offmag;
ctdloc(:,1)=ctdloc(:,1)+newoffsetx;  %sw x 
ctdloc(:,2)=ctdloc(:,2)+newoffsety;  %sw y
ctdloc(:,3)=ctdloc(:,3)+newoffsetx;  %nw x 
ctdloc(:,4)=ctdloc(:,4)+newoffsety;  %nw y
ctdloc(:,5)=ctdloc(:,5)+newoffsetx;  %ne x 
ctdloc(:,6)=ctdloc(:,6)+newoffsety;  %ne y
ctdloc(:,7)=ctdloc(:,7)+newoffsetx;  %se x 
ctdloc(:,8)=ctdloc(:,8)+newoffsety;  %se y
ctdloc=ctdloc/1000; % convert from meters to km

if pt_c==1,
  %rotate the utm location of a specific monitoring point into the model space
  xptmag=(xpt_utm-fmbox_center(1));
  yptmag=(ypt_utm-fmbox_center(2));
  ptmag=sqrt(xptmag^2+yptmag^2);
  ptangle=atan2(yptmag,xptmag)+(rot*pi/180); % =69.5+17.35=97o=almost directly up valley
  xpt=cos(ptangle)*ptmag/1000  % location of specific monitoring point in KM!
  ypt=sin(ptangle)*ptmag/1000
end;

if mofit==1,
  sim='mofl';
else,
  sim=input('Over what period would you like to run the puff model? [wasp,fmcm,fmcs,mz_1,mz12,mzcm, mofj, or rand]  ','s')
end;
if sim=='wasp',
 simstart=51685;	  % wasp2161 start
 simend=51826.7708333;    % wasp2161 end
elseif sim=='fmcm',
 simstart=51752.8333333;  % fm_cm start
 simend=51826.7916666667; % fm_cm end
elseif sim=='fmcs',
 simstart=51760;  % fm_cruise start
 %simstart=51763.8 % right before first abe ob
 simend=51777;    % fm_cruise end
elseif sim=='mz_1',
 % Atlantis II/Voyage 132/Leg 08 
 simstart=49866;  % First Tow:      May29   Jday149
 simend=49884;    % Last Tow:       June16  Jday167
elseif sim=='mz12',
 simstart=49862;  % mz_cruise start (both legs)
 simend=49988;    % mz_cruise end   (both legs)
elseif sim=='mzcm', 
 simstart=49866.65; % mz_cm start
 simend=49918.6;    % mz_cm end
elseif sim=='rand',
 simstart=51764.2;  % random start
 simend=51774;      % random end
elseif sim=='mofj',         %% Short duration run
  simstart=51757;           % Mofjeld prediction start 
  %simend=  51759;           % end of a few semi/durnal oscillations
  simend=  51761;           % end of ~6 semi/durnal oscillations 
elseif sim=='mofl',        %% long duration run
  simstart=51757;           % Mofjeld prediction start
  simend= 51772;           % approximate end of 1st spring/neap cycle (25 days) (~25min/<v> to run)
  %simend=  51816.979166;    % Mofjeld prediction end
end;

if mofit==1,
  currsrc='mofjd'
else,
  currsrc=input('What current record would you like to use? [synth,mofjd,wasp0,cm050,cm150,cm250,mz025,mz300...]  ','s')
end;

%% Define boundaries of model space and source location and properties
%emax=50; wmax=-150; nmax=100; smax=-100;  % km!!
%emax=5; wmax=-5; nmax=5; smax=-5;
emax=1.0; wmax=-1.0; nmax=1.0; smax=-1.0;  %1 km on a side box
near_multiplier=2;  % multiplied *emax to yield range at which puffs die
maxage=50;          % age in hours that determines if puffs die 
                    % (useful when mean flow is zero) Should be altered when diffusivity is changed
if currsrc=='mofjd',
  hres=0.05;   % 50m
  eside=+0.15; wside=-0.15; nside=0.35; sside=-0.35;  % use when generating puff animations 
  %eside=1; wside=-1; nside=0.4; sside=-0.4;  % use when getting accurate (entire) N/S flux
else,
  hres=0.1;     % 100m
  eside=+0.2; wside=-0.2; nside=0.4; sside=-0.4;  % MEF flux perimeter
  %eside=+0.15; wside=-0.15; nside=0.35; sside=-0.35;  % more exact MEF flux perimeter
end;
% #s below won't work unless you change the indexing scheme (much further below)
%eside=.15; wside=-.15; nside=.35; sside=-.35;  % MEF flux perimeter [km]
%eside=2*hres; wside=-2*hres; nside=4*hres; sside=-4*hres;  % MEF flux perimeter
zdeep=-2.300; zshoal=-1.700;  %km
vres=0.1;   %  100m
if (currsrc=='mz300'),
  tres=0.25   %  HOURS -- decrease this time step to smooth animation in high 
	      %  speed currents (in which displacement is bigger than init puff)
else,
  tres=0.5
end;

%cp=4.2e3;    % J/kg/oC
cp=3816;
%rho=1000;    % kg/m3
rho=1030;

xsrc=0; ysrc=0;    % source location
agesrc=tres;
dx=100;
dy=dx;	           % full width of initial top-hat puffs in METERS 
		   % (divided by 2 in gx,gy eqn so, e.g. the half width of initial top hat is 25m) 
dz=vres*1000;      % full height should equal vres (m) if fluxes are to be compared
%A=dx*dy;
%wi=0.1;
%Q=A*wi;             % = 10^4*10^-2= 10^2 m3/s, smaller would be 10^2*10^-1 = 10 m3/s 
%qi=Q*dTi*rho*cp;    % = 10^2*10*10^3*3.8*10^3 = 3.8*10^9 W = 3800 MW 
                     % or for smaller area... = 38 MW 
if (currsrc=='mz300' | currsrc=='cm250' | currsrc=='cm200'),
  qi=600*10^6
else,
  qi=300*10^6          % intermediate estimate of MEF diffuse source heat flux in W
end;
ci=qi*tres*3600/(dx*dy*dz*rho*cp)  %initial theta anom in oC, all units mks

xp=[xsrc];             % initialize puff locations and ages
yp=[ysrc];             % km
agep=[agesrc];
zp=300;

%kx=1.6;
kx=0.4;             % dispersion characteristics (diffusivity from Okubo in m2/s)
%kx=0.02;             % value used in Wetzler98
ky=kx;               % Consider changing maxage if you alter these... 
kz=0.1;

xi=[wmax:hres:emax];   % grid set up, all units are km
yi=[smax:hres:nmax];
%zi=[zdeep:vres:zshoal];
[x,y]=meshgrid(xi,yi);
z=zp*ones(size(x));
c=zeros(size(x));

%% get the subarrays for the sides of the flux box
% note that the box must contain the origin for this to work...
origin_xind=(-wmax/hres)+1;  % get the origin indices for the c matrix
origin_yind=(-smax/hres)+1;  % 11,11 if hres=0.1km and w/smax = 1km
eind=origin_xind+eside/emax/hres; 
wind=origin_xind-wside/wmax/hres;
nind=origin_yind+nside/nmax/hres; 
sind=origin_yind-sside/smax/hres;
nssub=[sind:1:nind];
ewsub=[wind:1:eind];
ceind=sub2ind(size(c),nssub,ones(size(nssub))*eind);
cwind=sub2ind(size(c),nssub,ones(size(nssub))*wind);
cnind=sub2ind(size(c),ones(size(ewsub))*nind,ewsub);
csind=sub2ind(size(c),ones(size(ewsub))*sind,ewsub);

%%% LOAD CURRENT METER RECORD
%%% Ensure that velocities are in m/s units!!
if ((currsrc=='wasp0')|(currsrc=='cm050')|(currsrc=='cm100')|(currsrc=='cm150')|(currsrc=='cm200')|(currsrc=='cm250')|(currsrc=='mz025')|(currsrc=='mz300')),  %% use field current data
 % Load the current meter record(s)
 data='real';
 if (currsrc=='wasp0'),  % use wasp data
  datasrc='WASP'
  load /snap/flowmow/currents/wasp/no_filter.hourly/wasp2161.mat
  fmwaspind=find(wasp2161.mjd>=simstart & wasp2161.mjd<=simend);
  u=wasp2161.u(fmwaspind); v=wasp2161.v(fmwaspind);
  meanu=mean(u)
  meanv=mean(v)
  currmjd=wasp2161.mjd(fmwaspind);
  curmjdstart=currmjd(1);
  cmres=24*(currmjd(2)-curmjdstart);  % sampling interval of curr meter in hrs
  curmjdend=currmjd(length(currmjd));
 elseif ((currsrc=='cm050')|(currsrc=='cm100')|(currsrc=='cm150')|(currsrc=='cm200')|(currsrc=='cm250')) 
  % use Flow Mow 5-meter mooring data
  load /snap/flowmow/currents/mooring/no_filter.hourly/all.mat;
  if (currsrc=='cm050'), % use deepest meter
   datasrc='CM50'
   fmmab50ind=find(mab50.mjd>=simstart & mab50.mjd<=simend);
   u=mab50.u(fmmab50ind); v=mab50.v(fmmab50ind);  
   meanu=mean(u)
   meanv=mean(v)
   currmjd=mab50.mjd(fmmab50ind);
   curmjdstart=mab50.mjd(fmmab50ind(1));
   curmjdend=mab50.mjd(fmmab50ind(length(fmmab50ind)));
  elseif (currsrc=='cm100'),
   datasrc='CM100'
   fmmab100ind=find(mab100.mjd>=simstart & mab100.mjd<=simend);
   u=mab100.u(fmmab100ind); v=mab100.v(fmmab100ind);
   meanu=mean(u)
   meanv=mean(v)
   currmjd=mab100.mjd(fmmab100ind);
   curmjdstart=mab100.mjd(fmmab100ind(1));
   curmjdend=mab100.mjd(fmmab100ind(length(fmmab100ind)));
  elseif (currsrc=='cm150'),
   datasrc='CM150'
   fmmab150ind=find(mab150.mjd>=simstart & mab150.mjd<=simend);
   u=mab150.u(fmmab150ind); v=mab150.v(fmmab150ind);
   meanu=mean(u)
   meanv=mean(v)
   currmjd=mab150.mjd(fmmab150ind);
   curmjdstart=mab150.mjd(fmmab150ind(1));
   curmjdend=mab150.mjd(fmmab150ind(length(fmmab150ind)));
  elseif (currsrc=='cm200'), 
   datasrc='CM200'
   fmmab200ind=find(mab200.mjd>=simstart & mab200.mjd<=simend);
   u=mab200.u(fmmab200ind); v=mab200.v(fmmab200ind);
   meanu=mean(u)
   meanv=mean(v)
   currmjd=mab200.mjd(fmmab200ind);
   curmjdstart=mab200.mjd(fmmab200ind(1));
   curmjdend=mab200.mjd(fmmab200ind(length(fmmab200ind)));
  elseif (currsrc=='cm250'), % use shallowest meter
   datasrc='CM250'
   fmmab250ind=find(mab250.mjd>=simstart & mab250.mjd<=simend);
   u=mab250.u(fmmab250ind); v=mab250.v(fmmab250ind);
   meanu=mean(u)
   meanv=mean(v)
   currmjd=mab250.mjd(fmmab250ind);
   curmjdstart=mab250.mjd(fmmab250ind(1));
   curmjdend=mab250.mjd(fmmab250ind(length(fmmab250ind)));
  end;
  cmres=24*(currmjd(2)-curmjdstart);  % sampling interval of curr meter in hrs
 elseif (currsrc=='mz025'), % use near-bottom data from 1995
  datasrc='mz025'
  load /shuksan/zephyr/currents/cm1995all.mat
  mzind=find(cm2175.mjd>=simstart & cm2175.mjd<=simend);
  u=cm2175.spd_mps(mzind).*cos(pi/180*(90-cm2175.dirT(mzind)));
  v=cm2175.spd_mps(mzind).*sin(pi/180*(90-cm2175.dirT(mzind)));
  meanu=mean(u)
  meanv=mean(v)
  currmjd=cm2175.mjd(mzind);
  curmjdstart=currmjd(1);
  cmres=24*(currmjd(2)-curmjdstart);  % sampling interval of curr meter in hrs
  curmjdend=currmjd(length(currmjd));
 elseif (currsrc=='mz300'), % use upper level data from 1995
  datasrc='mz300'
  load /shuksan/zephyr/currents/cm1995all.mat
  mzind=find(cm1900.mjd>=simstart & cm1900.mjd<=simend);
  u=cm1900.spd_mps(mzind).*cos(pi/180*(90-cm1900.dirT(mzind)));
  v=cm1900.spd_mps(mzind).*sin(pi/180*(90-cm1900.dirT(mzind)));
  meanu=mean(u)
  meanv=mean(v)
  currmjd=cm1900.mjd(mzind);
  curmjdstart=currmjd(1);
  cmres=24*(currmjd(2)-curmjdstart);  % sampling interval of curr meter in hrs
  curmjdend=currmjd(length(currmjd));
 end;
 tmax=24*(curmjdend-curmjdstart); 
 %% convert: rotate the u,v (trueN) components into the model space (axis perpendicular/parallel)
 %% where axis trend is defined by variable "rot"
 phi=atan2(v,u);
 mag=sqrt(u.^2+v.^2);
 uroti=cos(phi+rot*pi/180).*mag;
 vroti=sin(phi+rot*pi/180).*mag;
elseif (currsrc=='mofjd'),  %% use Hal Mojfeld's harmonic analysis time series
 data='real';
 datasrc=['Harmonic+' int2str(mofmean)]
 load /snap/flowmow/currents/predicted/pred.mat; 
 mofind=find(pred.mjd>=simstart & pred.mjd<=simend);
 maxpredv=max(pred.v(mofind)/100)
 minpredv=min(pred.v(mofind)/100)
 meanpredv=mean(pred.v(mofind)/100)
 %mofmean=input('What mean would you like to add to the Mofjeld time series? [0, 2, 4, etc; in cm/s!]  ','s')
 %mofmean=eval(mofmean);
 vroti=(mofmean+pred.v(mofind))/100;      % note pred.v is in cm/s, so convert here to m/s
 uroti=zeros(size(vroti));  % assume the predicted v is along axis...
 meanv=mean(vroti)
 currmjd=pred.mjd(mofind);
 curmjdstart=currmjd(1);
 cmres=24*(currmjd(2)-curmjdstart);  
 curmjdend=currmjd(length(currmjd));
 tmax=24*(curmjdend-curmjdstart);
else,  %% use synthetic current data
 data='synth';
 datasrc='Synthetic'
 umax  =0.0;     % remember: these u,v components will be rotated to be perpendicular,parallel to axis
 vmax  =0.025;   % m/s
 udrift=0.00; % specify in m/s
 vdrift=0.025;  
 period=12;
 %period=24;
 pi=3.14159;
 omega=2*pi/period;
 tmax=48;  %hr
end;

if (data=='real'),
 %% decimate urot,vrot arrays to the specified frequency
 %urot=interp1q([0:cmres:tmax]',uroti,[0:tres:tmax]');
 %vrot=interp1q([0:cmres:tmax]',vroti,[0:tres:tmax]'); 
 urot=interp1q(currmjd,uroti,[curmjdstart:tres/24:curmjdend]');
 vrot=interp1q(currmjd,vroti,[curmjdstart:tres/24:curmjdend]');
 % careful, the following 2 lines can cause NaNs and thereafter c-fields with Nans!
 %urot=interp1q(currmjd,uroti,[simstart:tres/24:simend]');
 %vrot=interp1q(currmjd,vroti,[simstart:tres/24:simend]');
end;

%%% Iterate from t=0 to tmax, calculating, surfacing, and displaying concentration (c) 
%% at each time step and waiting for user ok to go to next time step

if puffplot==1,
  figure(1)
  clf
  set(gcf,'menubar','none');
end;
if fluxplot==1,
 figure(2)
 clf
 %set(gcf,'Position',[701   522   560   420]);
end;
if vdthplot==1,
 figure(3)
 clf;
 %set(gcf,'Position',[702    13   560   420]);
end;

t_elapsed=[]; 
netflux=[]; n_totflux=[]; s_totflux=[]; e_totflux=[]; w_totflux=[]; vseries=[]; meancseries=[];    % initialize arrays
c_nseries=[]; c_sseries=[]; c_wseries=[]; c_eseries=[];
c_ptseries=[];
npuffs=1;
j=1; % counts number of time steps (circuits of main loop), increments current meter tseries
k=1; % used to count number of time steps and pause, if pause option is selected
for t=0:tres:tmax,
 %reset the concentration field and sum all puffs' contributions
 c=zeros(size(x));
 cpt=0;
 dt=tres*3600;
 for p=1:npuffs,
    % get concentration contribution of a puff #n for all grid points
    if (puffplot==1 | fluxplot==1 | fluxcalc==1),
      % top-hat --> gaussian profile function is:
      gx=erf((.5*dx+(x-xp(p))*1000)/sqrt(4*kx*agep(p)*dt))+erf((.5*dx-(x-xp(p))*1000)/sqrt(4*kx*agep(p)*dt));
      gy=erf((.5*dy+(y-yp(p))*1000)/sqrt(4*ky*agep(p)*dt))+erf((.5*dy-(y-yp(p))*1000)/sqrt(4*ky*agep(p)*dt));
      % hold vertical dispersion to ~zero
      %gz=erf( (0.5*dz+(zp-z)*1000)/sqrt(4*kz*.001) )+erf( (0.5*dz+(zp+z)*1000)/sqrt(4*kz*0.001) );
      % add it to the extant gridded concentration field
      c=c+ci.*(0.5*gx).*(0.5*gy);
    end;

    % Probably should put in a new if here for fluxcalc==1
    % This would get c at each point on a specified control surface at each time step
    % rather than dictating that the surface be exactly on the x,y grid, and of the 
    % same spatial resolution...
    % 	THIS SOLVES THE HELL OF GETTING SUBINDICES FOR ARBITRARILY PLACED SURFACES

    % now calculate concentration at a chosen point in the model region
    if pt_c==1,
      xdist=(xpt-xp(p))*1000;  % distance from each puff to the chosen point in METERS
      ydist=(ypt-yp(p))*1000; 
      gpx=erf((.5*dx+xdist)/sqrt(4*kx*agep(p)*dt))+erf((.5*dx-xdist)/sqrt(4*kx*agep(p)*dt));
      gpy=erf((.5*dy+ydist)/sqrt(4*ky*agep(p)*dt))+erf((.5*dy-ydist)/sqrt(4*ky*agep(p)*dt));
      cpt=cpt+ci.*(0.5*gpx).*(0.5*gpy);
    end;
 end;

 if (puffplot==1),
  figure(1)
  plotdate=curmjdstart+t/24; plotday=floor(plotdate); plothr=24*(plotdate-plotday);
  if plothr<9.99,
    datelabel=['Day: ' sprintf('%5.0f',plotday) ', hr: 0' sprintf('%2.1f',plothr)];
  else,
    datelabel=['Day: ' sprintf('%5.0f',plotday) ', hr: ' sprintf('%2.1f',plothr)];
  end;
  if ll==1, % first iteration
    % surface the concentration field and plot it...
    %colormap('hot');
    mapp=colormap('jet');
    mapp(1,:)=[1 1 1];                % make first row of colortable white (transparent)
    colormap(mapp);
    thmin=0; thinc=0.01; thmax=0.1; 
    thresh=0.01;                      % c<thresh is transparent -- defines plume edge
    thetas=[thmin:thinc:thmax];
    %img=imagesc(xi,yi,c);
    %set(gca,'YDir','normal')
    spp=subplot('position',[.1 .1 .75 .55]);
    img=imagesc(xi,yi,c);
    set(gca,'YDir','normal')
    %cf=contourf(xi,yi,c,thetas);
    %sfc=surface(xi,yi,c);
    axis equal;
    hold on;
    grid on;
    xlabel('X, km','FontSize',12);
    ylabel('Y, km','FontSize',12);
    %zlabel('Z, km');
    %   plot mef perimeter, puff locations, and source
    plot([wside eside eside wside wside],[nside nside sside sside nside],'k-.','LineWidth',1);
    %plot(xp,yp,'b.');
    plot(xsrc,ysrc,'r*');
    %datetx=text(wmax+(emax-wmax)/20,smax-(smax-nmax)*2/20,datelabel,'FontSize',14);
    title(datelabel,'FontSize',12,'FontWeight','bold','HorizontalAlignment','right');
    caxis([thmin thmax]);
    cb=colorbar('vert');
    set(cb,'YTick',thetas);
    set(cb,'YLim',[thresh thmax]);
    axis([wmax emax smax nmax]);
    set(gca,'YTick',[smax:(nmax-smax)/10:nmax]');
    set(gca,'XTick',[wmax:(emax-wmax)/10:emax]');
    if pt_c==1,
      % plot the location of a specified observation point
      plot(xpt,ypt,'kh')
    end;
 
    spt=subplot('position',[.15 .75 .4 .15]);
    %ptit=['Puff model: H_i=',int2str(qi/10^6),'MW, dx=',int2str(dx),'m, k=',sprintf('%2.1f',kx),')'];
    ptit=['Puff model: Net flux and \Delta\theta distribution'];
    tplot=plot(0,0,'b');
    hold on
    grid on
    %fax=axis;
    if (currsrc=='mz300' | currsrc=='cm250' | currsrc=='cm200'),
      fmax=5000; fmin=-fmax;  % MW
    else,
      fmax=1000; fmin=-fmax;  % MW
    end;
    axis([0 tmax fmin fmax])
    title(ptit,'FontSize',14);
    ylabel('Net Flux, MW')
    xlabel('Elapsed time, hr')
    
    spc=subplot('position',[.7 .75 .15 .15]);
    axis([-.15 .15 -.15 .15]);
    datalabel=['Current src: ',sprintf('%s',datasrc)];
    %text(wmax+(emax-wmax)/20,smax-(smax-nmax)*1/20,datalabel,'FontSize',10);
    % use typical current speed 0.05m/s 
    truenorthx=0.05*cos((90+rot)*pi/180); truenorthy=0.05*sin((90+rot)*pi/180);
    quiver(0,0,truenorthx,truenorthy,'b');
    axis equal;
    if (currsrc=='mz300'),
      axis([-.50 .50 -.50 .50]);
    else,
      axis([-.1 .1 -.1 .1]);
    end;
    title(datalabel)
    xlabel('u, m/s')
    ylabel('v, m/s')
    hold on
    grid on
    % plot the urot,vrot vector over the true north vector...
    qv=quiver(0,0,urot(j),vrot(j),'k.');  % . means no arrow head
    %text(emax-(emax-wmax)/3,smax-(smax-nmax)/20,'5cm/s True North','Color',[0 0 1],'FontSize',10);

    % Re-align plots!
    set(spp,'Position',[0.05 0.1 0.6 0.6],'DataAspectRatio',[1 1 1])
    set(spt,'Position',[0.15 0.8 .75 0.1])
    set(spc,'Position',[0.7 0.5 .15 0.15])
    set(cb,'Position',[0.75 0.1 0.05 0.25])
    ll=0;
  else,  % just reset data, rather than calling plot functions again... 
    subplot(spp),set(img,'CData',c),title(datelabel,'FontSize',12,'FontWeight','bold','HorizontalAlignment','right');
    subplot(spt),title(ptit,'FontSize',14);
    subplot(spt),set(tplot,'XData',t_elapsed,'YData',[netflux/1e6]);
    subplot(spc),set(qv(1),'XData',[0 urot(j)],'YData',[0 vrot(j)]);  % reset mag of current vector shaft...
    % realign subplots
    %set(spp,'Position',[0.1 0.1 0.65 0.65],'DataAspectRatio',[1 1 1])
    set(spp,'Position',[0.05 0.1 0.6 0.6],'DataAspectRatio',[1 1 1])
    set(spt,'Position',[0.15 0.8 .75 0.1])
    set(spc,'Position',[0.7 0.5 .15 0.15])
    set(cb,'Position',[0.75 0.1 0.05 0.25])
  end;

  % plot white points way beyond model region and matlab plot region to 
  % standarize cropping margins on all sides of .eps figures 
  % -- this stops jumpiness in animation...
  %xdistant=[100 100 -100 -100]; ydistant=[100 -100 100 -100];  % all in km
  %plot(xdistant,ydistant,'w.')

 if oblabel==1,
  % Check to see if an observation (CTD or ABE) was made during this time step
  % If so, pause the model and plot the observation location
  if (nn==1),  % if simulation is in first t-step, find # of next or present ctd ob
    previousstni=find(ctdstart<=curmjdstart+t/24);
    futurestni=find(ctdstart>=curmjdstart+t/24);
    presentstni=intersect(previousstni,find(ctdend>=curmjdstart+t/24));
    if (~isempty(presentstni)),    % set nn = num of ongoing station
      nn=presentstni; 
    elseif (futurestni(1)>nn),   % set nn = num of upcoming stn
      nn=futurestni(1);
    elseif (length(previousstni)==length(ctdlist)),  % simstart is past all ctd obs
      nn=length(previousabei);
    end;
  end;
  if (nn<=length(ctdlist) & curmjdstart+t/24>=ctdstart(nn) & curmjdstart+t/24<=ctdend(length(ctdlist)) ),
   ctdboxx=[ctdloc(nn,1) ctdloc(nn,3) ctdloc(nn,5) ctdloc(nn,7) ctdloc(nn,1)];
   ctdboxy=[ctdloc(nn,2) ctdloc(nn,4) ctdloc(nn,6) ctdloc(nn,8) ctdloc(nn,2)];
   ctdlabel=['CTD ',int2str(ctdlist(nn))];
   if (ll==0 | ll==3),
     subplot(spp),boxp=plot(ctdboxx,ctdboxy,'g','LineWidth',2,'Clipping','on');
     subplot(spp),boxt=text(ctdboxx(1),ctdboxy(1),int2str(ctdlist(nn)));
     subplot(spp),btl=text(emax-(emax-wmax)/4,-(nmax-smax)/2/2,ctdlabel);
     %plot points beyond the model domain to prevent resizing margins
     subplot(spp),plot([10*emax 10*wmax 0 0],[0 0 10*nmax 10*smax],'w.');
     if ll==0,  % then set flag to plot 1st ABE observation
       ll=2;
     else,      % an ABE ob has already been plotted, so set flag to plot 2ndary ABE obs 
       ll=4;
     end;
   else,
     %subplot(spp), hold off 
     subplot(spp), set(boxp,'XData',ctdboxx,'YData',ctdboxy);
     subplot(spp), set(boxt,'Position',[ctdboxx(1) ctdboxy(1)],'String',int2str(ctdlist(nn)) );
     subplot(spp), set(btl,'String',ctdlabel);
     %subplot(spp), hold on
   end;
   if (nn<length(ctdlist) & curmjdstart+(t+tres)/24>ctdend(nn)),
     nn=nn+1;
   end;
  elseif exist('boxp','var'),  % hide all lines and labels 
    %subplot(spp), hold off 
    subplot(spp), set(boxp,'XData',[],'YData',[]);
    subplot(spp), set(btl,'String',[]);
    subplot(spp), set(boxt,'String',[]);
    %subplot(spp), hold on
  end;

   if (mm==1),  % if simulation is in first t-step, find # of next or present abe ob
    previousabei=find(abe_mjd_loc(:,7)<=curmjdstart+t/24);
    futureabei=find(abe_mjd_loc(:,8)>=curmjdstart+t/24);
    abelist=abe_mjd_loc(:,1);
    presentabei=intersect(previousabei,futureabei);
    if (~isempty(presentabei)),
      mm=presentabei;
    elseif (futureabei(1)>mm),
      mm=futureabei(1);
    elseif (length(previousabei)==length(abelist)),  % simstart is past all abe obs
      mm=length(previousabei);
    end;
   end;
   if(mm<=length(abe_mjd_loc) & curmjdstart+t/24>=abeobstart(mm) & curmjdstart+t/24<=abeobend(length(abeobend))),
    divetype=abe_mjd_loc(mm,3);
    % eventually load dive nav here, plot ABE position at each time step, 
    % and simulate the ABE time/space T series for comparison with field series... 
    if (ll==0 | ll==2),       %This is the first plot of an ABE ob
     %subplot(spp), hold off
     switch(divetype)
      case 1, type='T'; subplot(spp), boxp2=plot([wside eside eside wside wside],[nside nside sside sside nside],'m','LineWidth',2);
      case 2, type='N'; subplot(spp), boxp2=plot([wside eside],[nside nside],'m','LineWidth',2);
      case 3, type='E'; subplot(spp), boxp2=plot([eside eside],[nside sside],'m','LineWidth',2);
      case 4, type='S'; subplot(spp), boxp2=plot([wside eside],[sside sside],'m','LineWidth',2);
      case 5, type='W'; subplot(spp), boxp2=plot([wside wside],[nside sside],'m','LineWidth',2);
     end;
     abelabel=['ABE ',int2str(abe_mjd_loc(mm,1)),' ',type];
     %btl2=text(wmax+(emax-wmax)/20,smax-(smax-nmax)*1.5/10,abelabel,'FontSize',14);
     subplot(spp),btl2=text(emax-(emax-wmax)/4,-(nmax-smax)/2/1.5,abelabel);
     %subplot(spp), hold on
     if ll==0,  % This was the first ob plotted, so set flag to plot initial CTD ob
       ll=3;
     else, 
       ll=4;
     end;
    else,
     %subplot(spp), hold off
     switch(divetype)
      case 1, type='T'; subplot(spp), set(boxp2,'XData',[wside eside eside wside wside],'YData',[nside nside sside sside nside]);
      case 2, type='N'; subplot(spp), set(boxp2,'XData',[wside eside],'YData',[nside nside]);
      case 3, type='E'; subplot(spp), set(boxp2,'XData',[eside eside],'YData',[nside sside]);
      case 4, type='S'; subplot(spp), set(boxp2,'XData',[wside eside],'YData',[sside sside]);
      case 5, type='W'; subplot(spp), set(boxp2,'XData',[wside wside],'YData',[nside sside]);
     end;
     abelabel=['ABE ',int2str(abe_mjd_loc(mm,1)),' ',type];
     subplot(spp), set(btl2,'String',abelabel);
     %subplot(spp), hold on
    end; 
    if(mm<length(abe_mjd_loc) & curmjdstart+(t+tres)/24>abeobend)
     mm=mm+1;
    end;
   elseif exist('boxp2','var'),  % hide lines and labels
     %subplot(spp), hold off;
     subplot(spp), set(boxp2,'XData',[],'YData',[]);
     subplot(spp), set(btl2,'String',[]); 
     %subplot(spp), hold on;
   end;
  end; % end oblabel if
 end; % end puffplot
 
 if (printframes==1),
  if mofit==1,
    dirnm=[sim,'_',currsrc,int2str(mofmean)];
  else,
    dirnm=[sim,'_',currsrc,int2str(simstart),'-',int2str(simend)];
  end;
  if ~exist(dirnm,'dir'),
    dirstr=['!mkdir ', dirnm];
    eval(dirstr);
  elseif rmflag==1,
      rmok=input('ok to erase eps files? [y,n] ','s')
      if rmok=='y',
        rmstr=['!rm ', dirnm, '/*.eps'];
        eval(rmstr)
      else,
        confirm='files not erased...'
      end;
      rmflag=0;
  end;
  if (j<10), 
   num=['000',int2str(j)]
  elseif (j<100), 
   num=['00',int2str(j)]
  elseif (j<1000),
   num=['0',int2str(j)]
  else
   num=int2str(j)
  end;
  mjd_this_step=curmjdstart+(j-1)*tres/24;
  fn=[dirnm,'/',currsrc,'_frm',num,'_',sprintf('%5.3f',mjd_this_step),'.eps'];
  print( gcf, '-depsc2', fn )
 end;

 % ADVECT THE PUFFS
 if (data == 'real'),
   unow=urot(j);
   vnow=vrot(j);
 else,
   % calculate u,v if not using real data
   unow=udrift+umax*sin(omega*t);
   vnow=vdrift+vmax*sin(omega*t);
 end;

 if (pt_c==1),
  c_ptseries=[c_ptseries; cpt];
  if mofit==1,
    dirnm=[sim,'_',currsrc,int2str(mofmean),'_',pt_name];
  else,
    dirnm=[sim,'_',currsrc,int2str(simstart),'-',int2str(simend),'_',pt_name];
  end;
 end;
 
 if (fluxcalc==1),
  t_elapsed=[0:tres:max(t)];
  % calculate simulated heat fluxes
  dA=dz*(hres*1000);  % m^2, not km^2!! ?? SHOULD dz=VRES or 2*dz?
  %% 12/03/01 NOTE: Need to convert u,v to volume orthogonal components...
  %% 01/21/02 NOTE: Did this conversion above, directly after loading cm data
  eflux=cp*rho*dA*c(ceind)*unow;  % recall u,v units are m/s
  wflux=cp*rho*dA*c(cwind)*-unow; % AND flux is positive OUTWARDS
  nflux=cp*rho*dA*c(cnind)*vnow;
  sflux=cp*rho*dA*c(csind)*-vnow;
  newflux=sum(eflux)+sum(wflux)+sum(nflux)+sum(sflux);
  n_newflux=sum(nflux);
  s_newflux=sum(sflux);
  w_newflux=sum(wflux);
  e_newflux=sum(eflux);
  netflux=[netflux newflux];
  n_totflux=[n_totflux n_newflux];
  s_totflux=[s_totflux s_newflux];
  e_totflux=[e_totflux e_newflux];
  w_totflux=[w_totflux w_newflux];
  c_nseries=[c_nseries; c(cnind)];
  c_sseries=[c_sseries; c(csind)];
  c_eseries=[c_eseries; c(ceind)];
  c_wseries=[c_wseries; c(cwind)];
  if mofit==1,
    dirnm=[sim,'_',currsrc,int2str(mofmean)];
  else,
    dirnm=[sim,'_',currsrc,int2str(simstart),'-',int2str(simend)];
  end;
  if ~exist(dirnm,'dir'),
    dirstr=['!mkdir ', dirnm];
    eval(dirstr);
  end;
  savestr=['save ', dirnm, '/', dirnm, '.mat netflux n_totflux s_totflux e_totflux w_totflux c_nseries c_sseries c_eseries c_wseries t_elapsed simstart tres mofmean urot vrot'];
  eval(savestr)
 end;
 

 if (fluxplot==1),
  % plot the net heat flux through the perimeter
  figure(2)
  t_elapsed=[0:tres:max(t)];
  mjd_elapsed=simstart+(t_elapsed/24);
  recent=find(t_elapsed>t_elapsed(length(t_elapsed))-24);
  plot(t_elapsed(recent),netflux(recent)/1e6,'b')
  grid on;
  hold on;
  plot(t_elapsed(recent),n_totflux(recent)/1e6,'m')
  plot(t_elapsed(recent),s_totflux(recent)/1e6,'c')
  plot([min(t_elapsed(recent)) max(t_elapsed(recent))],[qi/1e6 qi/1e6],'k')
  ylabel('Heat flux, MW')
  xlabel('Elapsed time, hr')
  legend('Net','North side','South side','Input')
  hold off;
 end;

 if (vdthplot==1),
  %% make a quick plot of v versus c
  figure(3)
  grid;
  vseries=[vseries vnow];
  meancseries=[meancseries mean(c(cnind))];
  plot(vseries,meancseries)
  plot(vseries,meancseries,'r.')
 end;

 %% find and eliminate the puffs that are too far away or old
 range=sqrt(xp.^2+yp.^2);
 live=find(range<=abs(emax*near_multiplier) & agep<=maxage);
 xp=[xp(live)+tres*3600*unow/1000; xsrc];   %recall that x,y units are KM!
 yp=[yp(live)+tres*3600*vnow/1000; ysrc];   % tres in hours! u,v in m/s
 agep=[agep(live)+tres; agesrc];
 npuffs=length(live)+1;

 if (displace==1),
  %% note positions of the most southerly, etc puff
  southest=[southest; min(yp-ysrc)];
  northest=[northest; max(yp-ysrc)];
  westest=[westest; max(xp-xsrc)];
  eastest=[eastest; min(xp-xsrc)];
 end;

 if 0,
  if (floor(k/20)>=1),
   ready=input('\n\nAre you ready to add another 20 puffs? (y/n) ','s')
   if (ready~='n'),
     k=1;
   else,
    kdbstop;
   end;
  else,
   k=k+1;
  end;
 end;

 j=j+1;
end;

if (displace==1),
figure(4)
clf
hold
grid
plot([1:tres:t+1],southest,'b')
plot([1:tres:t+1],northest,'c')
plot([1:tres:t+1],westest,'r')
plot([1:tres:t+1],eastest,'m')
legend('southward','northward','westward','eastward',3);
xlabel('Elapsed time, hrs')
ylabel('displacement, km')
title('Time series of displacement from puff source');
end;

if pt_c==1,
 figure(5);
 clf; 
 t_elapsed=[0:tres:max(t)];
 
 if (pt_name=='Wasp2161'),
 % t_elapsed=[0:tres:(tres*length(c_ptseries))-1*tres];
 compind=find(wasp2161.mjd>=simstart & wasp2161.mjd<=simstart+max(t_elapsed)/24);
 plot(wasp2161.mjd(compind),wasp2161.Te(compind),'g')
 hold on
 grid on
 Toffset=-mean(c_ptseries)+mean(wasp2161.Te(compind));
 offlabel=sprintf('%1.2g',Toffset);
 % plot model tseries - its mean + observed mean
 plot(simstart+t_elapsed/24,c_ptseries+Toffset,'b')
 axis([simstart simstart+max(t_elapsed)/24 min(wasp2161.Te) max(wasp2161.Te)]);
 relabelx(5,0)
 title('Modeled and observed T variations at the Wasp2161 mooring')
 xlabel('time (MJD)')
 ylabel('Temperature (^oC)')
 legend('Wasp2161 Expanded T',['Modeled anomaly + ' offlabel])
 print -depsc2 Wasp2161.dT.model.obs.tseries.eps2
 % plot scatter of model and field observations
 figure(50)
 clf
 wasp_c=interp1q(wasp2161.mjd(compind),wasp2161.Te(compind),(simstart+t_elapsed/24)');
 % could insert xcorr analysis here, get lag, and then shift signals prior to time series analysis...
 plot(wasp_c,c_ptseries,'b.')
 grid on
 hold on
 title('Scatter plot of Wasp2161 expanded T and modeled wasp \Delta\theta')
 [Tfit,Sfit]=polyfit(wasp_c,c_ptseries,1); % linear regression
 c_fit=polyval(Tfit,wasp_c,Sfit);
 plot(wasp_c,c_fit,'k')
 print -depsc2 Wasp2161.dT.model.obs.scatter.eps2
 print -djpeg Wasp2161.dT.model.obs.scatter.jpg
 % consider a xcorr plot: plot(xcorr(wasp_c,c_ptseries,'coeff')), but ensure interp1q did not produce NaNs!
 save Wasp2161.dT.model.obs.tseries.mat c_ptseries pt_name simstart tres currsrc
 end;

 if (pt_name=='NoMEFStn'),
 %plot(simstart+(t_elapsed(1:length(t_elapsed)-1)')/24,c_ptseries)
 plot(simstart+t_elapsed'/24,c_ptseries)
 %plot(simstart+t_elapsed/24,c_ptseries)
 relabelx(5,0);
 hold on; grid on;
 load /shuksan/htdocs/flowmow/processing/ctd/Depth_bins_of_dth/tseries.meandth.nomef.mat
 plot(tmeant,tmeandth,'r*')
 relabelx(5,0)
 title('Modeled and observed temperature anomalies North of the MEF')
 xlabel('time (MJD)')
 ylabel('Temperature (^oC)')
 legend('Modeled anomaly','Observed mean anomaly: Stns 3,5,14,28')
 print -depsc2 NoMEF.dTh.model.obs.tseries.eps2
 save NoMEF.dTh.model.obs.tseries.mat c_ptseries pt_name simstart tres currsrc
 end;

if (pt_name=='MZ_cm_95'),
 t_elapsed=[0:tres:(tres*length(c_ptseries))-1*tres];
 compind=find(cm2175.mjd>=simstart & cm2175.mjd<=simstart+max(t_elapsed)/24);
 plot(cm2175.mjd(compind),cm2175.expT(compind),'g')
 hold on
 grid on
 Toffset=-mean(c_ptseries)+mean(cm2175.expT(compind));
 offlabel=sprintf('%1.2g',Toffset);
 % plot model tseries - its mean + observed mean
 plot(simstart+t_elapsed/24,c_ptseries+Toffset,'b')
 axis([simstart simstart+max(t_elapsed)/24 min(cm2175.expT) max(cm2175.expT)]);
 relabelx(5,0)
 title('Modeled and observed T variations at the MZ 25mab mooring')
 xlabel('time (MJD)')
 ylabel('Temperature (^oC)')
 legend('MZcm 25mab Expanded T',['Modeled anomaly + ' offlabel])
 print -depsc2 MZcm.25mab.dT.model.obs.tseries.eps2
 save MZcm.dT.model.obs.tseries.mat c_ptseries pt_name simstart tres currsrc
 % Look for correlation
 % Te2175=interp1q(cm2175.mjd(compind),cm2175.expT(compind),(simstart+t_elapsed/24)');
 % figure; plot(c_ptseries,Te2175,'.')
 % xlabel('Modeled \Delta\theta'); ylabel('Observed Expanded T');
 % title('Modeled vs observed T variations at 1995 current meter, 25mab')
 end;

end;

if (endfluxplot==1),
  % plot the final simulated net/n/s heat fluxes and overlay ctd/abe coverage times
  figure(2)
  clf;
  tit=['Simulated MEF heat flux, W, ',currsrc,':~',int2str(simstart),'-',int2str(simend)];
  mjd_elapsed=simstart+(t_elapsed/24);
  subplot(3,1,1), plot(mjd_elapsed,netflux)
  hold on; grid on;
  title(tit);
  % plot_abe_times
  prevdive=0;
  for kk=1:length(abe_mjd_loc),
    divetype=abe_mjd_loc(kk,3);
    switch(divetype)
     case 1, type='T';
     case 2, type='N';
     case 3, type='E';
     case 4, type='S';
     case 5, type='W';
    end;
    divestart=abe_mjd_loc(kk,7);
    diveend  =abe_mjd_loc(kk,8);
    yaxmax=max(netflux);
    plot([divestart diveend],[0.05*yaxmax 0.05*yaxmax],'c')
    plot([divestart diveend],[0.05*yaxmax 0.05*yaxmax],'g.')
    if (abe_mjd_loc(kk,1)>prevdive),
      divenm=int2str(abe_mjd_loc(kk,1));
      prevdive=abe_mjd_loc(kk,1);
      text( (divestart+diveend)/2,0,divenm)
    end;
    text( (divestart+diveend)/2,0.1*yaxmax,[type int2str(abe_mjd_loc(kk,4))] )
  end;
  for jj=1:length(ctdlist),
    stnm=int2str(ctdlist(jj));
    plot([ctdstart(jj) ctdend(jj)],[0 0],'m')
    plot([ctdstart(jj) ctdend(jj)],[0 0],'r.')
    text( (ctdend(jj)+ctdstart(jj))/2,0,stnm)
  end;
  relabelx(5,0);
  ylabel('Net')

  subplot(3,1,2), plot(mjd_elapsed,n_totflux)
  hold on; grid on;
  prevdive=0;
  for kk=1:length(abe_mjd_loc),
    divetype=abe_mjd_loc(kk,3);
    switch(divetype)
     case 1, type='T';
     case 2, type='N';
     case 3, type='E';
     case 4, type='S';
     case 5, type='W';
    end;
    divestart=abe_mjd_loc(kk,7);
    diveend  =abe_mjd_loc(kk,8);
    yaxmax=max(n_totflux);
    plot([divestart diveend],[0.05*yaxmax 0.05*yaxmax],'c')
    plot([divestart diveend],[0.05*yaxmax 0.05*yaxmax],'g.')
    if (abe_mjd_loc(kk,1)>prevdive),
      divenm=int2str(abe_mjd_loc(kk,1));
      prevdive=abe_mjd_loc(kk,1);
      text( (divestart+diveend)/2,0,divenm)
    end;
    text( (divestart+diveend)/2,0.1*yaxmax,[type int2str(abe_mjd_loc(kk,4))] )
  end;
  for jj=1:length(ctdlist),
    stnm=int2str(ctdlist(jj));
    plot([ctdstart(jj) ctdend(jj)],[0 0],'m')
    plot([ctdstart(jj) ctdend(jj)],[0 0],'r')
    text((ctdend(jj)+ctdstart(jj))/2,0,stnm)
  end;
  relabelx(5,0);
 if 0,
  % overlay observed fluxes in Watts from N endzone
  A=(eside-wside)*1000*dz;
  qfact=cp*rho*A;
  if (currsrc=='cm050')
    load /shuksan/htdocs/flowmow/processing/ctd/Endzones/cm50.n_endzone.va.meanDth.meanMJD.mat 
    plot(mjdmeanall,v50all.*Dthmeanall*qfact,'y')
    plot(mjdmeanall,v50all.*Dthmeanall*qfact,'k*')
    %legend('Northern observed heat flux time series, cm50')
   elseif(currsrc=='wasp0'),
    load /shuksan/htdocs/flowmow/processing/ctd/Endzones/wasp2161.n_endzone.va.meanDth.meanMJD.mat
    plot(mjdmeanall,vwaspall.*Dthmeanall*qfact,'y')
    plot(mjdmeanall,vwaspall.*Dthmeanall*qfact,'k*')
    %legend('Northern observed heat flux time series, wasp')
  end;
 end;
  ylabel('North side')

  subplot(3,1,3), plot(mjd_elapsed,s_totflux)
  hold on; grid on;
  prevdive=0;
  for kk=1:length(abe_mjd_loc),
    divetype=abe_mjd_loc(kk,3);
    switch(divetype)
     case 1, type='T';
     case 2, type='N';
     case 3, type='E';
     case 4, type='S';
     case 5, type='W';
    end;
    divestart=abe_mjd_loc(kk,7);
    diveend  =abe_mjd_loc(kk,8);
    yaxmax=max(s_totflux);
    plot([divestart diveend],[0.05*yaxmax 0.05*yaxmax],'c')
    plot([divestart diveend],[0.05*yaxmax 0.05*yaxmax],'g.')
    if (abe_mjd_loc(kk,1)>prevdive),
      divenm=int2str(abe_mjd_loc(kk,1));
      prevdive=abe_mjd_loc(kk,1);
      text( (divestart+diveend)/2,0,divenm)
    end;
    text( (divestart+diveend)/2,0.1*yaxmax,[type int2str(abe_mjd_loc(kk,4))] )
  end;
  for jj=1:length(ctdlist),
    stnm=int2str(ctdlist(jj));
    plot([ctdstart(jj) ctdend(jj)],[0 0],'m')
    plot([ctdstart(jj) ctdend(jj)],[0 0],'r.')
    text((ctdend(jj)+ctdstart(jj))/2,0,stnm)
  end;
  relabelx(5,0);
  ylabel('South side')
  xlabel('Elapsed time, hr')
end;

%end; %mofmean loop

%mof_display