%irange = [1:250 252:2:500 505:5:1000 1005:20:ndepth];
%irange = [2500:2750];
%irange = [ 450:5:1000 1005:20:ndepth];% for LVS
irange = [240:2:650];

%jrange =  (nsta-15):nsta;
jrange = 1:nsta; % default, for all stations
%jrange = [76:80 82:2:112] - 75;
%jrange = [1 10 18 30 31];   % LVS

nj = length(jrange);
pp = p(irange, ones(nj,1));

subplot(2, 1, 1)
plot(S(irange,jrange), theta(S(irange,jrange),T(irange,jrange),pp, 0))
title(sprintf('I8S, stations %.0f to %.0f',jsta(jrange(1)),...
                                             jsta(jrange(nj))))
xlabel('S')
ylabel('theta') 
% text(0, 0, pctime, 'sc');    %print PC (current) time on the plot

subplot(2, 1, 2)
offset = ones(size(irange')) * (0:(nj-1)) * 0.02;
plot(S(irange,jrange)+offset, theta(S(irange,jrange),T(irange,jrange),pp, 0))
%plot(S(irange,jrange)+offset,T(irange,jrange))
ylabel('theta')
xlabel('S, successive offsets of 0.002 ppt')
title(sprintf('(Latitudes %.1f to %.1f)', y(jrange(1)),...
                                             y(jrange(nj))))
%   text(0, 0, pctime, 'sc');    %print PC (current) time on the plot

orient tall
print ts.ps



