Global Index (short | long) | Local contents | Local Index (short | long)
Do complex eof:
This script calls | |
---|---|
clear biff = 7; ind = 1; cd /home/disk/hayes2/dvimont/csiro/matlab_data if ind == 1; load 20degiso_yr101-550_temp_noave.mat lag = -60; ttit = ['101 - 550']; elseif ind == 2; load 20degiso_yr551-1000_temp_noave.mat lag = 90; ttit = ['551 - 1000']; end ctlim = [110 300 -60 60]; ntim = length(tim); nlat = length(lat); nlon = length(lon); temp3 = NaN * ones(ntim, nlat*nlon); temp3(:,kp20) = temp20; temp3 = detrend(temp3); temp3 = cosweight(temp3, lat); kp = find(~isnan(mean(temp3))); temp3 = temp3(:,kp); if biff == 1; yr1 = 10; yr2 = 10; tit = ['HP ( < 10 yr ) ']; elseif biff == 2; yr1 = 20; yr2 = 20; tit = ['HP ( < 20 yr ) ']; elseif biff == 3; yr1 = 10; yr2 = 60; tit = ['LP ( > 10 yr ) ']; elseif biff == 4; yr1 = 20; yr2 = 60; tit = ['LP ( > 20 yr ) ']; elseif biff == 5; yr1 = 8; yr2 = 20; tit = ['HP ( < 8 yr ) ']; elseif biff == 6; yr1 = 10; yr2 = 50; tit = ['BP ( 10-50 yr ) ']; elseif biff == 7; tit = ['Unfiltered Data ']; end [b1, a1] = butter(6, 2/yr1); [b2, a2] = butter(6, 2/yr2); if biff < 3 | biff == 5; temp3 = temp3 - filtfilt(b1, a1, temp3); disp('HP Filtered'); elseif biff > 2 & biff < 5; temp3 = filtfilt(b1, a1, temp3); disp('LP Filtered'); elseif biff ~= 7 temp3 = filtfilt(b1, a1, temp3) - filtfilt(b2, a2, temp3); disp('BP Filtered'); end [lam, lds, pcs, per] = eof_dan(temp3); %[lam, lds, pcs, per] = complex_eof(temp3); j = sqrt(-1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Get pcs from saved data clear cd /home/disk/hayes2/dvimont/csiro/matlab_data/20Deg_Isotherm/ ind = 1; ttit = ['101 - 550']; if ind == 1; load 20degiso_yr101-550_temp_noave.mat lag = -60; ttit = ['101 - 550']; elseif ind == 2; load 20degiso_yr551-1000_temp_noave.mat lag = 90; ttit = ['551 - 1000']; end ntim = length(tim); nlat = length(lat); nlon = length(lon); ctlim = [110 300 -60 60]; default_global % Get variables for regressions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% cd /home/disk/hayes2/dvimont/csiro/data filin = ['temp_A_L1-10.nc']; nc = netcdf(filin, 'nowrite'); lind = 1:10; depth = nc{'depth'}(lind); nc = close(nc); % Get regressions of fluxes: cd /home/disk/hayes2/dvimont/csiro/data nc = netcdf('zonal_stress_A_1000_years.nc', 'nowrite'); lat1 = nc{'latitude'}(:); lon1 = nc{'longitude'}(:); [xk, yk] = keep_var(FRAME, lon1, lat1); taux = nc{'smfzon'}(tim,yk,xk); mv = nc{'smfzon'}.missing_value(:); nc = close(nc); taux(find(taux == mv)) = NaN * ones(size(find(taux == mv))); nc = netcdf('merid_stress_A_1000_years.nc', 'nowrite'); tauy = nc{'smfmer'}(tim,yk,xk); mv = nc{'smfmer'}.missing_value(:); nc = close(nc); tauy(find(tauy == mv)) = NaN * ones(size(find(tauy == mv))); lat1 = lat1(yk); lon1 = lon1(xk); nc = netcdf('heat_flux_A_1000_years.nc', 'nowrite'); lat2 = nc{'latitude'}(:); lon2 = nc{'longitude'}(:); [xk, yk] = keep_var(FRAME, lon2, lat2); hflx = nc{'stfht'}(tim,yk,xk); mv = nc{'stfht'}.missing_value(:); nc = close(nc); hflx(find(hflx == mv)) = NaN * ones(size(find(hflx == mv))); lat2 = lat2(yk); lon2 = lon2(xk); nc = netcdf('temp_A_L1-10.nc', 'nowrite'); lat3 = nc{'latitude'}(:); lon3 = nc{'longitude'}(:); [xk, yk] = keep_var(FRAME, lon3, lat3); clear temp for i = 1:2; temp(:,i,:,:) = nc{'temp'}(tim,(4*(i-1)+1),yk,xk); end mv = nc{'temp'}.missing_value(:); nc = close(nc); temp(find(temp == mv)) = NaN * ones(size(find(temp == mv))); lat3 = lat3(yk); lon3 = lon3(xk); depth = depth([1 5]); [ntim, nlev, nlat, nlon] = size(temp); f = (2*7.292e-5)*sin(lat1*pi/180); clear curlt; for i = 1:ntim; temx = squeeze(taux(i,:,:)) ./ (f * ones(1, length(lon1))); temy = squeeze(tauy(i,:,:)) ./ (f * ones(1, length(lon1))); curlt(i,:,:) = sph_curl(temx, temy, lon1, lat1); % temx = squeeze(taux(i,:,:)) / 10; % temy = squeeze(tauy(i,:,:)) / 10; % curlt(i,:,:) = sph_div(temx, temy, lon1, lat1); end yk = keep_var3([-2 2], lat1); curlt(:,yk,:) = NaN; curlt = (1 / (10 * 1000)) * curlt; % This is in units m/s (Ekman vertical velocity) hflx = reshape(hflx, ntim, size(hflx,2)*size(hflx,3)); curlt = reshape(curlt, ntim, size(curlt,2)*size(curlt,3)); taux = reshape(taux, ntim, size(taux,2)*size(taux,3)); temp = reshape(temp, ntim, nlev, size(temp,3)*size(temp,4)); [hflx, clim] = remove_mean(hflx); [taux, clim] = remove_mean(taux); kpx = find(~isnan(clim)); [curlt, clim] = remove_mean(curlt); kpy = find(~isnan(clim)); [temp, clim] = remove_mean(temp); [temp20, clim] = remove_mean(temp20); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Load the pcs %%%%%%%%%%%%%%%%%%%%%%%%%%% for biff = [7]; biff cd /home/disk/hayes2/dvimont/csiro/matlab_data/20Deg_Isotherm/ if biff == 1; load 20deg_CEOF_HP10_yr101-550.mat yr1 = 10; yr2 = 10; tit = ['HP ( < 10 yr ) ']; ptit = 'HP10'; elseif biff == 3; load 20deg_CEOF_LP10_yr101-550.mat yr1 = 10; yr2 = 60; tit = ['LP ( > 10 yr ) ']; ptit = 'LP10'; elseif biff == 5; load 20deg_CEOF_HP8_yr101-550.mat yr1 = 5; yr2 = 20; tit = ['HP ( < 8 yr ) ']; ptit = 'HP8'; elseif biff == 7; load 20deg_CEOF_RAW_yr101-550.mat tit = ['Unfiltered Data ']; ptit = 'RAW'; end % Plot the regression maps: timeseries = sqrt(2)*pcs(:,1)./std(pcs(:,1)); cint1 = [0.05 0.05 0.05 1 1 1]; lag = 0; lg = lag*pi/180; num = 1; nfrm = 6; for i = 4:4; figure(i); clf; figure_orient(1); if i <= 2; tem1 = squeeze(temp(:,i,:)); cint = cint1(i); clev = [-1:cint:1]; xlab = [num2str(cint) ' K (std)^-^1']; XAX = lon3; YAX = lat3; nlat = length(YAX); nlon = length(XAX); tit2 = [num2str(depth(i)/100) 'm Temperature']; elseif i == 3; XAX = lon; YAX = lat; nlat = length(YAX); nlon = length(XAX); tem1 = NaN * ones(ntim, nlat*nlon); tem1(:,kp20) = temp20; cint = cint1(i); clev = [-1:cint:1]; xlab = [num2str(cint) ' K (std)^-^1']; tit2 = ['20 deg T''']; elseif i == 4; tem1 = 1e7*curlt; % tem1 = 1e9*curlt; cint = cint1(i); clev = [-10:cint:10]; xlab = [num2str(cint) ' x 10^-^7 m s^-^1']; % xlab = [num2str(cint) ' x 10^-^9 N m^-^3']; XAX = lon1; YAX = lat1; nlat = length(YAX); nlon = length(XAX); tit2 = ['W_e ( Ekman Upwelling )']; % tit2 = ['Div ( Wind Stress )']; elseif i == 5; tem1 = 1e2*taux; cint = cint1(i); clev = [-20:cint:20]; XAX = lon1; YAX = lat1; nlat = length(YAX); nlon = length(XAX); tit2 = ['Zonal Wind Stress']; xlab = [num2str(cint) ' x 10^-^3 N m^-^2']; elseif i == 6; tem1 = -1*hflx; cint = cint1(i); clev = [-20:cint:20]; xlab = [num2str(cint) ' W m^-^2 (std)^-^1']; XAX = lon2; YAX = lat2; nlat = length(YAX); nlon = length(XAX); tit2 = ['TOT Heat Flux (Pos. Up)']; end; % Filter Data: [b1, a1] = butter(6, 2/yr1); [b2, a2] = butter(6, 2/yr2); if biff < 3 | biff == 5; tem1 = tem1 - filtfilt(b1, a1, tem1); disp('HP Filtered'); elseif biff > 2 & biff < 5; tem1 = filtfilt(b1, a1, tem1); disp('LP Filtered'); elseif biff ~= 7 tem1 = filtfilt(b1, a1, tem1) - filtfilt(b2, a2, tem1); disp('BP Filtered'); end for ind = 1:nfrm; wgt = conj(exp(j * ((ind-1) * pi/nfrm + lg) )); temtim = real(wgt .* timeseries); tem = temtim'*tem1./ntim; % [ccoef, tstat] = corr_sig(temtim, tem1); ccoef = corr_sig(temtim, tem1); tem = reshape(tem, nlat, nlon); % tstat = reshape(tstat, nlat, nlon); tstat = reshape((abs(ccoef) >= 0.5), nlat, nlon); subplot(3,2,ind); cla; gshade2(abs(tstat), 0.9); hold on; gcont(tem, clev); hold off; dc; title(['Phase = ' num2str((ind-1)*180/nfrm + lag)]); if ind >= 5; xlabel(['Cont. Int.: ' xlab]); end end subplot(3,2,3); ylabel([tit ': ' tit2 ' regressed on CPC1 of 20deg T'', years 101-550; ( r >= 0.5 shaded )']); if 0; cd /home/disk/tao/dvimont/matlab/CSIRO/Plots5/yr101-550 if i == 1; eval(['print -dps2 ' ptit '_T1_CPC1_20deg_iso.ps']); elseif i == 2; eval(['print -dps2 ' ptit '_T5_CPC1_20deg_iso.ps']); elseif i == 3; eval(['print -dps2 ' ptit '_T20prime_CPC1_20deg_iso.ps']); elseif i == 4; eval(['print -dps2 ' ptit '_Wekman_CPC1_20deg_iso.ps']); elseif i == 5; eval(['print -dps2 ' ptit '_TauX_CPC1_20deg_iso.ps']); elseif i == 6; eval(['print -dps2 ' ptit '_HFlx_CPC1_20deg_iso.ps']); end end end end