clear all; close all;
offset=-0.000001;
%make the precip from the trace experiment
dir='F:\trace\PRECL\';
fil={'trace.01.22000-20001BP.cam2.h0.PRECL.0000101-0200012.nc';  
'trace.02.20000-19001BP.cam2.h0.PRECL.0200101-0300012.nc';  
'trace.03.19000-18501BP.cam2.h0.PRECL.0300101-0350012.nc';  
'trace.04.18500-18401BP.cam2.h0.PRECL.0350101-0360012.nc';  
'trace.05.18400-17501BP.cam2.h0.PRECL.0360101-0450012.nc';  
'trace.06.17500-17001BP.cam2.h0.PRECL.0450101-0500012.nc';  
'trace.07.17000-16001BP.cam2.h0.PRECL.0500101-0600012.nc';  
'trace.08.16000-15001BP.cam2.h0.PRECL.0600101-0700012.nc';  
'trace.09.15000-14901BP.cam2.h0.PRECL.0700101-0710012.nc';  
'trace.10.14900-14351BP.cam2.h0.PRECL.0710101-0765012.nc';  
'trace.11.14350-13871BP.cam2.h0.PRECL.0765101-0813012.nc';  
'trace.12.13870-13101BP.cam2.h0.PRECL.0813101-0890012.nc';  
'trace.13.13100-12901BP.cam2.h0.PRECL.0890101-0910012.nc';  
'trace.14.12900-12501BP.cam2.h0.PRECL.0910101-0950012.nc';  
'trace.15.12500-12001BP.cam2.h0.PRECL.0950101-1000012.nc'; 
'trace.16.12000-11701BP.cam2.h0.PRECL.1000101-1030012.nc'; 
'trace.17.11700-11301BP.cam2.h0.PRECL.1030101-1070012.nc';  
'trace.18.11300-10801BP.cam2.h0.PRECL.1070101-1120012.nc';  
'trace.19.10800-10201BP.cam2.h0.PRECL.1120101-1180012.nc';  
'trace.20.10200-09701BP.cam2.h0.PRECL.1180101-1230012.nc';  
'trace.21.09700-09201BP.cam2.h0.PRECL.1230101-1280012.nc';  
'trace.22.09200-08701BP.cam2.h0.PRECL.1280101-1330012.nc';  
'trace.23.08700-08501BP.cam2.h0.PRECL.1330101-1350012.nc';  
'trace.24.08500-08001BP.cam2.h0.PRECL.1350101-1400012.nc';  
'trace.25.08000-07601BP.cam2.h0.PRECL.1400101-1440012.nc';  
'trace.26.07600-07201BP.cam2.h0.PRECL.1440101-1480012.nc';  
'trace.27.07200-06701BP.cam2.h0.PRECL.1480101-1530012.nc';  
'trace.28.06700-06201BP.cam2.h0.PRECL.1530101-1580012.nc';  
'trace.29.06200-05701BP.cam2.h0.PRECL.1580101-1630012.nc';  
'trace.30.05700-05001BP.cam2.h0.PRECL.1630101-1700012.nc';  
'trace.31.05000-04001BP.cam2.h0.PRECL.1700101-1800012.nc';  
'trace.32.04000-03201BP.cam2.h0.PRECL.1800101-1880012.nc';  
'trace.33.03200-02401BP.cam2.h0.PRECL.1880101-1960012.nc';  
'trace.34.02400-01401BP.cam2.h0.PRECL.1960101-2060012.nc';  
'trace.35.01400-00401BP.cam2.h0.PRECL.2060101-2160012.nc';  
'trace.36.400BP-1990CE.cam2.h0.PRECL.2160101-2204012.nc';};
off=0;
NF=length(fil);
for n=1:NF;
    tic
