sortomp2.m

用MATLAB进行海洋水团的模拟程序

%  SORTOMP2.M: sorts water type results and residuals,
%  latitude and longitude into increasing pressure and
%  writes the data as sections into files.
%
%  NOTES: 1. 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.
%         3. This routine assumes that you have a directory/folder
%            "output" in the current work directory and places all
%            files into that directory/folder.
%
%---------------------------------------------
% 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]
%--------------------------------------------
%
%  This routine produces ASCII output in the following format:
%
%  <number of depth levels> <number of stations>
%  <maxvar> <minvar>
%  <depth level 1> <depth level 2> <depth level 3> ... etc.
%  <station location 1> <station location 2> <station location 3> ... etc.
%  <data station 1 depth 1> <data station 2 depth 1> <data station 3 depth 1> ... etc.
%  <data station 1 depth 2> <data station 2 depth 2> <data station 3 depth 2> ... etc.
%  <data station 1 depth 3> <data station 2 depth 3> <data station 3 depth 3> ... etc.
%  ... etc.
%
%  This format is read automatically by Transform (Spyglass/Fortner) and Noesys (Fortner).
%
%  Notes: 1. "Station location" means "Distance from first station" in km.
%         2. maxvar and minvar are the maximum and minimum value for the data. It is set
%            to zero, which forces Transform to calculate these values. This is done because
%            in most cases max(max(data)) and min(min(data)) turns out to be NaN.
%            If you are sure that there are no NaNs in your data you can include the
%            calculation of maxvar and minvar from the data in the code.
%

%name variables
varname = [''];
varname(1,1:11) = ['   latitude'];
varname(2,1:11) = ['  longitude'];
varname(3,1:11) = ['   pressure'];
varname(4,1:11) = ['   salinity'];
varname(5,1:11) = [' pot. temp.'];
varname(6,1:11) = ['     oxygen'];
varname(7,1:11) = ['  phosphate'];
varname(8,1:11) = ['    nitrate'];
varname(9,1:11) = ['   silicate'];
varname(10,1:11) = [' pot. vort.'];
varname(11,1:11) = ['mass consn.'];

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


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

clear ptsort
maxvar = 0;
minvar = 0;

%warning to check existing files
m = 1;
disp('  ')
disp('WARNING! This program will create the following files:')
for i = 4:10
	if esx(i) == 1 disp([filname(i,:) '      (for ' varname(i,1:11) ')']); m = m+1; end
end
disp('  ')
for i = 4:10
	if esx(i) == 1 disp([errname(i,:) '      (for ' varname(i,1:11) ' residual)']); end
end
disp([errname(11,:) '      (for ' varname(11,1:11) ' residual)']);
disp('  ')

for i = 1:nr_of_wm
	disp([wmname(i,1:7) '       (for ' tit_index(5*(i-1)+1:5*(i-1)+5) ')']);
end
want = 'n';
drswitch('It will place the files into');
incontrol = input('Do you want to continue (y/n)? [n] ','s');
if incontrol== 'y' want = 'y'; end

if want == 'n' return; end

%main block, only run if want == 'y'
%find start and end of stations
statind=[find(diff(press)<0) length(press)]; 

%calculate distance between stations
[dist,phaseangle] = sw_dist(lat,long,'km');
cumdist=[0 cumsum(dist)];
cumdist = cumdist(statind);
statind = [0 statind];
lcd = length(cumdist);

%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 = 4:10
if esx(nn) == 1
	disp(['working on ' varname(nn,1:11) ' 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 = filname(nn,1:8);
 	sout = sprintf('save %s llv lcd maxvar minvar levels cumdist ptsort -ascii',fout);
 	eval(sout);
	disp(['file ' filname(nn,1:8) ' created and saved'])
	disp('  ')
	ptsort = ptsort';
end %esx
end %nn

for nn = 1:7
if esx(nn+3) == 1
	disp(['working on ' varname(nn+3,1:11) ' 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('save %s llv lcd maxvar minvar levels cumdist ptsort -ascii',fout);
	eval(sout);
	disp(['file ' errname(nn+3,1:8) ' created and saved'])
	disp('  ')
	ptsort = ptsort';
end %esx
end %nn

disp(['working on ' varname(11,1:11) ' 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('save %s llv lcd maxvar minvar levels cumdist ptsort -ascii',fout);
eval(sout);
disp(['file ' errname(11,1:8) ' created and saved'])
disp('  ')
ptsort = ptsort';

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('save %s llv lcd maxvar minvar levels cumdist ptsort -ascii',fout);
	eval(sout);
	disp(['file ' wmname(nn,1:7) ' created and saved'])
	disp('  ')
	ptsort = ptsort';
end %esx
end %nn
clear ptsort levels cumdist dist phaseangle llv lcd maxvar minvar nn fout sout want
clear ans lat long oxy ph ni si press ptemp pvort temp sal tit_index statind