%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% M-file to process buoy m"x" ADCP data.
% Simple bounds checks have already been
% applied to the data which have been output
% from the perl script.  This step is
% to concattenate the data and to create a
% much large files and do more sophisticated (sp?)
% processing of the data.  At this point we will
% rotate the data for magnetic declination and
% attempt to remove mooring motion (lateral) from
% the computed velocities.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define an array
fpath=['/data/span9/bahr/perlprogs/'];
mond=[0,31,59,90,120,151,181,212,243,273,305,334];
monn=[31,28,31,30,31,30,31,31,30,31,30,31];
monc=['jan',
      'feb',
      'mar',
      'apr',
      'may',
      'jun',
      'jul',
      'aug',
      'sep',
      'oct',
      'nov',
      'dec'];
dtt=10.0/60.0/24.0;
ylat=[36.717 36.783
      36.630 36.717
      36.517 36.617
      36.100 36.200];
ylon=[-122.05 -122.00
      -122.45 -122.37
      -123.00 -122.95
      -122.47 -122.40];
% open file with mooring number which we will be working on
mid=input('Enter the mooring number to append ADCP data to: ');
year=input('Enter the year of the data to append: ');
month=input('Enter the month of the data to append: ');
fnpre=input('Enter the filename prefix to use: ','s');
% create filename for the current month
pathm1=['/data/span9/bahr/m',num2str(mid),'/adcp/',num2str(year),'/'];
fnm1=['m',num2str(mid),'_adcp_uv_',monc(month,:),'_',num2str(year)];
fnm2=['m',num2str(mid),'_adcp_sc_',monc(month,:),'_',num2str(year)];
% see if the file exists
fid1=fopen([pathm1,fnm1,'.mat'],'r+');
fid2=fopen([pathm1,fnm2,'.mat'],'r+');
% if the file exists then fid will be greater than 2.  If the file
% doesn't exist we want to create it but we also want to load
% the previous months file.  We also have to careful around the year
% change.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The adcp_uv file will contain the following arrays:
% un-->U array vn-->V array en-->Error velocity array
% cmn-> Correlation magnitude array.
% gn--> AGC array gpsn-->GPS data array (time,lat,lon,stat)
% tn--> Time vector for ADCP data (NOT GPS data).
% fu--> Flag array for U velocities
% fv--> Flag array for V velocities
% fe--> Flag array for error velocity
% fg--> Flag array for AGC
% fgps--> Flag array for GPS info.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if fid1 > 2,
	fclose(fid1);
	eval(['load ',pathm1,fnm1,'.mat']);
% load velocity arrays
	eval(['u=un;']);
	eval(['v=vn;']);
	eval(['e=en;']);
% load correlation magnitude arrays (multidimensional)
	eval(['cm=cmn;']);
% load echo intensity (agc) array (multidimensional)
	eval(['g=gn;']);
% load percent 4 beam good solution array
	eval(['pg4=pg4n;']);
% load percent 3 beam goos solution array
	eval(['pg3=pg3n;']);
% load time vector
	eval(['t=tn;']);
% load gps data array
	eval(['gps=gpsn;']);
% load velocity flags
	eval(['uf=fu;']);
	eval(['vf=fv;']);
% load gps data flags
	eval(['gpsf=fgps;']);
else,
	u=[];
	v=[];
	e=[];
	cm=[];
	g=[];
	pg4=[];
	pg3=[];
	t=[];
	gps=[];
	uf=[];
	vf=[];
	gpsf=[];
end
% Now load the scalar information
% add flag to this for transducer temperature
if fid2 > 2,
	fclose(fid2);
	eval(['load ',pathm1,fnm2,'.mat']);
	eval(['c=adcpscn;']);
	eval(['cf=adcpfsc;']);
else,
	c=[];
	cf=[];
