function y=contourt(tc,y,u,v,year,mid)
% Program to write out the script files for Plotplus
fpath=['/data/span9/bahr/perlprogs/'];
x=tc;
mond=[0,32,60,91,121,152,182,213,244,274,305,335,366];
if(rem(year,4)==0),
	leap=[0 0 1 1 1 1 1 1 1 1 1 1 1];
	ileap=1;
else,
	leap=[0 0 0 0 0 0 0 0 0 0 0 0 0];
	ileap=0;
end
mond=mond+leap;
% before we get going to far check to see if wind file is there
if mid ~=5,
	eval('load mooringid.dat');
	mid=mooringid;
end
ffd=fopen([fpath,'wind',num2str(mid),'.dat'],'r');
if ffd > 2,
	fclose(ffd);
	fwin=1;
else,
%	fclose(ffd);
	fwin=0;
end
fid=fopen([fpath,'contouruv.ppl'],'w');
fprintf(fid,'pltype,0\n');
fprintf(fid,'pltnme,%scontouruv.meta\n',fpath(:));
% determine using flags above, wether to plot wind or not
if fwin==0,
	fprintf(fid,'multplt,1,6\n');
	fprintf(fid,'6.2\n');
	fprintf(fid,'2.0,2.0,2.0,2.0,2.0,2.0\n');
	fprintf(fid,'1\n');
	fprintf(fid,'-2.0,-2.0,0.5,-2.0,-2.0,1\n');
else,
	fprintf(fid,'multplt,1,7\n');
	fprintf(fid,'6.2\n');
	fprintf(fid,'2.0,2.0,2.0,2.0,2.0,2.0,1.5\n');
	fprintf(fid,'1\n');
	fprintf(fid,'-2.0,-2.0,0.5,-2.0,-2.0,0.25,1\n');
end
fprintf(fid,'dfltfnt,@tr\n');
fprintf(fid,'size,11.0,11.0\n');
if fwin==1,
	fprintf(fid,'limits,-99,XLE,on\n');
	fprintf(fid,'limits,-99,YLE,on\n');
end
fprintf(fid,'window,on\n');
fprintf(fid,'axlabp,0,-1\n');
fprintf(fid,'axlint,1,1\n');
fprintf(fid,'axnmtc,1,1\n');
fprintf(fid,'axset,0,0,1,1\n');
fprintf(fid,'ylab,@trDepth(m)\n');
% set up the yaxis for depth from surface to 300 m
	ymax=120; ymin=0; dy=-50;
yt=max(y);
if yt < ymax-10,
	yt=yt/10;
	yt=ceil(yt);
	yf=rem(yt,10);
	da=floor(yf/5);
	yh=floor(yt/10);
	if da==1,
		yh=yh+1;
		yf=0;
	else,
		yf=5;
	end
	ymax=yh*100+yf*10;
end
if yt > ymax+10,
	yt=yt/10;
	yt=ceil(yt);
	yf=rem(yt,10);
	da=floor(yf/5);
	yh=floor(yt/10);
	if da==1,
		yh=yh+1;
		yf=0;
	else,
		yf=5;
	end
	ymax=yh*100+yf*10;
end
if ymax > 300,
	dy=-100;
end
fprintf(fid,'yaxis,%d,%d,%d\n',ymax,ymin,dy);
% we have to be careful that the x limits match those
% we have to use for the taxis command
ztmin=min(tc);
ztmax=max(tc);
ztmaxr=ztmax;
if (ztmax-ztmin) < 7,
	ztmax=ztmin+7;
end
if ztmin < 0,
	dayst=ceil(ztmin);
else,
	dayst=floor(ztmin);
end
daysp=floor(ztmax);
hrst=rem(ztmin,1)*24;
hrsp=rem(ztmax,1)*24;
hrst=floor(hrst);
hrsp=ceil(hrsp);
if hrsp > 24,
	daysp=daysp+1;
	hrsp=hrsp-24;
end
xmin=dayst+hrst/24;
xmax=daysp+hrsp/24;
tlow=xmin;
thig=xmax;
fprintf(fid,'xaxis,%9.5f,%9.5f,1\n',xmin,xmax);
% compute contour interval stuff
ll=u > -3000;
tu=u(ll);
ll=tu < 3000;
tu=tu(ll);
zumin=min(min(tu));
zumax=max(max(tu));
ll=v > -3000;
tv=v(ll);
ll=tv < 3000;
tv=tv(ll);
zvmin=min(min(tv));
zvmax=max(max(tv));
if zvmin <= zumin,
	zzmin=zvmin;
