%start to compose the atmospheric energy budget year by year

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; close all;


dir='E:\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; %this is going to be used to define the start index for deacde averages in each file
NF=length(fil);
for n=1:2;%NF;
    tic
FIL=[dir char(fil(n,:))]; %read in file
ncdisp(FIL);

t=ncread(FIL, 'time', [1], [Inf]); NT=length(t); %syntax is ( file name as a character, [start index], [ number to grab], [number to skip]);
lon=ncread(FIL, 'lon', [1], [Inf]); N=length(lon);
lat=ncread(FIL,'lat', [1], [Inf]); NN=length(lat);
precl=permute(ncread(FIL, 'PRECL', [1 1 1], [Inf Inf Inf]), [3 2 1]); %permute changes the order, [ start lon index, start lat index, start time index]
year=floor(t*1000);%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

t_step=10; %^to take decade averages
yrs=min(year):t_step:max(year); NY=length(yrs);
% cut=20; %this code was for finding ITCZ location - not relevant here
% don=find(lat>=0 & lat <= cut);
% dos=find(lat<=0 & lat >= -cut);

for ny=1:NY; %indexes in the deacde avaerages -n referenced to the years in yrs
    do = find( year >= yrs(ny) & year < yrs(ny)+t_step ); %year vector is referenced to monthly data
    length_do(ny)= length(do);
    PRECL(off+ny,1:NN, 1:N)= (squeeze(mean(precl(do,:,:))));
    
end
off=off+NY;
[min(year) max(year) off]
clear precl
toc
end

