Documentation of curlt_eof

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

Help text

  Do complex eof:

Cross-Reference Information

Listing of script curlt_eof

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);

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
ntim = length(tim); nlat = length(lat); nlon = length(lon);
default_global



%%%%%  Do the same for CURLT: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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);

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,:,:));
  temy = squeeze(tauy(i,:,:));
  curlt(i,:,:) = sph_div(temx, temy, lon1, lat1);
end
%yk2 = keep_var3([-10 10], lat1);
%curlt(:,yk2,:) = NaN;
curlt = (10 / 1000) * curlt; % This is in units m/s (Ekman vertical velocity)
[curlt, clim] = remove_mean(curlt);

curlt1 = cosweight(curlt, lat1);
curlt1 = reshape(curlt1, ntim, size(curlt1,2)*size(curlt1,3));
kp = find(~isnan(mean(curlt1)));
curlt1 = curlt1(:,kp);

[lam2, lds2, pcs2, per2] = eof_dan(curlt1);

curlt = reshape(curlt, ntim, size(curlt,2)*size(curlt,3));

%%%%%%%%% Plot the data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


sp(1);
XAX = lon; YAX = lat;
nlat = length(YAX); nlon = length(XAX);
temtim = pcs(:,1)./std(pcs(:,1));
tem = NaN*ones(1, nlat*nlon);
tem(kp20) = temtim' * temp20 / ntim;
tem = reshape(tem, nlat, nlon);
gcont(tem, [-0.5:.1:0.5]); dc;

sp(2);
XAX = lon1; YAX = lat1;
nlat = length(YAX); nlon = length(XAX);
temtim = pcs2(:,1)./std(pcs2(:,1));
tem = 1e5*temtim' * curlt / ntim;
tem = reshape(tem, nlat, nlon);
gcont(tem, [-10:10]); dc;