statns2.m

% STATNS2.M: sorts water type results from omp2.m and residuals,
%  latitude and longitude into increasing pressure and
%  produces new arrays which can be studied as individual stations.
%
%  NOTES:    This routine reads output from a file. It does
%            not require output from omp2.m in the workspace
%            but takes output from omp2.m from the file which was
%            generated by that program. It clears the workspace
%            before it starts.
%
%---------------------------------------------
% This program is part of the OMP package from:
% Institut fuer Meereskunde
% J. Karstensen 
% Troplowitzstr. 7
% 22529 Hamburg
% Germany
%
%  It was written by Matthias Tomczak.
%
% BUGS: [e-mail]
%--------------------------------------------
%

clear all

%name variables
varname = [''];
varname(1,1:8) = ['latitude'];
varname(2,1:8) = ['longitud'];
varname(3,1:8) = ['pressure'];
varname(4,1:8) = ['salinity'];
varname(5,1:8) = ['potltemp'];
varname(6,1:8) = ['dvoxygen'];
varname(7,1:8) = ['phosphat'];
varname(8,1:8) = ['dvnitrat'];
varname(9,1:8) = ['silicate'];
varname(10,1:8) = ['potlvort'];
varname(11,1:8) = ['masscons'];

filname = [''];
filname(1,1:8) = ['        '];
filname(2,1:8) = ['        '];
filname(3,1:8) = ['        '];
filname(4,1:8) = ['filsalin'];
filname(5,1:8) = ['filptemp'];
filname(6,1:8) = ['filoxyge'];
filname(7,1:8) = ['filphosp'];
filname(8,1:8) = ['filnitra'];
filname(9,1:8) = ['filsilic'];
filname(10,1:8) = ['filpvort'];
filname(11,1:8) = ['filtempe'];

errname(1,1:8) = ['        '];
errname(2,1:8) = ['        '];
errname(3,1:8) = ['        '];
errname(4,1:8) = ['errsalin'];
errname(5,1:8) = ['errptemp'];
errname(6,1:8) = ['erroxyge'];
errname(7,1:8) = ['errphosp'];
errname(8,1:8) = ['errnitra'];
errname(9,1:8) = ['errsilic'];
errname(10,1:8) = ['errpvort'];
errname(11,1:8) = ['errmcons'];

%switch to output directory/folder for input and output
drswitch('This program reads data from');

%define input file
incontrol = input('Which file for data input?  [result]  ','s');
if length(incontrol) > 0
	dataset = incontrol;
else
	dataset = 'result';
end
eval(['load ' dataset]);

m = 1;
for i = 4:10
	if esx(i) == 1 m = m+1; end
end

for i = 1:nr_of_wm
	wmname(i,1:7) = ['wamass' num2str(i)];
end

%find start and end of stations
if size(press,1)>1
 press=press';
end
 
statind=[find(diff(press)<0) length(press)]; 

%define interpolation levels
incontrol = input('Which interpolation file do you want to use?  [testlevels]  ','s');
disp('  ')
if length(incontrol) > 0
	dataset = incontrol;
else
	dataset = 'testlevels';
end
eval(['load ' dataset]);
llv = length(levels);

%check for duplicates and separate them by 0.5 m.
for i=2:length(press)
	if press(i-1) == press(i) press(i-1) = press(i-1) - 0.5; end
end

%initialize data and result matrices
ptsort(length(statind)-1,length(levels)) = NaN;
for i=1:length(statind)-1
	for j=1:length(levels)
		ptsort(i,j) = NaN;
	end
end

%sort data into station format
for nn = 1:10
if esx(nn) == 1
	disp(['working on ' varname(nn,1:8) ' extraction'])
	for i=2:length(statind)
		j = 1;
		while press(statind(i-1)+1)>levels(j) & j<length(levels)
			j = j+1;
		end
		k = length(levels);
		while press(statind(i))<levels(k) & k>1
			k = k-1;
		end
 		sortwrt %stores data in ptsort; see sortwrt.m
	end
 	ptsort = ptsort';
	fout = varname(nn,1:8);
 	sout = sprintf(' %s = ptsort;',fout);
 	eval(sout);
	ptsort = ptsort';
end %esx
end %nn
disp('  ')

for nn = 1:7
if esx(nn+3) == 1
	disp(['working on ' varname(nn+3,1:8) ' residual extraction'])
	for i=2:length(statind)
		j = 1;
		while press(statind(i-1)+1)>levels(j) & j<length(levels)
			j = j+1;
		end
		k = length(levels);
		while press(statind(i))<levels(k) & k>1
			k = k-1;
		end
		if j < k
			ptsort(i-1,j:k) = interp1(press(statind(i-1)+1:statind(i)), ...
				err(nn,statind(i-1)+1:statind(i)),levels(j:k))';
		end
	end
	ptsort = ptsort';
	fout = errname(nn+3,1:8);
 	sout = sprintf(' %s = ptsort;',fout);
	eval(sout);
	ptsort = ptsort';
end %esx
end %nn

disp(['working on ' varname(11,1:8) ' residual extraction'])
for i=2:length(statind)
	j = 1;
	while press(statind(i-1)+1)>levels(j) & j<length(levels)
		j = j+1;
	end
k = length(levels);
	while press(statind(i))<levels(k) & k>1
		k = k-1;
	end
	if j < k
		ptsort(i-1,j:k) = interp1(press(statind(i-1)+1:statind(i)), ...
			err(m,statind(i-1)+1:statind(i)),levels(j:k))';
	end
end
ptsort = 100 * ptsort';
fout = errname(11,1:8);
sout = sprintf(' %s = ptsort;',fout);
eval(sout);
ptsort = ptsort';
disp('  ')
for nn = 1:nr_of_wm
	disp(['working on ' tit_index(5*(nn-1)+1:5*(nn-1)+5) ' extraction'])
	for i=2:length(statind)
		j = 1;
		while press(statind(i-1)+1)>levels(j) & j<length(levels)
			j = j+1;
		end
		k = length(levels);
		while press(statind(i))<levels(k) & k>1
			k = k-1;
		end
		if j < k
			ptsort(i-1,j:k) = interp1(press(statind(i-1)+1:statind(i)), ...
				A(nn,statind(i-1)+1:statind(i)),levels(j:k))';
		end
	end
	ptsort = 100 * ptsort';
	fout = wmname(nn,1:7);
 	sout = sprintf(' %s = ptsort;',fout);
	eval(sout);
	ptsort = ptsort';
end %esx
end %nn
clear ptsort nn fout sout want

disp('  ')
disp('Extraction finished.')
disp('Run displvar.m to see individual stations.')

clear ans lat long oxy ph ni si press ptemp pvort temp sal tit_index statind