%add_bias
%routine to add a bias to salt over specified pres bounds

stas=input('Which station do you want to process?  ');

 [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]);


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


eval(['load ',matfmtfile]);

thetave=sw_ptmp(s1,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,'k');

jstart=input(' Enter the first scan number of the bad C interval ');
jstop=input('   Now enter the final scan no. of the bad C interval ');

bias=input('  Input the amount to bias salt by ');

for n=jstart:jstop
s_ave(n)=s_ave(n)+bias;
cave(n)=sw_cndr(s_ave(n),tave(n),pave(n))*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']);
         



