%  Routine for defining the PDO.  
%
%  Note that you will have to download some .m files from:
%  http://www.aos.wisc.edu/~dvimont/matlab
%  
%  Also, you will need to have the NetCDF toolboxes installed for MATLAB,
%  which is currently being done on the 14th floor computers.  I've
%  deposited the SST and SLP data in ~dvimont/matlab/Data on the 14th
%  floor computers if you would like to do this yourself.
%

clear

%  Start by loading SST data.  Define region:
lims = [120 240 20 60];
lev = 1;
tim = get_time(1948, 2003, 1948);

%  Load the data - note that this uses getnc2, which is avialable via the
%  website above.  You need the netcdf toolboxes to use this, so you can
%  also load the data from ~dvimont/matlab/Data

cd /Volumes/Data/NCEP/monthly   %  You will have to replace this line
[sst, lat, lon] = getnc2('skt.mon.mean.nc', 'skt', lims, 1, tim);
land = getnc2('land.sfc.gauss.nc', 'land', lims, 1, 1);

%  Now, I'll save the data to ~dvimont/matlab/Data (don't do this)
%
%  cd ~dvimont/matlab/Data
%  save sst.pdo.mat sst lat lon lims lev tim land
%
%  If you can't run the getnc2 command above, then load the sst data from
%  my home account.

%  Now, get rid of land points and ice points
lm = find(land == 1);
icept = find(min(sst) < 0);
lm = union(lm, icept);
land(lm) = 1;
ocept = find(land == 0);

%  Remove the annual cycle
[sst, sstc] = remove_anncyc(sst);

%  Weight data by sqrt(cos(lat)):
sstw = cosweight(sst, lat);

%  Now, perform EOF analysis
[lam, lds, pcs, per] = eof_dan(sstw(:,ocept));

%  We'll define the first PC as the PDO.

%  Now, generate regression map onto that PC
[rtem, ctem] = regress_eof(sst(:,ocept), -pcs(:,1));
regmap = repmat(NaN, [length(lat) length(lon)]);
regmap(ocept) = rtem;

%  And finally, plot it up nicely.  Again, you'll need a number of
%  routines from the website above.

%  Set up figure window
figure_tall(1); clf;

%  Define globals - YAX, XAX, and FRAME.  These are used by my mapping
%  routines.  
global_latlon(lat, lon, lims);

%  Now, set up subplot in the current window.
global_axes(6, 2.5, .5, .5, 1.5);
subplot2(1,1); cla;

%  Set up a map axis for the mapping toolbox - this makes pretty plots.
map_axis('giso');
map_contour_pn(regmap, 0.1, 'zero');

%  OK!  finish it up:
fill_landmap('over');
setm(gca, 'MLineLocation', 10, 'PLineLocation', 15);
gridm on; tightmap;
xlabel('Contour:  0.1 ^oC', 'fontsize', 12);
title('\bf PDO regression map', 'fontsize', 14);

%  Print this to a file
cd ~/Class/aos575/Homework/HW4
print -depsc2 -painters PDO_regmap.eps