Global Index (short | long) | Local contents | Local Index (short | long)
temp = nc{'temp'}(tim, yk, xk);
This script calls | |
---|---|
clear cd /home/disk/hayes2/dvimont/csiro/data tim = [101:1000]; filin = ['temp_A_L1-10.nc']; nc = netcdf(filin, 'nowrite'); depth = nc{'depth'}(1:7); lat = nc{'latitude'}(:); lon = nc{'longitude'}(:); ctlim = [-.01 360 -90 90]; [xk, yk] = keep_var(ctlim, lon, lat); temp = nc{'temp'}(tim,1,yk,xk); mv = nc{'temp'}.missing_value(:); nc = close(nc); [ntim, nlev, nlat, nlon] = size(temp); temp(find(temp == mv)) = NaN * ones(size(find(temp == mv))); lat = lat(yk); lon = lon(xk); get_global; default_global; FRAME = ctlim; ctlim = [180 270 -6 6]; [xk, yk] = keep_var(ctlim, lon, lat); ct = squeeze(mean2(mean2(shiftdim(temp(:,:,yk,xk), 1)))); ct = detrend(ct); temp = reshape(temp, ntim, nlat*nlon); temp = detrend(temp); clim = mean(temp); kp = find(~isnan(clim)); temp = temp(:,kp); tim = [101:550]; ctem = detrend(ct(tim)); ctem = ctem./std(ctem); ctpat = ctem' * detrend(temp(tim, :)) ./ ntim; ct2 = ctpat * detrend(temp(tim,:))' ./ (length(kp)); ct2 = ct2./std(ct2); ctpat2 = ct2 * detrend(temp(tim, :)) ./ ntim; get_global; default_global; figure(1) sp(1) tem = NaN * ones(1, nlat*nlon); tem(kp) = ctpat; tem = reshape(tem, nlat, nlon); gcont(tem, [-.5:.05:.5]); dc sp(2) tem = NaN * ones(1, nlat*nlon); tem(kp) = ctpat2; tem = reshape(tem, nlat, nlon); gcont(tem, [-.5:.05:.5]); dc for i = 1:3; if i == 1; ctann = pcs1(:,1)/std(pcs1(:,1)); tit = 'PC1, Years 101 - 550'; elseif i == 2; ctann = pcs2(:,1)/std(pcs2(:,1)); tit = 'PC1, Years 551 - 1000'; elseif i == 3; ctann = pcs3(:,1)/std(pcs3(:,1)); tit = 'PC1, Years 101 - 1000'; end nx = length(ctann); nfft = 128; noverlap = nfft/2; [p, f] = spectrum(ctann, nfft, noverlap); n2 = nfft/2; f2 = 0.5 * (0:n2)/n2; rho = (corr(ctann, ctann, 1) + sqrt(corr(ctann, ctann, 2))) / 2 % + ... % (corr(ctann, ctann, 3)^(1/3))) / 3; rn = (1-rho^2) ./ (1-2*rho*cos((pi/n2)*(0:n2))+rho^2); rn = rn / mean(rn); rn = rn * mean(p(:,1)); if i ~= 3; rner1 = rn * 2.64; rner5 = rn * 2.01; else rner1 = rn * 2.05; rner5 = rn * 1.68; end dofx = nx / n2; figure(2); figure_orient; subplot(3,1,i) semilogy(f2, p(:,1), 'b-', f2, rn, 'b-', ... f2, rner5, 'b--', f2, rner1, 'b-.') % tim = 2:(n2+1); f2 = f2(tim); p = p(tim,:); rn = rn(tim); % rner1 = rner1(tim); rner5 = rner5(tim); % plot(1*log(f2), f2'.*p(:,1), 'b-', 1*log(f2), f2.*rn, 'b-', ... % 1*log(f2), f2.*rner5, 'b--', 1*log(f2), f2.*rner1, 'b-.') axis([0 0.5 .05 10]) grid % title(['Power Spectrum for PC' num2str(num) ', 30s to 30n, ' tit]); % title(['Power Spectrum for CT' num2str(fig) ', ' tit]); title(tit) xlabel(['Frequency: yr^-^1; Degrees of Freedom: ' ... num2str(round(dofx)) ';' ... ' Bandwidth: 7.8 x 10^-^3 yr^-^1']) % l=legend('CT Spectrum', 'Red Noise', '95% Confidence', '99% Confidence'); end cd /home/disk/tao/dvimont/matlab/CSIRO/Plots3 % Check out EOF of entire globe's SST temp = squeeze(shiftdim(temp, 1)); for i = 1:nlat/2 indy = 2*(i-1)+[1:2]; for j = 1:nlon/2 indx = 2*(j-1)+[1:2]; sst(:,i,j) = squeeze(mean(mean(temp(indy,indx,:)))); end end lat = mean([lat(1:2:nlat)'; lat(2:2:nlat)']); lon = mean([lon(1:2:nlon)'; lon(2:2:nlon)']); tim = [451:900]; tim = [1:450]; tim = [1:900]; temp = sst(tim,:,:); [ntim, nlat, nlon] = size(temp); temp = cosweight(temp, lat); temp = reshape(temp, ntim, nlat*nlon); clim = mean(temp); kp = find(~isnan(clim)); temp = temp(:,kp); temp = detrend(temp); [lam, per, lds, pcs] = eof_dan(temp); %pcs1 = pcs; %pcs2 = pcs; pcs3 = pcs; ctem = pcs1(:,1); ct2 = pcs2(:,1);