Documentation of gr_nmc

Global Index (short | long) | Local contents | Local Index (short | long)

Help text

  Load and configure NMC.REANAL data

Cross-Reference Information

Listing of script gr_nmc

cd /home/disk/hayes2/dvimont/ccm/ccm3.6/data

%filin = ['sst_nov1981-jul1998.nc'];
filin = ['sst.mon.mean.nc'];
nc = netcdf(filin, 'nowrite');
  ts = nc{'air'}(:,:,:);
  datets = nc{'time'}(:);
  latts = nc{'lat'}(:);
  lonts = nc{'lon'}(:);
nc = close(nc)
[ntim, nlat, nlon] = size(ts);

%  Load landmask

filin = ['landmask.nc'];
nc = netcdf(filin, 'nowrite');
  lm = nc{'lsmask'}(:);
  lmlat = nc{'lat'}(:);
  lmlon = nc{'lon'}(:);
nc = close(nc);

%  Get CCM coordinates

filin = ['1983-12.nc'];
nc = netcdf(filin, 'nowrite');
  clat = nc{'lat'}(:);
  clon = nc{'lon'}(:);
nc = close(nc);

%  Interpolate everything to CCM coordinates

for i = 1:ntim
  cts(i,:,:) = interp2(lonts, latts, squeeze(ts(i,:,:)), clon, clat');
end
ts = cts;
clear cts

clm = interp2(lmlon, lmlat, lm, clon, clat');
clear lm

[ntim, nlat, nlon] = size(ts);

%  Reshape and get rid of land points

ts = reshape(ts, ntim, nlat*nlon);
clm = reshape(clm, 1, nlat*nlon);
ocpts = find(clm == 1);
ts = ts(:,ocpts);

%  Get SSTA's instead of SST's

[ts, clim] = annave(ts);

%  Calculate low pass filtered SSTA's

tslp = myrunning_ave(ts, 25);
tslp = myrunning_ave(tslp, 37);
tshp = ts - tslp;

%  Calculate ctstar and such

ctlim = [180 270 -6 6];
lonct = find(clon >= ctlim(1) & clon <= ctlim(2));
latct = find(clat >= ctlim(3) & clat <= ctlim(4));

tsnew = NaN*ones(ntim, nlat*nlon);
tsnew(:, ocpts) = tshp;
tsnew = reshape(tsnew, ntim, nlat, nlon);
ctstar = squeeze(mean(mean(shiftdim(tsnew(:,latct,lonct),1))));

tsnew = NaN*ones(ntim, nlat*nlon);
tsnew(:, ocpts) = tslp;
tsnew = reshape(tsnew, ntim, nlat, nlon);
ctlp   = squeeze(mean(mean(shiftdim(tsnew(:,latct,lonct),1))));

%clear tslp tshp
clear tsnew

%  Calculate GR time series
%    1.  Recalculate mean for ts

[ts, clim] = annave(ts, clim, 1);

%    2.  Calculate regression coefficients

[ts, clim] = annave(ts(37:(ntim-36),:));
%ntim = length(37:(ntim-36));
ctemp = ctstar(37:(ntim-36))-mean(ctstar(37:(ntim-36)));
a1 = (1/(ctemp'*ctemp))*ctemp'*ts;

%    3.  Calculate EOF of residual time series

[ntim, nkeep] = size(ts);

resid = ts - ctemp*a1;
tsnew = zeros(ntim, nlat*nlon);
tsnew(:,ocpts) = resid;
tsnew = reshape(tsnew, ntim, nlat, nlon);
tsnew = cosweight(tsnew, clat);
tsnew = reshape(tsnew, ntim, nlat*nlon);
resid = tsnew(:,ocpts);
%clear tsnew

%  Try this on the correlation matrix:
%resid = ones(ntim, 1)*(1./std(resid)) .* resid;

%  Try this by eliminating the poles:

tsnew = zeros(ntim, nlat*nlon);
tsnew(:,ocpts) = resid;
tsnew = reshape(tsnew, ntim, nlat, nlon);
latmask = find(clat <= -60 | clat >= 60);
tsnew(:, latmask, :) = zeros(ntim, length(latmask), nlon);
tsnew = reshape(tsnew, ntim, nlat*nlon);
resid = tsnew(:,ocpts);

c = (resid*resid')/nkeep;
[lam, pc, per] = eof(c);

pc = ones(ntim,1) * (1./(std(pc))) .* pc;

lds = pc(:,1:10)'*tsnew/ntim;
num = 1;
tem = zeros(1, nlat*nlon);
tem = lds(num,:);
tem = reshape(tem, nlat, nlon);

lims = [0 360 -90 90];

subplot_dan(1)
gcont(-1.*tem, lims,[-4:.1:4]);
dc
title('GR -- Covariance Matrix (Latitudes between 60S and 60N)')
xlabel('Contour Interval -- 0.1 K/std')
subplot_dan(2)
plot(-1.*pc(:,num))
title('GR')
xlabel('Year')
ylabel('STD')
set(gca, 'XTick',49:60:ntim,'XTickLabel',[1965:5:2000])
axis([1 408 -2.5 2.5])
grid

gr = pc(:,1);

ctpat = (1/std(ctemp)).*ctemp'*ts./ntim;
tem = zeros(1, nlat*nlon);
tem(ocpts) = ctpat;
tem = reshape(tem, nlat, nlon);

subplot_dan(1)
gcont(tem, lims, [-1:.2:2])
dc
title('CTSTAR -- Regression Pattern')
xlabel('Contour Interval -- 0.2 K/std')
subplot_dan(2)
plot(1/std(ctemp)*ctemp)
title('CTSTAR')
xlabel('Year')
ylabel('STD')
set(gca, 'XTick',49:60:ntim,'XTickLabel',[1965:5:2000])
axis([1 408 -3 3.5])
grid

%  NOTE:  std(ctemp) = 0.6815