%fix_c_big
%routine to interpolate over large segments of bad conductivity data
%  uses neighboring profiles and theta-potc relationships

stas=input('Which station do you want to process?  ');
sta1=input('  Enter the first reference station ');
sta2=input('  Enter the second reference station ');

 [directory]=input('Enter the directory with the input grd mat file:  ','s');
 [outraw]=input('Now enter the directory where you want the new raw mat file to go:  ','s');


matfmtfile=[directory,'/grd',sprintf('%4.4i',stas),'.mat'];
newfile=[outraw,'/grd',sprintf('%4.4i',stas),'.mat'];

disp(['   working on file ',matfmtfile]);

matfmtfile1=[directory,'/grd',sprintf('%4.4i',sta1),'.mat'];
matfmtfile2=[directory,'/grd',sprintf('%4.4i',sta2),'.mat'];

eval(['load ',matfmtfile1]);
th1=thetave;
s1=s_ave;

eval(['load ',matfmtfile2]);
th2=thetave;
s2=s_ave;

eval(['load ',matfmtfile]);

s=(s1+s2)/2;
thetave=sw_ptmp(s,tave,pave,0);

disp([' ']);
disp([' Zoom in to find the scan numbers bracketing the spike ']);
disp([' ']);
figure(1)
clf
plot(s_ave)
grid
zoom on
hold on
title('Salinity')

figure(2)
clf
plot(s_ave,thetave);
grid
zoom on
hold on
xlabel('Salinity');
ylabel('Pot. Temperature');
plot(s1,th1,'r');
plot(s2,th2,'m');

jstart=input(' Enter the scan number just prior to the C spike ');
jstop=input('   Now enter the first good C scan no. after the spike ');

for j=start:stop
diff1=abs(th1-thetave(j));
diff2=abs(th2-thetave(j));

min1=min(diff1);
min2=min(diff2);

k1=find(diff1-min1);
k2=find(diff2-min2);

s_ave(j)=(s1(k1)+s2(k2))/2;
cave(j)=sw_cndr(s_ave(j),tave(j),pave(j))*sw_c3515;

end;


thetave=sw_ptmp(s_ave,tave,pave,0);
sigthave=sw_pden(s_ave,tave,pave,0)-1000.;


figure(1)
plot([jstart:jstop],s_ave(jstart:jstop),'r');
figure(2)
plot(s_ave(jstart:jstop),thetave(jstart:jstop),'r');

eval(['save ',newfile,' ascan1 file pstop thetave ascan2 pave s_ave ', ...
'uave cave pedate sigthave vave cscan1 pgrid startdaytime wave cscan2 ', ...
'psdate stopdaytime dpdtave pstart tave']);
         



