%* slice levitus data
%  draw filled contour transects through the world ocean

% levitus94_sal.mat		mean annual salinity (p.s.u.), 1 deg graticule, 
%                       Z depths
% levitus94_temp.mat	mean annual temperature (K above triple point), 
%                       1 deg graticule, Z depths
%
% Z (-ve):  0 10 20 30 50 75 100 125 150 200 250 300 400 500 600 700 800 900 
%			1000 1100 1200 1300 1400 1500 1750 2000 2500 3000 3500 4000 4500 
%			5000 5500

%* load the data
load levitus94_sal
load levitus94_temp
load coast


%* make some ploting variables
%  -ve is W lon, S lat
%  +ve is E lon, N lat

maplat=[-90:1:89];
maplon=[-180:1:179];


%* SET slice longitude
SL=-20;                       
slicelon=find(maplon==SL);      %get index number for selected longitude

%* set starting figure numbers, program uses FN, FN+1, and FN+2
FN=111;                         %starting figure number, change to save oldies
FNref=100;                      %reference map, don't change or overwrite


%* ranges for color scales
Csmin=min(min(min(sal)));							
Csmax=max(max(max(sal)));

Ctmin=min(min(min(temp)));
Ctmax=max(max(max(temp)));

%* define contour values for filled contours
cvT=[-2:2:30];					% temperature
cvS=[30:0.25:42];                % salinity

%* map plot windows
wind=[min(maplon) max(maplon) min(maplat) max(maplat)]; 
windslice=[min(maplat) max(maplat) min(Z) max(Z)];



%*****
% reference map to show location of meridional slice
%*****
figure(FNref)
%clf
axis(wind), axis image, hold on
plot(long, lat, 'k-')
plot(SL*ones(180,1), maplat, 'b-')



%*********
% meridional slices through the world ocean
%*********

%* salinity
salslice=sal(:,slicelon,:);			%extract values
salslice=squeeze(salslice);			%get rid of singleton column
salslice=salslice';					%transpose for ease of plotting

%* temperature
tempslice=temp(:,slicelon,:);		%extract values
tempslice=squeeze(tempslice);		%get rid of singleton column
tempslice=tempslice';				%rotate for ease of plotting



figure(FN)
clf
caxis([min(cvS) max(cvS)])
axis(windslice), hold on
[csS, hS]=contourf(maplat,Z, salslice, cvS);
clabel(csS, hS, 'manual')             % use return values to lable contours
xlabel('latitude')
ylabel('depth (m)') 
title('mean annual salinity in a meridional slice (ci = 0.25 p.s.u.)')
text(100, 0, int2str(SL))
colorbar('horiz')



figure(FN+1)
clf
caxis([Ctmin Ctmax])
axis(windslice), hold on
[csT, hT]=contourf(maplat,Z, tempslice, cvT);
clabel(csT, hT, 'manual')			% use return values to lable contours
xlabel('latitude')
ylabel('depth (m)') 
title('mean annual temperature in a meridional slice (ci = 2 deg C)')
text(100, 0, int2str(SL))
colorbar('horiz')




%* salinity and temperature slices on the same plot
%
%figure(FN+2)
%clf
%subplot(2,1,1)
%caxis([min(cvS) Csmax])
%axis(windslice), hold on
%pcolor(maplat,Z, salslice)
%shading flat
%ylabel('depth (m)')
%title('mean annual salinity in a meridional slice (p.s.u.)')
%text(100, 0, int2str(maplon(slicelon)))
%colorbar('horiz')

%subplot(2,1,2)
%caxis([Ctmin Ctmax])
%axis([min(maplat) max(maplat) min(Z) max(Z)]), hold on
%pcolor(maplat,Z, tempslice)
%shading flat
%xlabel('latitude')
%ylabel('depth (m)') 
%title('mean annual temperature in a meridional slice (deg C)')
%text(100, 0, int2str(SL))
%colorbar('horiz')