% This script get profiles and save each one in a different file
% The name of the file will contain all the information of the 
% data, position, depth, code number for the run and time (the
% code will ask you for a reference time).
%
%  The information for extractin the profiles have to be 
%  in the pro.dat file, the first column will be the time, 
%  the second the x position and the third the y position. 
%  All this information is in grid-coordinates.
%  You need to edit this file and modificy the filein in case 
%  it is necesary.
%
%  Pedro Velez (UIB), July 1999
%
%  Note: This is not a function, just a script

filein='pe_out.nc';
filepro='pro.dat' ;

disp('>>>>> By default the model file is pe_model.nc')
disp('>>>>> and pro.dat contains the information')
disp('>>>>> about the profiles you want to extract')
disp('>>>>> The following inputs have to be numbers')

modelrun=input('>>>>>Number of model run: ');
absyear=input('>>>>>Write the reference year:  ');
absmonth=input('>>>>>Write the reference month:  ');
absday=input('>>>>>Write the reference day:  ');
abshour=input('>>>>>Write the reference hour:  ');

load(filepro)

%Loop for the different profiles
for n=1:size(pro,1)

time=pro(n,1);
ix=pro(n,2);
iy=pro(n,3);


%Reading the data
x=getnc(filein,'tgrid3',[1 -1 1 1],[1 -1 1 1],1,1,2);
y=getnc(filein,'tgrid3',[-1 1 1 2],[-1 1 1 2],1,1,2);
depth=getnc(filein,'tgrid3',[iy ix -1 3],[iy ix -1 3],1,1,2);
temp=getnc(filein,'temp',[time iy ix -1],[time iy ix -1],1,1,2);
salt=getnc(filein,'salt',[time iy ix -1],[time iy ix -1],1,1,2);
souvel=svel80(salt,temp,-depth);



%Setting the lon and lat
lon=x(ix);
lat=y(iy);
glat= floor(abs(lat));
glon= floor(abs(lon));
mlat = floor(60*(abs(lat)-glat));
mlon = floor(60*(abs(lon)-glon));
slon=round((abs(lon)-glon)*3600.-(60.*mlon));
slat=round((abs(lat)-glat)*3600.-(60.*mlat));

if lat > 0 
ns='N';
else 
ns='S';
end
if lon > 0 
we='E';
else 
we='W';
end

strlat=strrep(sprintf('%s%2d%2d%2d',ns,glat,mlat,slat),' ','0');
strlon=strrep(sprintf('%s%2d%2d%2d',we,glon,mlon,slon),' ','0');


%Setting the absolute time
ntlast=getnc(filein,'nlast');
ntsout=getnc(filein,'ntsout');
dtuv=getnc(filein,'dtuv');
timemax=ntlast/ntsout+1;
timeinterval=(dtuv*ntsout)/3600. ;
hora=((time-1)*timeinterval)+abshour;
sfecha=datestr(datenum(absyear,absmonth,absday,hora,0,0),0);
dia=sfecha(1:2);
mes=sfecha(4:6);
hora=sfecha(13:14);
strfecha=sprintf('%s%s%s',dia,mes,hora);


fileout=strrep(sprintf('%2d_%s_%s_%s.asc',modelrun,strfecha,strlat,strlon),' ','0');

subplot(1,3,1) 
	plot(souvel,depth,'k');xlabel('Sound Velocity');ylabel('Depth');grid on
subplot(1,3,2) 
	plot(temp,depth,'k');xlabel('Temperature');ylabel('Depth');grid on;title(fileout,'interpreter','none')
subplot(1,3,3) 
	plot(salt,depth,'k');xlabel('Salinity');ylabel('Depth');grid on
drawnow


%Writing the data 
fid=fopen(fileout,'w');
for i=1:size(depth,1)
fprintf(fid,'%10.3f  %10.3f  %10.3f  %10.3f \n',depth(i),temp(i),salt(i),souvel(i));
end


end
