function varargout=trsp2(year_base, x_var, gname, cruise_id, plot_title)
% ------------------------------------------------------------------
% function []=trsp2(year_base, x_var, gname, cruise_id, plot_title)
%
% Load transport data from `adcpsect' and plot vertically integrated
% and horizontally cumulated transport.
%
% PARAMETERS:
%
%     year_base
%        4 digit year base, NOT 1900-base.
%
%     x_var
%        String specifying the x axis variable. Possible values are:
%
%        'lon_180' : longitude [-180.0, 180.0)
%        'lon_360' : longitude [0.0, 360.0)
%        'lat'     : latitude
%        'year_day': decimal day
%     
%     gname
%        Generic name of mat files from `adcpsect'. This is the same as
%        the `output' parameter from the `adcpsect'¨configuration file.
%
%     cruise_id
%        String containing the cruise identification code.
%
%     plot_title
%        String containing the title for the plot.
%
%     pos
%        String for positions. It can take the following values:
%
%          'lon_180' : convert longitude to [-180, 180)
%          'lon_360' : convert longitude to [0, 360)
%
%
% The plot is printed in the file `trsp.ps'.
% ------------------------------------------------------------------
%Time:<Last updated on 98/11/30 at 14:23:18 by jaccard@pc34.gfi.uib.no>
%File:</home1/jaccard/data1/codas3/gfi/matlab/matlab50/trsp2.m>

% ------------------------------------------------------------------
% Initialization:
% ------------------------------------------------------------------

max_sig = 5.0;

% ------------------------------------------------------------------
% Load data:
% ------------------------------------------------------------------
file_tr = strcat(gname, '_tr');
file_uv = strcat(gname, '_uv');
file_xy = strcat(gname, '_xy');
ps_file = [ gname '.ps' ];
load(file_tr);
load(file_xy);
	
% ------------------------------------------------------------------
% Calculate the transport at each depth normal to track:
%
% This is given by (udy - vdx), with positive values 
% being transport toward starboard.
%
% This assumes that that the `transport' matrix contains
% elements (udy) and (-vdx). Otherwise, transport will be
% positive along path.
% ------------------------------------------------------------------
imax = size(transport,2);
kmax = size(transport,1);
zmax= size(z,1);

iu = 1:2:imax;
iv = 2:2:imax;
ik = 1:1:kmax;

Mxz = transport(ik,iu) + transport(ik,iv);

%------------------------------------------------------------------
% Smooth the transport data:
%
%    If too few points have been used to calculate the transport,
%    some values might be much too high. In order too avoid editing
%    the transport file, this part tries to edit out bad values, in
%    a similar way as `btcaluv'.
%------------------------------------------------------------------
imax = size(Mxz,2);
kmax = size(Mxz,1);
good_mask = ones(kmax, imax);

%--- Flag bad data, and fill gaps along x-axis:
for k=1:kmax,
	[Mxz(k,:), good_mask(k,:)] = flag_bad(Mxz(k,:), max_sig);
	Mxz(k,:) = fillgaps(Mxz(k,:), 0, 1);
	%--- Complete rows of NaN are set zero:
	ib = find(isnan(Mxz(k,:)));
	if length(ib) > 0,
		Mxz(k,ib) = 0;
	end
end

%------------------------------------------------------------------
% Integrate transport over depth
%------------------------------------------------------------------
M = sum(Mxz,1);

% ------------------------------------------------------------------
% Convert positions:
% ------------------------------------------------------------------
switch x_var
	case 'lon_180',
		ii = find(xyt(1,:) >= 180.0);
		if length(ii),
			xyt(1,ii) = xyt(1,ii) - 360.0;
		end
		x = 1;
	case 'lon_360',
		ii = find(xyt(1,:) < 0.0);
		if length(ii),
			xyt(1,ii) = xyt(1,ii) + 360.0;
		end
		x = 1;
	case 'lat',
		x = 2;
	case 'year_day',
		x = 3;
	otherwise,
		error([ 'Unknown x variable string ' pos ]);
end
x_var = x;
xmax  = size(xyt(x_var,:),2);

% ------------------------------------------------------------------
% Plot transport:
% ------------------------------------------------------------------
figure;
orient landscape;
hold on;
subplot('position', [0.1, 0.1, 0.8, 0.7]);
% ---------------------------------
plot(xyt(x_var,:), M(:), '-');
% ---------------------------------
if x_var == 1,
	xlabel('longitude');
elseif x_var == 2,
   xlabel('latitude');
elseif x_var == 3,
	str0 = num2str(year_base);
	str1 = cat(2, str0, ' decimal day');
	xlabel(str0);
end
% ---------------------------------
ylabel('cumulated transport (m^3/s)');

% ------------------------------------------------------------------
% Put crosses for bad points:
% ------------------------------------------------------------------
gm_ymax = max(M);
gm_ymin = min(M);
gm_zmax = max(z);
gm_zmin = min(z);
gm_m    = (gm_ymax - gm_ymin) / (gm_zmin - gm_zmax);
gm_h    = (gm_ymax*gm_zmax - gm_ymin*gm_zmin) / (gm_zmax - gm_zmin);
gm_dz   = (gm_zmax - gm_zmin)/(kmax - 1);
gm_y    = gm_m*(gm_zmin:gm_dz:gm_zmax) + gm_h;
gm_x    = xyt(x_var,:);
for i=1:xmax
	for j=1:kmax
		if good_mask(j,i) == 0,
			plot(gm_x(i), gm_y(j), 'x', 'MarkerSize', 3);
		end
	end
end
grid on

% ------------------------------------------------------------------
% Write additionnal stuff:
% ------------------------------------------------------------------
note_str = [ - 'cruise_id: ' cruise_id ' - dbname: ' dbname ...
      ' - regrid: ' regrid_method ' ' regrid_var ' - trsp: ' direction];
subtitle = [ to_date(year_base, xyt(3,1)) ' to ' ...
      to_date(year_base, xyt(3,size(xyt,2))) ];
if strcmp(direction, 'ALONG'),
	info = 'positive along track';
else
	info = 'positive across track, starbord';
end
info = [ info ' - ' ...
      z_label ' from ' num2str(z(1)) ' to ' num2str(z(zmax))];
inst = sprintf('(%08.4f %07.4f) to (%08.4f %07.4f)', ...
    xyt(1,1), xyt(2,1), xyt(1,size(xyt,2)), xyt(2,size(xyt,2)));
plotstuf('-title', plot_title, '-subtitle', subtitle, '-note', note_str, ...
    '-info', info, '-institute', inst, '-format', 'landscape');


eval([ 'print -dps ' ps_file ]);


%% 1999/03/02 Jaccard Pierre 
%%            Added possibility to return calculated transport to calling
%%            process. 

if nargout > 0,
  varargout{1} = M;
end
if nargout > 1,
  varargout{2} = xyt;
end