YEARS=-22000:t_step:39;
%save('F:\trace\PRECL_dec.mat', 'YEARS', 'lon','lat', 'PRECL');
break
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%start to compose the atmospheric energy budget year by year

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; close all;

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);
precc=permute(ncread(FIL, 'PRECC', [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

t_step=10;
yrs=min(year):t_step:max(year); NY=length(yrs);
cut=20;
don=find(lat>=0 & lat <= cut);
dos=find(lat<=0 & lat >= -cut);

for ny=1:NY;
    do = find( year >= yrs(ny) & year < yrs(ny)+t_step );
    PRECC(off+ny,1:NN, 1:N)= (squeeze(mean(precc(do,:,:))));
    
end
off=off+NY;
[min(year) max(year) off]
clear precc
toc
end

YEARS=-22000:t_step:39;
save('F:\trace\PRECC_dec.mat', 'YEARS', 'lon','lat', 'PRECC');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;close all;
load('F:\trace\PRECC_dec.mat', 'YEARS', 'lon','lat', 'PRECC');

load('F:\trace\PRECL_dec.mat', 'YEARS', 'lon','lat', 'PRECL');

SF=1000*3600*24; %converts from m/s to mm/day
PRECC(find(PRECC*SF>100))=NaN; %some messed up values snuck in there
P=PRECC+PRECL;

[NT, NN,N]=size(P);
ZP=squeeze(mean(shiftdim(P,2)));
cut_lat=20;
flat=-cut_lat:.05:cut_lat; %fine mesh of latitudes in the tropics
for nt=1:NT;
    FL=squeeze(ZP(nt,:));
    fP=interp1(lat,FL, flat);
    cdf=cumsum(fP)/sum(fP);
    pcent(nt)=flat(find2(cdf,0.5));
end

gray=[0.5 0.5 0.5];
LWD=1.5; LWL=2.5; LWC=2.4;
FSz=24; FN='Arial';
dgreen=[34 139 34]/255;
low_pcent=low_filter(pcent, 1,20); %200 year period low pass filter
figure(1);clf;
hold on; plot(YEARS, pcent, 'Color',gray, 'LineWidth',LWD); 
hold on; plot(YEARS, low_pcent, 'Color','k', 'LineWidth',LWL); 
xticl=-20:5:0;
xtic=1000*xticl;
xlabel('Years X1000', 'FontSize', FSz, 'FontName', FN)
ylabel('Precipitation Centroid (^o)', 'FontSize', FSz, 'FontName', FN);
set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2.5, 'XTick', xtic, 'XTickLabel', xticl);
box on;

ts=low_pcent;
ts=(ts-mean(ts))/std(ts);
num_std=2;
for nn=1:NN;
    for n=1:N;
        P_map(nn,n) = mean(squeeze(P(:,nn,n)).*ts');
    end
end

sphcolor2(num_std*SF*double(P_map),lat,lon,3);
[LON LAT]=meshgrid(lon,lat);
FL=SF*double(squeeze(mean(P)));
hold on; contourm(LAT,LON, FL, [ 4  8 12 16], 'Color', dgreen,'LineWidth', LWC) 
load coast;
H=colorbar('LineWidth', 2.5);
hold on; plotm(lat,long, 'k-', 'LineWidth', 2.5)
load('F:\trace\PRECC_dec.mat', 'YEARS', 'lon','lat', 'PRECC');

%and zonal means
CP=SF*squeeze(mean(ZP));
PP=SF*(squeeze(mean(ZP))+num_std*mean(P_map'));
NP=SF*(squeeze(mean(ZP))-num_std*mean(P_map'));

LW=2.4;
xmx=6;
figure(6);clf;
subplot(1,2,1);
hold on; plot(CP,lat,'k-', 'LineWidth', LW);
hold on; plot(PP,lat,'r-', 'LineWidth', LW);
hold on; plot(NP,lat,'b-', 'LineWidth', LW);
axis([0 xmx -90 90]);
set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2.4, 'YTick', [-80:20:80]); box on;
hold on; plot([0 xmx], [0 0], 'm-', 'LineWidth', 2.4)

%define the local precipitation centroid and the precip centroid for
%positive and negative anomalies
FL=squeeze(mean(P));
PFL=squeeze(mean(P)) + num_std*P_map;
NFL=squeeze(mean(P)) - num_std*P_map;

for n=1:N;
    zp=squeeze(FL(:,n));
    pzp=squeeze(PFL(:,n));
    nzp=squeeze(NFL(:,n));
    
    fP=interp1(lat,zp,flat);
    fPP=interp1(lat, pzp, flat);
    fNP=interp1(lat, nzp, flat);
    
    cdf=cumsum(fP)/sum(fP);
    pcdf=cumsum(fPP)/sum(fPP);
    ncdf=cumsum(fNP)/sum(fNP);
    
    loc_PCENT(n) = flat(find2(cdf,0.5));
    P_loc_PCENT(n) = flat(find2(pcdf,0.5));
    N_loc_PCENT(n) = flat(find2(ncdf,0.5));
end


figure(3); hold on;
hold on; plotm(loc_PCENT',lon, 'k--', 'LineWidth', LW);
hold on; plotm(P_loc_PCENT',lon, 'r--', 'LineWidth', LW);
hold on; plotm(N_loc_PCENT', lon,'b--', 'LineWidth', LW);
LWA=1.8;
hold on; plotm(0*ones(size(lon)),lon,'m-','LineWidth', LWA);
hold on; plotm(20*ones(size(lon)),lon,'m--','LineWidth', LWA);
hold on; plotm(-20*ones(size(lon)),lon,'m--','LineWidth', LWA);

%lastly, the precipitation centroids
zp=mean(squeeze(mean(P))');
pzp=mean( (squeeze(mean(P)) + num_std*P_map)');
nzp=mean( (squeeze(mean(P)) - num_std*P_map)');

fP=interp1(lat,zp,flat);
fPP=interp1(lat, pzp, flat);
fNP=interp1(lat, nzp, flat);

cdf=cumsum(fP)/sum(fP);
pcdf=cumsum(fPP)/sum(fPP);
ncdf=cumsum(fNP)/sum(fNP);
    
 P_CENT = flat(find2(cdf,0.5));
 P_P_CENT = flat(find2(pcdf,0.5));
 N_P_CENT = flat(find2(ncdf,0.5));

 %figure(6);
 %hold on; plot([0 xmx], P_CENT*[1 1],'k--', 'LineWidth', LW);
 %hold on; plot([0 xmx], P_P_CENT*[1 1],'r--', 'LineWidth', LW);
 %hold on; plot([0 xmx], N_P_CENT*[1 1],'b--', 'LineWidth', LW);