function y=contourt(tc,y,z,year)
fpath=['/data/span9/bahr/perlprogs/'];
x=tc;
mond=[0,32,60,91,121,152,182,213,244,275,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;
fid=fopen([fpath,'contours.ppl'],'w');
fprintf(fid,'pltype,0\n');
fprintf(fid,'pltnme,%scontours.meta\n',fpath(:));
fprintf(fid,'multplt,1,3\n');
fprintf(fid,'5.0\n');
fprintf(fid,'2.0,2.0,2.0\n');
fprintf(fid,'1.5\n');
fprintf(fid,'-2.0,-2.0,1.5\n');
fprintf(fid,'dfltfnt,@tr\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,300,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
dayst=floor(ztmin);
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
zzmin=min(z);
zzmax=max(z);
% base increment on the difference between zzmin and zzmax
dz=zzmax-zzmin;
% set zst and zsp to the nearest 1/10
zst=floor(zzmin*10)/10;
zsp=ceil(zzmax*10)/10;
fprintf(fid,'lev,arg()\n');
if (dz >= 0.8)&(dz < 1.6),
	ccinc=0.2;
	ks=rem(zst,1)*10;
	kp=rem(zsp,1)*10;
	if (rem(ks,2)) == 1,
		zst=zst-0.1;
	end
	if (rem(kp,2)) == 1,
		zsp=zsp+0.1;
	end
	zc=[zst:ccinc:zsp];
	nc=length(zc);
	fprintf(fid,'lev,arg(%4.1f,%4.1f,%4.1f,1)\n',zst,zsp,ccinc);
	fprintf(fid,'cbaxis,%4.1f,%4.1f,%4.1f\n',zst,zsp,ccinc);
	fprintf(fid,'cblint,1\n');
	fprintf(fid,'cbfor,(f4.1)\n');
end
if dz < 0.8,
	ccinc=0.1;
	zc=[zst:ccinc:zsp];
	nc=length(zc);
	fprintf(fid,'lev,arg(%4.1f,%4.1f,%4.1f,1)\n',zst,zsp,ccinc);
	fprintf(fid,'cbaxis,%4.1f,%4.1f,%4.1f\n',zst,zsp,ccinc);
	fprintf(fid,'cblint,1\n');
	fprintf(fid,'cbfor,(f4.1)\n');
end
if (dz >= 1.6)&(dz < 4),
	ccinc=0.5;
	ks=rem(zst,1)*10;
	kp=rem(zsp,1)*10;
	if rem(ks,5) > 0.1,
		ks=ks-rem(ks,5);
		zst=floor(zst);
		zst=zst+ks/10;
	end
	if rem(kp,5) > 0.1,
		kp=kp+(5-rem(kp,5));
		zsp=floor(zsp);
		zsp=zsp+kp/10;
	end
	zc=[zst:ccinc:zsp];
	nc=length(zc);
	fprintf(fid,'lev,arg(%4.1f,%4.1f,%4.1f,1)\n',zst,zsp,ccinc);
	fprintf(fid,'cbaxis,%4.1f,%4.1f,%4.1f\n',zst,zsp,ccinc);
	fprintf(fid,'cblint,1\n');
	fprintf(fid,'cbfor,(f4.1)\n');
end
if dz >= 4,
	ccinc=1;
	zst=floor(zzmin); zsp=ceil(zzmax);
	zc=[zst:ccinc:zsp];
	nc=length(zc);
	fprintf(fid,'lev,arg(%d,%d,%d,0)\n',zst,zsp,ccinc);
	fprintf(fid,'cbaxis,%d,%d,%d\n',zst,zsp,ccinc);
	fprintf(fid,'cblint,1\n');
	fprintf(fid,'cbfor,(i3)\n');
end
fprintf(fid,'loadcmap,%smycolor.16\n',fpath(:));
% 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,300,1,2\n',xmax+0.01);
% 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);
ymin=10;
ymax=300;
% 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,%sms.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,'area\n');
fprintf(fid,'contour\n');
fprintf(fid,'axlabp,-1,0\n');
fprintf(fid,'axset,1,1,0,0\n');
fprintf(fid,'taxis,10,on\n');
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=rem(year,100);
yy=year;
if tlow < 0,
	if ileap==1,
		tlow=tlow+367;
	else,
		tlow=tlow+366;
	end
end
for i=1:12,
	if (thig >mond(i))&(thig <= mond(i+1)),
		moh=i;
		dah=fix(thig-mond(i))+1;
		hoh=rem(thig,1);
		hoh=fix(hoh*24);
	end
	if (tlow > mond(i))&(tlow <= mond(i+1)),
		mol=i;
		dal=fix(tlow-mond(i))+1;
		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;
%	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,Salinity\n');
fprintf(fid,'exit\n');
y=0;