else,
	zzmin=zumin;
end
if zvmax >= zumax,
	zzmax=zvmax;
else,
	zzmax=zumax;
end
% base increment on the difference between zzmin and zzmax
dz=zzmax-zzmin;
% set zst and zsp to the nearest integers
zst=floor(zzmin);
zsp=ceil(zzmax);
fprintf(fid,'lev,arg()\n');
% figure out what to use for a contour interval
% Since we are working with velocites in the cm/s range with
% concievable ranges from -50 to 50  
fprintf(fid,'loadcmap,%smycolor.16\n',fpath(:));
if (dz >= 50),
	ccinc=10;
	zst=floor(zst/10)*10;
	zsp=ceil(zsp/10)*10;
	zc=[zst:ccinc:zsp];
	nc=length(zc);
	fprintf(fid,'lev,arg(%d,%d,%d,-1)\n',zst,zsp,ccinc);
	fprintf(fid,'cbaxis,%d,%d,%d\n',zst,zsp,ccinc);
	fprintf(fid,'cblint,1\n');
	fprintf(fid,'cbfor,(i4)\n');
end
if (dz < 50),
	ccinc=5;
	ks=rem(zst,10);
	kp=rem(zsp,10);
	ks=ks-rem(ks,5);
	zst=floor(zst/10);
	zst=zst*10-ks;
	kp=kp+(5-rem(kp,5));
	zsp=floor(zsp/10);
	zsp=zsp*10+kp;
	zc=[zst:ccinc:zsp];
	nc=length(zc);
	fprintf(fid,'lev,arg(%d,%d,%d,-1)\n',zst,zsp,ccinc);
	fprintf(fid,'cbaxis,%d,%d,%d\n',zst,zsp,ccinc);
	fprintf(fid,'cblint,1\n');
	fprintf(fid,'cbfor,(i4)\n');