end
adcpsc=[c];
fsc=[cf];
adt=t;
% we should now have the old matlab process adcp files loaded.
% The next step is to load the new files with data to be processed.
% up to this point we have not decide how the arrays are stored in
% the ADCP files.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% So we can handle changes in number of velocities put out by the
% ADCP's we will create arrays which grow in time and are n columns
% wide (where n is the greater bin number) so we will have:
% Depth vector       U array       V array     E array      AGC array by 4
% D1--->D(n)      U1--->U(n)      V1--->V(n)   E1--->E(n)   A1--->A(n)
%                 |               |            |            |  
%                 |               |            |            |
%                 v               v            v            v
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now also load GPS data for mooring motion if the files exist......
% What do we do if there is no GPS data for a mooring....
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ff1=[fnpre(1,:),'_uv.dat'];
fs1=[fnpre(1,:),'_sc.dat'];
fg1=[fnpre(1,:),'_gps.dat'];
fod1=fopen([pathm1,ff1],'r');
fsd1=fopen([pathm1,fs1],'r');
fog1=fopen([pathm1,fg1],'r');
if fod1 > 2,
	fclose(fod1);
	eval(['load ',pathm1,ff1(1,:)]);
	eval(['m1uv=',ff1(1,[1:length(ff1)-4]),';']);
else,
	m1uv=[];
end
if fsd1 > 2,
	fclose(fsd1);
	eval(['load ',pathm1,fs1(1,:)]);
	eval(['m1sc=',fs1(1,[1:length(fs1)-4]),';']);
else,
	m1sc=[];
end
if fog1 > 2,
	fclose(fog1);
	eval(['load ',pathm1,fg1(1,:)]);
	eval(['mgps=',fg1(1,[1:length(fg1)-4]),';']);
else,
	mgps=[];
end
% we have a problem here at the year change.  We need to worry about.
% since the files only look at year day and not have the year in it.
m1=[m1uv];
% m2 contains the ADCP scalar data corresponding to the vector data in m1.
m2=[m1sc];
% mg will contain the gps information
mg=[mgps];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Next we want to get the GPS data sorted and checked out
% At this point we want to get GPS data cleaned up and ready 
% to compute velocities from.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (isempty(mg) == 0),
	dtg=diff(mg(:,1));
	[inegg]=find(dtg < 0);
	if isempty(inegg)==0,
		[newtg indg]=sort(mg(:,1));
		mg=[newtg mg(indg,2:size(mg,2))];
	end
	tmin=min(mg(:,1));
	if isempty(gps)==0,
		tmax=max(gps(:,1));
		cmin=min(gps(:,1));
	else
		tmax=0;
		cmin=370;
	end
	lgp=size(mg,1);
	ig=zeros(lgp,1);
	if isempty(gps)==0,
		for i=1:lgp,
			if mg(i,1) > tmax,
				i=lgp;
			end
			igz=find(gps(:,1)==mg(i,1));
			if isempty(igz)==0,
				ig(i)=1;
			end
		end
	end
	mg(logical(ig),:)=[]; % remove identical times.
end
mymax=ylat(mid,2);
mymin=ylat(mid,1);
mxmax=ylon(mid,2);
mxmin=ylon(mid,1);
if(isempty(mg)==0),
% Now we eliminate data based upon some criteria from the above data.
% Add arrays earlier defining the max/min lat/lon for the mooring positions
% Then flag the data.
% Define a flag array for new gps data.
	fmg=zeros(size(mg));
% Check which values are out of bounds and flag them.
	ll=((mg(:,2) < mymin)|(mg(:,2) > mymax));
	fmg(ll,2)=1;
	ll=((mg(:,3) < mxmin)|(mg(:,3) > mxmax));
	fmg(ll,3)=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ==>fill in here determination if data is good.<==
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now fill in the gps array with the new data....
	i=0; j=0;
	lgp=size(gps,1);
	gmask=zeros(lgp,1);
	lnew=size(mg,1);
	k=1;
	if(cmin==370),
		gps=mg;
		gpsf=fmg(:,[2 3]);
		gmask=[ones(lnew,1)];
		k=lnew+1;
	else,
		if(mg(k,1) < cmin),
			while((mg(k,1) < cmin)&(k < size(mg,1))),
				k=k+1;
			end
			gps=[mg(1:k-1,:); gps];
			gpsf=[fmg(1:k-1,[2 3]); gpsf];
			gmask=[ones(length([1:k-1]),1)
			       gmask];
		end
	end
	while(k <= lnew),
		lgp=size(gps,1);