FIL=[dir char(fil(n,:))];
%ncdisp(FIL);
t=ncread(FIL, 'time', [1], [Inf]); NT=length(t);
lon=ncread(FIL, 'lon', [1], [Inf]); N=length(lon);
lat=ncread(FIL,'lat', [1], [Inf]); NN=length(lat);
PL=permute(ncread(FIL, 'PRECL', [1 1 1], [Inf Inf Inf]), [3 2 1]);
year=rounddown((t+offset)*1000);
month=round(12*((t+offset)*1000-(year)));
if n == 36; %the times are just a little fubarred
    year=g2c((-400:39)'*ones(1,12));
    month=g2c( ((1:12)'*ones(1,length(-400:39)))');
end
yrs=min(year):max(year); NY=length(yrs);
SF=1000*24*3600;
for ny=1:NY;
    for mn=1:12;
    do = find( year == yrs(ny) & month == mn);
    ZPL(off+ny,mn, 1:NN)= mean(squeeze(squeeze(PL(do,:,:)))');
    end
end
off=off+NY;
clear PL
[min(year) max(year) off]
toc
end

YEARS=-2200:40;
save('F:\trace\zonal_PL_seas.mat', 'YEARS', 'lat', 'ZPL');

break
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; close all;
%make the precip from the trace experiment
offset=-0.000001;
dir='F:\trace\PRECC\';
fil={'trace.01.22000-20001BP.cam2.h0.PRECC.0000101-0200012.nc'  
'trace.02.20000-19001BP.cam2.h0.PRECC.0200101-0300012.nc';  
'trace.03.19000-18501BP.cam2.h0.PRECC.0300101-0350012.nc';  
'trace.04.18500-18401BP.cam2.h0.PRECC.0350101-0360012.nc'; 
'trace.05.18400-17501BP.cam2.h0.PRECC.0360101-0450012.nc';  
'trace.06.17500-17001BP.cam2.h0.PRECC.0450101-0500012.nc';  
'trace.07.17000-16001BP.cam2.h0.PRECC.0500101-0600012.nc';  
'trace.08.16000-15001BP.cam2.h0.PRECC.0600101-0700012.nc';  
'trace.09.15000-14901BP.cam2.h0.PRECC.0700101-0710012.nc';  
'trace.10.14900-14351BP.cam2.h0.PRECC.0710101-0765012.nc';  
'trace.11.14350-13871BP.cam2.h0.PRECC.0765101-0813012.nc';  
'trace.12.13870-13101BP.cam2.h0.PRECC.0813101-0890012.nc';  
'trace.13.13100-12901BP.cam2.h0.PRECC.0890101-0910012.nc';  
'trace.14.12900-12501BP.cam2.h0.PRECC.0910101-0950012.nc';  
'trace.15.12500-12001BP.cam2.h0.PRECC.0950101-1000012.nc';  
'trace.16.12000-11701BP.cam2.h0.PRECC.1000101-1030012.nc';  
'trace.17.11700-11301BP.cam2.h0.PRECC.1030101-1070012.nc';  
'trace.18.11300-10801BP.cam2.h0.PRECC.1070101-1120012.nc';  
'trace.19.10800-10201BP.cam2.h0.PRECC.1120101-1180012.nc';  
'trace.20.10200-09701BP.cam2.h0.PRECC.1180101-1230012.nc';  
'trace.21.09700-09201BP.cam2.h0.PRECC.1230101-1280012.nc';  
'trace.22.09200-08701BP.cam2.h0.PRECC.1280101-1330012.nc';  
'trace.23.08700-08501BP.cam2.h0.PRECC.1330101-1350012.nc';  
'trace.24.08500-08001BP.cam2.h0.PRECC.1350101-1400012.nc';  
'trace.25.08000-07601BP.cam2.h0.PRECC.1400101-1440012.nc';  
'trace.26.07600-07201BP.cam2.h0.PRECC.1440101-1480012.nc';  
'trace.27.07200-06701BP.cam2.h0.PRECC.1480101-1530012.nc';  
'trace.28.06700-06201BP.cam2.h0.PRECC.1530101-1580012.nc';  
'trace.29.06200-05701BP.cam2.h0.PRECC.1580101-1630012.nc';  
'trace.30.05700-05001BP.cam2.h0.PRECC.1630101-1700012.nc';  
'trace.31.05000-04001BP.cam2.h0.PRECC.1700101-1800012.nc';  
'trace.32.04000-03201BP.cam2.h0.PRECC.1800101-1880012.nc';  
'trace.33.03200-02401BP.cam2.h0.PRECC.1880101-1960012.nc';  
'trace.34.02400-01401BP.cam2.h0.PRECC.1960101-2060012.nc';  
'trace.35.01400-00401BP.cam2.h0.PRECC.2060101-2160012.nc';  
'trace.36.400BP-1990CE.cam2.h0.PRECC.2160101-2204012.nc';}   
off=0;
NF=length(fil);
for n=1:NF;
    tic
FIL=[dir char(fil(n,:))];
%ncdisp(FIL);
t=ncread(FIL, 'time', [1], [Inf]); NT=length(t);
lon=ncread(FIL, 'lon', [1], [Inf]); N=length(lon);
lat=ncread(FIL,'lat', [1], [Inf]); NN=length(lat);
PC=permute(ncread(FIL, 'PRECC', [1 1 1], [Inf Inf Inf]), [3 2 1]);
year=rounddown((t+offset)*1000);
month=round(12*((t+offset)*1000-(year)));
if n == 36; %the times are just a little fubarred
    year=g2c((-400:39)'*ones(1,12));
    month=g2c( ((1:12)'*ones(1,length(-400:39)))')
end
yrs=min(year):max(year); NY=length(yrs);
SF=1000*24*3600;
for ny=1:NY;
    for mn=1:12;
    do = find( year == yrs(ny) & month==mn );
    MNS(off+ny, 1:length(do)) = month(do);
    ZPC(off+ny,mn, 1:NN)= mean(squeeze(squeeze(PC(do,:,:)))');
    end
end
off=off+NY;
clear PL
[min(year) max(year) off]
toc
end

YEARS=-22000:40;
save('F:\trace\zonal_PC_seas.mat', 'YEARS', 'lat', 'ZPC');

break
% 
% 
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% clear all; close all;
% %make the precip from the trace experiment
% FIL='F:\trace\trace.36.400BP-1990CE.cam2.h0.LANDFRAC.2160101-2204012.nc';
% %ncdisp(FIL);
% LAND=double(squeeze(mean(permute(ncread(FIL, 'LANDFRAC', [1 1 1], [Inf Inf Inf]), [3 2 1])))); %time mean land fraction over the last 500 ish years
% lat=ncread(FIL, 'lat', [1], [Inf]); NN=length(lat);
% lon=ncread(FIL, 'lon', [1], [Inf]); N=length(lon);
% sphcolor2jet(LAND, lat,lon,3);
% OMASK=round(1-LAND);
% LMASK=round(LAND);
% OMASK2=OMASK; OMASK2(find(OMASK2<0.5))= NaN;
% LMASK2=LMASK; LMASK2(find(LMASK2<0.5))= NaN;
% 
% dir='F:\trace\TS\';
% fil={'trace.01.22000-20001BP.cam2.h0.TS.0000101-0200012.nc';     
% 'trace.02.20000-19001BP.cam2.h0.TS.0200101-0300012.nc';     
% 'trace.03.19000-18501BP.cam2.h0.TS.0300101-0350012.nc';     
% 'trace.04.18500-18401BP.cam2.h0.TS.0350101-0360012.nc';     
% 'trace.05.18400-17501BP.cam2.h0.TS.0360101-0450012.nc';     
% 'trace.06.17500-17001BP.cam2.h0.TS.0450101-0500012.nc';     
% 'trace.07.17000-16001BP.cam2.h0.PRECL.0500101-0600012.nc';  
% 'trace.07.17000-16001BP.cam2.h0.TS.0500101-0600012.nc';     
% 'trace.09.15000-14901BP.cam2.h0.TS.0700101-0710012.nc';     
% 'trace.10.14900-14351BP.cam2.h0.TS.0710101-0765012.nc';     
% 'trace.11.14350-13871BP.cam2.h0.TS.0765101-0813012.nc';     
% 'trace.12.13870-13101BP.cam2.h0.TS.0813101-0890012.nc';     
% 'trace.13.13100-12901BP.cam2.h0.TS.0890101-0910012.nc';     
% 'trace.14.12900-12501BP.cam2.h0.TS.0910101-0950012.nc';     
% 'trace.15.12500-12001BP.cam2.h0.TS.0950101-1000012.nc';     
% 'trace.16.12000-11701BP.cam2.h0.TS.1000101-1030012.nc';     
% 'trace.17.11700-11301BP.cam2.h0.TS.1030101-1070012.nc';     
% 'trace.18.11300-10801BP.cam2.h0.TS.1070101-1120012.nc';     
% 'trace.19.10800-10201BP.cam2.h0.TS.1120101-1180012.nc';     
% 'trace.20.10200-09701BP.cam2.h0.TS.1180101-1230012.nc';     
% 'trace.21.09700-09201BP.cam2.h0.TS.1230101-1280012.nc';     
% 'trace.22.09200-08701BP.cam2.h0.TS.1280101-1330012.nc';     
% 'trace.23.08700-08501BP.cam2.h0.TS.1330101-1350012.nc';     
% 'trace.24.08500-08001BP.cam2.h0.TS.1350101-1400012.nc';     
% 'trace.25.08000-07601BP.cam2.h0.TS.1400101-1440012.nc';     
% 'trace.26.07600-07201BP.cam2.h0.TS.1440101-1480012.nc';     
% 'trace.27.07200-06701BP.cam2.h0.TS.1480101-1530012.nc';     
% 'trace.28.06700-06201BP.cam2.h0.TS.1530101-1580012.nc';     
% 'trace.29.06200-05701BP.cam2.h0.TS.1580101-1630012.nc';     
% 'trace.30.05700-05001BP.cam2.h0.TS.1630101-1700012.nc';     
% 'trace.31.05000-04001BP.cam2.h0.TS.1700101-1800012.nc';     
% 'trace.32.04000-03201BP.cam2.h0.TS.1800101-1880012.nc';     
% 'trace.33.03200-02401BP.cam2.h0.TS.1880101-1960012.nc';     
% 'trace.34.02400-01401BP.cam2.h0.TS.1960101-2060012.nc';     
% 'trace.35.01400-00401BP.cam2.h0.TS.2060101-2160012.nc';     
% 'trace.36.400BP-1990CE.cam2.h0.TS.2160101-2204012.nc'};   
% off=0;
% NF=length(fil);
% for n=1:NF;
%     tic
% FIL=[dir char(fil(n,:))];
% %ncdisp(FIL);
% t=ncread(FIL, 'time', [1], [Inf]); NT=length(t);
% lon=ncread(FIL, 'lon', [1], [Inf]); N=length(lon);
% lat=ncread(FIL,'lat', [1], [Inf]); NN=length(lat);
% TS=permute(ncread(FIL, 'TS', [1 1 1], [Inf Inf Inf]), [3 2 1]);
% year=rounddown(t*1000);
% month=round(12*(t*1000-(year)));
% if n == 36; %the times are just a little fubarred
%     year=g2c((-400:39)'*ones(1,12));
%     month=g2c( ((1:12)'*ones(1,length(-4000:39)))')
% end
% 
% 
% yrs=min(year):max(year); NY=length(yrs);
% SF=1000*24*3600;
% cut=20;
% don=find(lat>=0 & lat <= cut);
% dos=find(lat<=0 & lat >= -cut);
% 
% for ny=1:NY;
%     do = find( year == yrs(ny) );
%     ZTS(off+ny, 1:NN)= mean(squeeze(mean(TS(do,:,:)))');
%     ZSST(off+ny, 1:NN)= nanmean((squeeze(mean(TS(do,:,:))).*OMASK2)');
%     ZLTS(off+ny, 1:NN)= nanmean((squeeze(mean(TS(do,:,:))).*LMASK2)');
%     
%     grad_SST(off+ny) = nanmean(g2c(squeeze(mean(TS(do,don,:))).*OMASK2(don,:))) - nanmean(g2c(squeeze(mean(TS(do,dos,:))).*OMASK2(dos,:))); %tropical SST gradient
%     
%     
% end
% off=off+NY;
% clear PL
% [min(year) max(year) off]
% toc
% end
% 
% YEARS=-22000:39;
% save('F:\trace\TS.mat', 'YEARS', 'lat', 'ZTS', 'ZSST', 'ZLTS', 'grad_SST');
% 
clear all;close all;
do=1:22000;
load('F:\trace\zonal_PC.mat', 'YEARS', 'lat', 'ZPC');
load('F:\trace\zonal_PL.mat', 'YEARS', 'lat', 'ZPL');
load('F:\trace\TS.mat', 'YEARS', 'lat', 'ZTS', 'ZSST', 'ZLTS', 'grad_SST');
grad_SST=double(grad_SST);
YEARS=-22000:39;
P=double(ZPC+ZPL);
clear ZPC ZPL
NY=length(YEARS);
flat=-20:0.05:20;
gray=[0.5 0.5 0.5];
lred=[0.6 0 0];
for ny=1:NY;
    FL=squeeze(P(ny,:));
    fP=interp1(lat,FL,flat);
    cdf=cumsum(fP)/sum(fP);
    PCENT(ny)=flat(find2(cdf,0.5));
end

LW=0.7; LWL=2;
FSz=22; FN='Arial';

figure(1); clf;
plot(YEARS(do), PCENT(do), '-','Color', gray, 'LineWidth', LW);
%low pass filter
cut_per=200;
low_PCENT = low_filter(PCENT, 1, cut_per);
hold on; plot(YEARS(do), low_PCENT(do), 'k', 'LineWidth', LWL)
SF_SST=2;
hold on; plot(YEARS(do), SF_SST*grad_SST(do), 'Color',lred,'LineWidth',LW);
low_SST = low_filter(grad_SST, 1, cut_per);
hold on; plot(YEARS(do), SF_SST*low_SST(do), 'r','LineWidth',LWL);

xlabel('Thousand Years','FontSize', FSz, 'FontName', FN)
ylabel('ITCZ location (P_C_E_N_T in ^o)', 'FontSize', FSz, 'FontName', FN);
xticl=-20:5:5;
xtic=xticl*1000;
set(gca, 'FontSize', FSz, 'FontName',FN, 'XTick', xtic, 'XTickLabel', xticl, 'LineWidth', 2.4)
xmn=-22000; xmx=2000;
ymn=-3.4; ymx=1;
title('CCSM trace timeseries of ITCZ position', 'FontSize', FSz, 'FontName', FN);
axis([xmn xmx ymn ymx]); box on;