end
% for at test we will put in where the data points are
%fprintf(fid,'datpt,3,3\n');
% we need to compute where to put the color bar
fprintf(fid,'colorbar,%f,%d,1,2\n',xmax+0.01,ymax);
% need to compute the levels for the contour intervals we
% plan to use.
% compute values for plot min,max and number of triplets
% to read
lx=size(x,1)*size(u,2);
% set up contour grid (7days*144pts/day, 2*10pts in depth)
% actually try to compute grid size.  Currently hardwire
% y values to 10pts in vertical.  x<time> we can set...
dx=ztmaxr-xmin;
xgrid=floor(dx*24);
fprintf(fid,'skp,0,%suv.dat\n',fpath(:));
fprintf(fid,'format,free\n');
fprintf(fid,'vars,1,1,2,3\n');
fprintf(fid','rd,%d,20\n',xgrid);
% number of points to read
fprintf(fid,'%d\n',lx);
% we will have to carefully define xmin and xmax to correspond
% here with what will be used later by the taxis command...
fprintf(fid,'%9.5f,%9.5f,%d,%d\n',xmin,ztmaxr,ymin,ymax);
fprintf(fid,'labs,\n');
zrat=(-10*ymax)/120;
fprintf(fid,'labs,1,%f,%f,0,Eastward Velocity (cm/s)\n',(xmin+xmax)/2,zrat);
fprintf(fid,'area\n');
fprintf(fid,'contour\n');
fprintf(fid,'axlabp,-1,0\n');
fprintf(fid,'axset,1,1,0,0\n');
fprintf(fid,'labs,\n');
fprintf(fid,'taxis,10,on\n');
fprintf(fid,'txnmtc,-1\n');
fprintf(fid,'txlabp,0\n');
fprintf(fid,'txlint,24,2\n');
fprintf(fid,'txtype,DAY,HR\n');
% Now compute the start and stop time limits for the plot
% yy=rem(year,100);
% New Y2K compliant version of pplus allows for 4 digit year.
yy=year;
if tlow < 0,
	if ileap==1,
		tlow=tlow+366;
	else,
		tlow=tlow+365;
	end
end
for i=1:12,
	if (thig >mond(i))&(thig <= mond(i+1)),
		moh=i;
		if i==1,
			dah=fix(thig-mond(i));
		else,
			dah=fix(thig-mond(i))+1;
		end
		hoh=rem(thig,1);
		hoh=fix(hoh*24);
	end
	if (tlow > mond(i))&(tlow <= mond(i+1)),
		mol=i;
		if (i==1)&(dah >= 8),
			dal=fix(tlow-mond(i));
		else,
			dal=fix(tlow-mond(i))+1;
		end
		hol=rem(tlow,1);
		mil=hol;
		hol=fix(hol*24);
		mil=mil*24-hol;
		mil=fix(mil*60);
	end
end
if(moh==1)&(dah < 7),
	yo=year-1;
% Leave above alone and remove below because of Y2K 
% compliant results
%	yo=rem(yo,100);
else,
	yo=yy;
end
if mol < 10,
	mol=['0',num2str(mol)];
else,
	mol=num2str(mol);
end
if moh < 10,
	moh=['0',num2str(moh)];
else,
	moh=num2str(moh);
end
if dal < 10,
	dal=['0',num2str(dal)];
else,
	dal=num2str(dal);
end
if dah < 10,
	dah=['0',num2str(dah)];
else,
	dah=num2str(dah);
end
if hol < 10,
	hol=['0',num2str(hol)];
else,
	hol=num2str(hol);
end
if hoh < 10,
	hoh=['0',num2str(hoh)];
else,
	hoh=num2str(hoh);
end
if mil < 10,
	mil=['0',num2str(mil)];
else,
	mil=num2str(mil);
end
tstar=['W',num2str(yo),mol,dal,hol,'00'];
tstop=['W',num2str(yy),moh,dah,hoh,'00'];
tstart=['W',num2str(yo),mol,dal,hol,mil];
fprintf(fid,'time,%s,%s,%s\n',tstar(:),tstop(:),tstart(:));
fprintf(fid,'format,free\n');
%fprintf(fid,'ylab,@tr Depth (m)\n');
fprintf(fid,'skp,0,%sdummy.dat\n',fpath(:));
fprintf(fid,'vars,1,1,2\n');
fprintf(fid,'rd\n');
fprintf(fid,'plot\n');
%fprintf(fid,'exit\n');
fprintf(fid,'taxis,off\n');
fprintf(fid,'window,on\n');
fprintf(fid,'axlabp,0,-1\n');
fprintf(fid,'axlint,1,1\n');
fprintf(fid,'axnmtc,1,1\n');
fprintf(fid,'axset,0,0,1,1\n');
fprintf(fid,'ylab,@trDepth(m)\n');
% set up the yaxis for depth from surface to 300 m
fprintf(fid,'yaxis,%d,%d,%d\n',ymax,ymin,dy);
%fprintf(fid,'yaxis,120,0,-50\n');
% we have to be careful that the x limits match those
% we have to use for the taxis command
ztmin=min(tc);
ztmax=max(tc);
ztmaxr=ztmax;
if (ztmax-ztmin) < 7,
	ztmax=ztmin+7;
end
if ztmin < 0,
	dayst=ceil(ztmin);
else,
	dayst=floor(ztmin);
end
daysp=floor(ztmax);
hrst=rem(ztmin,1)*24;
hrsp=rem(ztmax,1)*24;
hrst=floor(hrst);
hrsp=ceil(hrsp);
if hrsp > 24,
	daysp=daysp+1;
	hrsp=hrsp-24;
end
xmin=dayst+hrst/24;
xmax=daysp+hrsp/24;
tlow=xmin;
thig=xmax;
fprintf(fid,'xaxis,%9.5f,%9.5f,1\n',xmin,xmax);
% compute contour interval stuff
%ll=v > -3000;
%tv=v(ll);
%zzmin=min(min(tv));
%zzmax=max(max(tv));
% base increment on the difference between zzmin and zzmax
%dz=zzmax-zzmin;
% set zst and zsp to the nearest integers
%zst=floor(zzmin);
%zsp=ceil(zzmax);
fprintf(fid,'lev,arg()\n');
% figure out what to use for a contour interval
% Since we are working with velocites in the cm/s range with
% concievable ranges from -50 to 50 
if (dz >= 50),
	ccinc=10;
	zst=floor(zst/10)*10;
	zsp=ceil(zsp/10)*10;
	zc=[zst:ccinc:zsp];
	nc=length(zc);
	fprintf(fid,'lev,arg(%d,%d,%d,-1)\n',zst,zsp,ccinc);
	fprintf(fid,'cbaxis,%d,%d,%d\n',zst,zsp,ccinc);
	fprintf(fid,'cblint,1\n');
	fprintf(fid,'cbfor,(i4)\n');
end
if (dz < 50),
	ccinc=5;
	ks=rem(zst,10);
	kp=rem(zsp,10);
%	if rem(ks,5) > 0,
		ks=ks-rem(ks,5);
		zst=floor(zst/10);
		disp([ks zst])
		zst=zst*10-ks;
%	end
%	if rem(kp,5) > 0,
		kp=kp+(5-rem(kp,5));
		zsp=floor(zsp/10);
		zsp=zsp*10+kp;
%	end
	zc=[zst:ccinc:zsp];
	nc=length(zc);
	fprintf(fid,'lev,arg(%d,%d,%d,-1)\n',zst,zsp,ccinc);
	fprintf(fid,'cbaxis,%d,%d,%d\n',zst,zsp,ccinc);
	fprintf(fid,'cblint,1\n');
	fprintf(fid,'cbfor,(i4)\n');
end
% for at test we will put in where the data points are
%fprintf(fid,'datpt,3,3\n');
% we need to compute where to put the color bar
fprintf(fid,'colorbar,%f,%d,1,2\n',xmax+0.01,ymax);
% need to compute the levels for the contour intervals we
% plan to use.
% compute values for plot min,max and number of triplets
% to read
lx=size(x,1)*size(v,2);
%ymin=min(y);
%ymax=max(y);
% set up contour grid (7days*144pts/day, 2*10pts in depth)
% actually try to compute grid size.  Currently hardwire
% y values to 10pts in vertical.  x<time> we can set...
% xmin=min(x);
% xmax=max(x);
dx=ztmaxr-xmin;
xgrid=floor(dx*24);
fprintf(fid,'skp,0,%suv.dat\n',fpath(:));
fprintf(fid,'format,free\n');
fprintf(fid,'vars,1,1,2,,3\n');
fprintf(fid','rd,%d,20\n',xgrid);
% number of points to read
fprintf(fid,'%d\n',lx);
% we will have to carefully define xmin and xmax to correspond
% here with what will be used later by the taxis command...
fprintf(fid,'%9.5f,%9.5f,%d,%d\n',xmin,ztmaxr,ymin,ymax);
fprintf(fid,'labs,\n');
fprintf(fid,'labs,1,%f,%f,0,Northward Velocity (cm/s)\n',(xmin+xmax)/2,zrat);
fprintf(fid,'area\n');
fprintf(fid,'contour\n');
if fwin==0,
	fprintf(fid,'axlabp,-1,0\n');
else,
	fprintf(fid,'axlabp,0,0\n');
end
fprintf(fid,'axset,1,1,0,0\n');
fprintf(fid,'taxis,10,on\n');
fprintf(fid,'labs,\n');
if fwin==0,
	fprintf(fid,'txlabp,-1\n');
else,
	fprintf(fid,'txlabp,0\n');
end
fprintf(fid,'txnmtc,-1\n');
fprintf(fid,'txlint,24,2\n');
fprintf(fid,'txtype,DAY,HR\n');
% Now compute the start and stop time limits for the plot
yy=year;
%Add above line which is y2k compliant below isn't.
%yy=rem(year,100);
% tlow=min(tc);
if tlow < 0,
	if ileap==1,
		tlow=tlow+366;
	else,
		tlow=tlow+365;
	end
end
% thig=max(tc);
for i=1:12,
	if (thig >mond(i))&(thig <= mond(i+1)),
		moh=i;
		if i==1,
			dah=fix(thig-mond(i));
		else,
			dah=fix(thig-mond(i))+1;
		end
		hoh=rem(thig,1);
		hoh=fix(hoh*24);
	end
	if (tlow > mond(i))&(tlow <= mond(i+1)),
		mol=i;
		if (i==1)&(dah >= 8),
		  dal=fix(tlow-mond(i));
		else,
			dal=fix(tlow-mond(i))+1;
		end
		hol=rem(tlow,1);
		mil=hol;
		hol=fix(hol*24);
		mil=mil*24-hol;
		mil=fix(mil*60);
	end
end
if(moh==1)&(dah < 7),
	yo=year-1;
% remove non-y2k compliant code
%	yo=rem(yo,100);
else,
	yo=yy;
end
if mol < 10,
	mol=['0',num2str(mol)];
else,
	mol=num2str(mol);
end
if moh < 10,
	moh=['0',num2str(moh)];
else,
	moh=num2str(moh);
end
if dal < 10,
	dal=['0',num2str(dal)];
else,
	dal=num2str(dal);
end
if dah < 10,
	dah=['0',num2str(dah)];
else,
	dah=num2str(dah);
end
if hol < 10,
	hol=['0',num2str(hol)];
else,
	hol=num2str(hol);
end
if hoh < 10,
	hoh=['0',num2str(hoh)];
else,
	hoh=num2str(hoh);
end
if mil < 10,
	mil=['0',num2str(mil)];
else,
	mil=num2str(mil);
end
tstar=['W',num2str(yo),mol,dal,hol,'00'];
tstop=['W',num2str(yy),moh,dah,hoh,'00'];
tstart=['W',num2str(yo),mol,dal,hol,mil];
fprintf(fid,'time,%s,%s,%s\n',tstar(:),tstop(:),tstart(:));
fprintf(fid,'format,free\n');
%fprintf(fid,'ylab,@tr Depth (m)\n');
fprintf(fid,'skp,0,%sdummy.dat\n',fpath(:));
fprintf(fid,'vars,1,1,2\n');
fprintf(fid,'rd\n');
fprintf(fid,'plot\n');
% now to add all the wind file specific plot stuff
if fwin==1,
	eval(['load ',fpath,'wind',num2str(mid),'.dat']);
	eval(['ww=wind',num2str(mid),';']);
% find starting time
	wint=min(ww(:,1));
	if wint < 0,
		if ileap==1,
			wint=wint+366;
		else,
			wint=wint+365;
		end
	end
	for i=1:12,
		if (wint > mond(i))&(wint <= mond(i+1)),
			mow=i;
			if i==1,
				daw=fix(wint-mond(i));
			else,
				daw=fix(wint-mond(i))+1;
			end
			how=rem(wint,1);
			miw=how;
			how=fix(how*24);
			miw=miw*24-how;
			miw=fix(miw*60);
		end
	end
% Now convert numbers above to strings
	if daw < 10,
		daw=['0',num2str(daw)];
	else,
		daw=num2str(daw);
	end
	if mow < 10,
		mow=['0',num2str(mow)];
	else,
		mow=num2str(mow);
	end
	if how < 10,
		how=['0',num2str(how)];
	else,
		how=num2str(how);
	end
	if miw < 10,
		miw=['0',num2str(miw)];
	else,
		miw=num2str(miw);
	end
	wstart=['W',num2str(yo),mow,daw,how,miw];
% Now to find the y limits for the plot
	amax=10;
	amin=-10;
	wl=ww(:,2) > -99;
	wx=ww(wl,:);
	mau=max(wx(:,2));
	mav=max(wx(:,3));
	miu=min(wx(:,2));
	miv=min(wx(:,3));
	maxvel=sqrt(mau^2+mav^2);
	minvel=sqrt(miu^2+miv^2);
	if maxvel > amax,
		amax=maxvel;
	end
	if minvel < amin,
		amin=minvel;
	end
	ad=rem(amax,5);
	if ad ~=0,
		amax=amax+(5-ad);
	end
	ad=rem(amin,5);
	if ad ~=0,
		amin=amin-(5+ad);
	end
	fprintf(fid,'axlabp,-1,-1\n');
	fprintf(fid,'axset,1,1,1,1\n');
	fprintf(fid,'taxis,10,on\n');
	fprintf(fid,'txlabp,-1\n');
	fprintf(fid,'yaxis,%f,%f,5\n',amin,amax);
	fprintf(fid,'time,%s,%s,%s\n',tstar(:),tstop(:),wstart(:));
	fprintf(fid,'ylab,@trWind Speed (m/s)\n');
	fprintf(fid,'skp,0,%swind%d.dat\n',fpath(:),mid);
	fprintf(fid,'vars,1,,1,2\n');
	fprintf(fid,'rd\n');
	fprintf(fid,'plotv\n');
end
fprintf(fid,'exit\n');
fclose(fid);
fy=0;