% matlab matrix is empty so all data is new and fill it in
		if lgp==0,
			gps=mg;
			gpsf=fmg(:,[2 3]);
			k=lnew+1;
			gmask=[ones(lnew,1)];
		else,
% matlab matix is not empty and fill it in correct spots
			[index]=find(gps(:,1) < mg(k,1));
			idmz=max(index);
% check if we are at the end of the array if so fill to the
% end
			if idmz==lgp,
				gps=[gps
				     mg([k:lnew],:)];
				gpsf=[gpsf
				      fmg([k:lnew],[2:3])];
				gmask=[gmask
				       ones(length([k:lnew]),1)];
% we just have one point in the middle to put in
			else,
				gps=[gps([1:idmz],:)
		                     mg(k,:)
                                     gps([idmz+1:lgp],:)];
				gpsf=[gpsf([1:idmz],:)
				      fmg(k,[2:3])
				      gpsf([idmz+1:lgp],:)];
				gmask=[gmask([1:idmz])
 				       1
                                       gmask([idmz+1:lgp])];
			end
		end
		k=k+1;
	end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%>>>>****Ended here until figured out what to do with 
% GPS data (keep it in a file with ADCP or seperate
%%%%%>>>>
%% Sort out the ADCP data now.
% we now have the data.  Now we need to check times and sort times?
% since the scalar data comes from the same as the vector data we
% should process them at the same time here.  Both should have the
% same time stamps.
% ADCP array formats are given above.  We will have a Time vector and
% Depth vector with arrays of U,V, Error Vel, AGC
% Note m1 and m2 should be the same length.  They are both computed from the
% same files.  If for some wierd reason they are not the same we need to worry.
if ((isempty(m1) == 0)&(isempty(m2) == 0)),
	dt1=diff(m1(:,1));
	dt2=diff(m2(:,1));
%check for neg values to see if we have times out of order.
	[ineg1]=find(dt1 < 0);
	[ineg2]=find(dt2 < 0);
% we want to do all the time check stuff as vetor/matrix checks and
% try to avoid costly loops and redoing stuff that we don't need to
% do more than once.
	if isempty(ineg1)==0,
% we need to do a sort of the file
		[newt1 idex1]=sort(m1(:,1));
		m1=[newt1 m1(idex1,2:size(m1,2))];
	end
	if isempty(ineg2)==0,
		[newt2 idex2]=sort(m2(:,1));
		m2=[newt2 m2(idex2,2:size(m2,2))];
	end
% again below we have to worry about the year change.  This could
% be a recipe for disaster.  But try to get something that runs now
% and patch up the year bug later.
% data is now sorted.  Now we need to find which values are in the
% cleaned data file already and remove get send them to the nowhere.
% we only want to process what we have to.  If I can think up a vectorized
% method for doing this, then this should be done before the sort to make
% the sort faster (fewer points sort).
% find the minimum time in the currently captured files
	tmin1=min(m1(:,1));
	tmin2=min(m2(:,1));
	if isempty(adt)==0,
		tmax1=max(adt);
		cmin1=min(adt);
	else,
		tmax1=0;
		cmin1=370;
	end
	if isempty(adcpsc)==0,
		tmax2=max(adcpsc(:,1));
		cmin2=min(adcpsc(:,1));
	else,
		tmax2=0;
		cmin2=370;
	end
% Note: Only m1 is beyond tmax of file.
% Now a hard check to do.  Look for equal times and remove them
% from the m1 file because we have already dealt with them.
	lm1=size(m1(:,1),1);
	lm2=size(m2(:,1),1);
	im1=zeros(lm1,1);
	im2=zeros(lm2,1);
	if isempty(adt)==0,
		for i=1:lm1,
			if m1(i,1) > tmax1,
				i=lm1;
			end
			id=find(adt(:,1)==m1(i,1));
			if isempty(id)==0,
				im1(i)=1;
			end
		end
	end
	if isempty(adcpsc)==0,
		for i=1:lm2,
			if m2(i,1) > tmax2,
				i=lm2;
			end
			id=find(adcpsc(:,1)==m2(i,1));
			if isempty(id)==0,
				im2(i)=1;
			end
		end
	end
	m1(logical(im1),:)=[]; % Remove points we already have in the files.
	m2(logical(im2),:)=[];
end
if ((isempty(m1) == 0)&(isempty(m2) == 0)),
%
% Now we need to stuff our values in.  We also want to keep a "binary"
% mask of where the values are put in for computing our filters later.
% compare against current day -1 and see how much of old file
% we need to retain before stuffing in our values in a slow loop.
% 
	i=0; j=0;
% adt=acoustic doppler time (vector)
	ladt=size(adt,1);
	ladcpsc=size(adcpsc,1);
% even though we have seperate arrays for U,V,and such we only need one
% mask for point insertion since we insert for all arrays at the same point.
% we do however need a seperate mask for masking the insertion of good 
% adcp x-ducer temps.
% (we should only need one mask, theoretically the sc file and the uv file
% should be the same sizes, use 2 for now and change later?)
	bmask1=zeros(ladt,1);
	bmask2=zeros(ladcpsc,1);
	lnew1=size(m1,1);
	lnew2=size(m2,1);
	k=1;
% Format of m1 data file is:
% time,u(1)-->u(n),9999,v(1)-->v(n),9999,e(1)-->e(n),9999,cm1(1)-->cm1(n),9999,
% cm2(1)->cm2(n),9999,cm3(1)->cm3(n),9999,...9999,agc1(1)->agc1(n),9999,
% agc2(1)->agc2(n),9999,...,9999,pg3(1)->pg3(n),9999,pg4(1)->pg4(n).
% Format of m2 data file is:
% time,head,stdhead,pitch,stdpitch,roll,stdroll,xtemp,nbin,nping,
% bin_length,blank_length,depth,ver,rev,sspd,tpmin,tpsec,tphun,
% ex,ez,wa
% We only want to flag heading, stdhead, pitch,stdpitch, roll, stdroll, and xtemp.
%
% before we start the while loop we need to know which array will be
% bigger n or the old arrays.
	nbin=m2(1,9);
	narr=size(u,2);
% This is the usual case arrays are the same size and we don't
% have to pad either to get them to fit.
	if nbin==narr,
		lpad=0;
	end
% This case old array had more bins so we need to pad new data to
% fill in old array size.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NOTE: At this point we need to recompute the depth vector
% if nbin < narr we don't need a new vector, BUT it could be
% possible that the deployment has changed the depth at which
% the ADCP is at.  Which will change all the bins.
% if nbin > narr: we need to add to the depth vector.  Also
% the depth of the ADCP could change and cause all types of
% problems.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NOTE: Proper depths can be gotten for all bins by computing 
% using the values in the scalar ADCP data file.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	if nbin < narr,
		lpad=narr-nbin;
	end
% This case old arrays are too small and we need to expand to fill
% in the data.
	if nbin > narr,
		lpad=nbin-narr;
		lx=size(u,1);
		pad=ones(lx,lpad).*NaN;
		pad3=ones(lx,lpad,4).*NaN;
		u=[u pad];
		v=[v pad];
		e=[e pad];
		g=[g pad3];
		cm=[cm pad3];
		pg3=[pg3 pad];
		pg4=[pg4 pad];
		uf=[uf pad];
		vf=[vf pad];
		lpad=0;
	end
%*************************************************************
%************************************************************
% Warning!! We can't just stuff velocities in.  The stored
% velocities will be corrected for mooring motion AND rotated
% to degrees true (take into account magnetic declination)
%************************************************************
%************************************************************
% Unfortuneatly, we can't quite use multidimensional arrays like
% the other arrays.  If we append to 1 "page" the others get a
% row of zeros which we don't want.  So we pull out the pages
% temporarily as we append data and then recreate our multi-dim
% array.
	if(isempty(cm)==0),
		c1=squeeze(cm(:,:,1));
		c2=squeeze(cm(:,:,2));
		c3=squeeze(cm(:,:,3));
		c4=squeeze(cm(:,:,4));
		g1=squeeze(g(:,:,1));
		g2=squeeze(g(:,:,2));
		g3=squeeze(g(:,:,3));
		g4=squeeze(g(:,:,4));
	else,
		c1=[]; c2=[]; c3=[]; c4=[];
		g1=[]; g2=[]; g3=[]; g4=[];
	end
	clear cm; clear g;
	if(cmin1==370),
		adt=m1(:,1);
		k=lnew1+1;
		bmask1=[bmask1
		       	ones(lnew1,1)];
		u=m1(:,[2:(nbin+1)]).*100;
		v=m1(:,[(nbin+3):(2*nbin+2)]).*100;
		e=m1(:,[(2*nbin+4):(3*nbin+3)]).*100;
		c1=m1(:,[(3*nbin+5):(4*nbin+4)]);
		c2=m1(:,[(4*nbin+6):(5*nbin+5)]);
		c3=m1(:,[(5*nbin+7):(6*nbin+6)]);
		c4=m1(:,[(6*nbin+8):(7*nbin+7)]);
		g1=m1(:,[(7*nbin+9):(8*nbin+8)]);
		g2=m1(:,[(8*nbin+10):(9*nbin+9)]);
		g3=m1(:,[(9*nbin+11):(10*nbin+10)]);
		g4=m1(:,[(10*nbin+12):(11*nbin+11)]);
		pg3=m1(:,[(11*nbin+13):(12*nbin+12)]);
		pg4=m1(:,[(12*nbin+14):(13*nbin+13)]);
		uf=zeros(lnew1,nbin);
		vf=zeros(lnew1,nbin);
	else,
		if(m1(k,1) < cmin1),
			while((m1(k,1) < cmin1)&(k < size(m1,1))),
				k=k+1;
			end
			adt=[m1(1:k-1,1); adt];
			bmask1=[ones(length([1:k-1]),1)
		       	        bmask1];
			u=[m1(1:k-1,[2:(nbin+1)]).*100; u];
			v=[m1(1:k-1,[(nbin+3):(2*nbin+2)]).*100; v];
			e=[m1(1:k-1,[(2*nbin+4):(3*nbin+3)]).*100; e];
			c1=[m1(1:k-1,[(3*nbin+5):(4*nbin+4)]); c1];
			c2=[m1(1:k-1,[(4*nbin+6):(5*nbin+5)]); c2];
			c3=[m1(1:k-1,[(5*nbin+7):(6*nbin+6)]); c3];
			c4=[m1(1:k-1,[(6*nbin+8):(7*nbin+7)]); c4];
			g1=[m1(1:k-1,[(7*nbin+9):(8*nbin+8)]); g1];
			g2=[m1(1:k-1,[(8*nbin+10):(9*nbin+9)]);g2];
			g3=[m1(1:k-1,[(9*nbin+11):(10*nbin+10)]);g3];
			g4=[m1(1:k-1,[(10*nbin+12):(11*nbin+11)]);g4];
			pg3=[m1(1:k-1,[(11*nbin+13):(12*nbin+12)]); pg3];
			pg4=[m1(1:k-1,[(12*nbin+14):(13*nbin+13)]); pg4];
			uf=[zeros(length([1:k-1]),nbin);uf];
			vf=[zeros(length([1:k-1]),nbin);vf];
		end
	end	
	while(k <= lnew1),
		ladt=size(adt,1);
		if ladt==0,
% for this case all arrays are empty and we just fill in
			adt=m1(:,1);
			k=lnew1+1;
			bmask1=[bmask1
		       	       ones(lnew1,1)];
			u=m1(:,[2:(nbin+1)]).*100;
			v=m1(:,[(nbin+3):(2*nbin+2)]).*100;
			e=m1(:,[(2*nbin+4):(3*nbin+3)]).*100;
			c1=m1(:,[(3*nbin+5):(4*nbin+4)]);
			c2=m1(:,[(4*nbin+6):(5*nbin+5)]);
			c3=m1(:,[(5*nbin+7):(6*nbin+6)]);
			c4=m1(:,[(6*nbin+8):(7*nbin+7)]);
			g1=m1(:,[(7*nbin+9):(8*nbin+8)]);
			g2=m1(:,[(8*nbin+10):(9*nbin+9)]);
			g3=m1(:,[(9*nbin+11):(10*nbin+10)]);
			g4=m1(:,[(10*nbin+12):(11*nbin+11)]);
			pg3=m1(:,[(11*nbin+13):(12*nbin+12)]);
			pg4=m1(:,[(12*nbin+14):(13*nbin+13)]);
			uf=zeros(lnew1,nbin);
			vf=zeros(lnew1,nbin);
		else,
% here the arrays are not empty and we have to be carefull.
			[index]=find(adt(:,1) < m1(k,1));
			idmz=max(index);
			if idmz==ladt,
				adt=[adt
                     	     	     m1([k:lnew1],1)];
				bmask1=[bmask1
			               ones(length(k:lnew1),1)];
				pad=ones(length(k:lnew1),lpad).*NaN;
				uf=[uf
				    zeros(length(k:lnew1),nbin) pad];
				vf=[vf
				    zeros(length(k:lnew1),nbin) pad];
				u=[u
                                   m1([k:lnew1],[2:(nbin+1)]).*100 pad];
				v=[v
				   m1([k:lnew1],[(nbin+3):(2*nbin+2)]).*100 pad];
				e=[e
                                   m1([k:lnew1],[(2*nbin+4):(3*nbin+3)]).*100 pad];
				c1=[c1
				   m1([k:lnew1],[(3*nbin+5):(4*nbin+4)]) pad];
				c2=[c2
				   m1([k:lnew1],[(4*nbin+6):(5*nbin+5)]) pad];
				c3=[c3
				   m1([k:lnew1],[(5*nbin+7):(6*nbin+6)]) pad];
				c4=[c4
				   m1([k:lnew1],[(6*nbin+8):(7*nbin+7)]) pad];
				g1=[g1
				   m1([k:lnew1],[(7*nbin+9):(8*nbin+8)]) pad];
				g2=[g2
				   m1([k:lnew1],[(8*nbin+10):(9*nbin+9)]) pad];
				g3=[g3
				   m1([k:lnew1],[(9*nbin+11):(10*nbin+10)]) pad];
				g4=[g4
				   m1([k:lnew1],[(10*nbin+12):(11*nbin+11)]) pad];
				pg3=[pg3
				   m1([k:lnew1],[(11*nbin+13):(12*nbin+12)]) pad];
				pg4=[pg4
				   m1([k:lnew1],[(12*nbin+14):(13*nbin+13)]) pad];
			 	k=lnew1;
			else
				adt=[adt([1:idmz],:)
		     	     	     m1(k,1)
		     	             adt([idmz+1:ladt],:)];
				bmask1=[bmask1([1:idmz])
				       1
			               bmask1([idmz+1:ladt])];
				pad=ones(1,lpad).*NaN;
				uf=[uf([1:idmz],:)
				    zeros(1,nbin) pad
				    uf([idmz+1:ladt],:)];
				vf=[vf([1:idmz],:)
				    zeros(1,nbin) pad
				    vf([idmz+1:ladt],:)];
				u=[u([1:idmz],:)
				   m1(k,[2:(nbin+1)]).*100 pad
			           u([idmz+1:ladt],:)];
				v=[v([1:idmz],:)
				   m1(k,[(nbin+3):(2*nbin+2)]).*100 pad
			           v([idmz+1:ladt],:)];
				e=[e([1:idmz],:)
				   m1(k,[(2*nbin+4):(3*nbin+3)]).*100 pad
				   e([idmz+1:ladt],:)];
				c1=[c1([1:idmz],:)
				   m1(k,[(3*nbin+5):(4*nbin+4)]) pad
				   c1([idmz+1:ladt],:)];
				c2=[c2([1:idmz],:)
				   m1(k,[(4*nbin+6):(5*nbin+5)]) pad
				   c2([idmz+1:ladt],:)];
				c3=[c3([1:idmz],:)
				   m1(k,[(5*nbin+7):(6*nbin+6)]) pad
				   c3([idmz+1:ladt],:)];
				c4=[c4([1:idmz],:)
				   m1(k,[(6*nbin+8):(7*nbin+7)]) pad
				   c4([idmz+1:ladt],:)];
				g1=[g1([1:idmz],:)
				   m1(k,[(7*nbin+9):(8*nbin+8)]) pad
				   g1([idmz+1:ladt],:)];
				g2=[g2([1:idmz],:)
				   m1(k,[(8*nbin+10):(9*nbin+9)]) pad
				   g2([idmz+1:ladt],:)];
				g3=[g3([1:idmz],:)
				   m1(k,[(9*nbin+11):(10*nbin+10)]) pad
				   g3([idmz+1:ladt],:)];
				g4=[g4([1:idmz],:)
				   m1(k,[(10*nbin+12):(11*nbin+11)]) pad
				   g4([idmz+1:ladt],:)];
				pg3=[pg3([1:idmz],:)
				   m1(k,[(11*nbin+13):(12*nbin+12)]) pad
				   pg3([idmz+1:ladt],:)];
				pg4=[pg4([1:idmz],:)
				   m1(k,[(12*nbin+14):(13*nbin+13)]) pad
				   pg4([idmz+1:ladt],:)];
			end
		end
		k=k+1;
	end
	cm(:,:,1)=c1;
	cm(:,:,2)=c2;
	cm(:,:,3)=c3;
	cm(:,:,4)=c4;
	clear c1; clear c2; clear c3; clear c4;
	g(:,:,1)=g1;
	g(:,:,2)=g2;
	g(:,:,3)=g3;
	g(:,:,4)=g4;
	clear g1; clear g2; clear g3; clear g4;
	k=1;
	if(cmin2==370),
		adcpsc=m2;
		k=lnew2+1;
		bmask2=[bmask2
			ones(lnew2,1)];
		fsc=zeros(lnew2,size(adcpsc,2)-1);
	else
		if(m2(k,1) < cmin2),
			while((m2(k,1) < cmin2)&(k < size(m2,1))),
				k=k+1;
			end
			adcpsc=[m2([1:k-1],:); adcpsc];
			fsc=[zeros(length([1:k-1]),size(adcpsc,2)-1); fsc];
			bmask2=[ones(length([1:k-1]),1); bmask2];
		end
	end
	while(k <= lnew2),
		ladcpsc=size(adcpsc,1);
		if ladcpsc==0,
			adcpsc=m2;
			k=lnew2+1;
			bmask2=[bmask2
		       	       ones(lnew2,1)];
			fsc=zeros(lnew2,size(adcpsc,2)-1);
		else,
			[index]=find(adcpsc(:,1) < m2(k,1));
			idmz=max(index);
			if idmz==ladcpsc,
				adcpsc=[adcpsc
                     	     	        m2([k:lnew2],:)];
				bmask2=[bmask2
			               ones(length(k:lnew2),1)];
				fsc=[fsc
				     zeros(length(k:lnew2),size(m2,2)-1)];
			 	k=lnew2;
			else
				adcpsc=[adcpsc([1:idmz],:)
		     	     	     	m2(k,:)
		     	             	adcpsc([idmz+1:ladcpsc],:)];
				bmask2=[bmask2([1:idmz])
				       1
			               bmask2([idmz+1:ladcpsc])];
				fsc=[fsc([1:idmz],:)
				     zeros(1,size(m2,2)-1)
				     fsc([idmz+1:ladcpsc],:)];
			end
		end
		k=k+1;
	end
% Do a simple test to put marked bad velocities (by the adcp software)
% to have flags of 3 (no data).
	for yi=1:nbin,
		yu=u(:,yi)<=-3276.8;
		yv=v(:,yi)<=-3276.8;
		uf(yu,yi)=3;
		vf(yv,yi)=3;
	end
% If the above works, bmask should be a file with 1's where
% we have inserted data and should be the same length as
% the array met with the inserted data.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now we can apply our gradient filter to the data (?)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For now, no filter on adcp velocities........
%====> put filter here <====%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Need to save the appropriate data here!
% we need to be aware of year and month breaks
% in our data if we want to save it correctly.
	tn=adt;
	un=u;
	vn=v;
	en=e;
	cmn=cm;
	gn=g;
	pg3n=pg3;
	pg4n=pg4;
	fu=uf;
	fv=vf;
	gpsn=gps;
	fgps=gpsf;
	adcpscn=adcpsc;
	adcpfsc=fsc;
	eval(['save ',pathm1,fnm1,' tn un vn en cmn gn pg3n pg4n fu fv gpsn fgps']);
	eval(['save ',pathm1,fnm2,' adcpscn adcpfsc']);
end % end of im m1 has data
%quit
