% compute the toa atmospheric energy budegt for the 24 mm slab ocean
% aquaplanet
clear all; close all;

%constants
a=6.36*10^6;
g=9.81;
cp=1004;
L=2.51*10^6; %latent heat of evaporation in J/kg

%model data
startyear=1987;
endyear=1991;
for yr=startyear:endyear;
dir = '/home/disk/eos10/dargan/aquafv/output/seas24mml/history/'; 
fil =[ num2str(yr), '0101.atmos_month.nc'];
FIL=[dir fil];

nc=netcdf(FIL);
%ncdump(FIL);

lon = nc{'lon'}(:); N=length(lon);
lat = nc{'lat'}(:); NN=length(lat);
Phalf = nc{'phalf'}(:); 
Pfull = nc{'pfull'}(:); NL=length(Pfull);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  TOA %%%%%%%%%%%%%%%%%%%%%%%%%%
t(yr-startyear+1,1:12,1:NL,1:NN,1:N) = nc{'temp'}(:,:,:,:); %4d temp
ts(yr-startyear+1,1:12,1:NN,1:N) = nc{'t_surf'}(:,:,:); %surface temp temp

end



TS=squeeze(mean(ts));
T=squeeze(mean(t));
save('T.mat', 'T','lat', 'lon', 'Pfull')
break

cut=max(find(Phalf < 200)); %tropopause index
%define the mass weigthing vector (normalized)
W=(Phalf(2:end) - Phalf(1:end-1)) / (max(Phalf - min(Phalf)));
W(1:(cut-1))=0; %leave out stratosphere
W=W/sum(W); % renormalize

W9 = zeros(size(W)); W9(18:19) = .5; %weighting function for low level temp
W8 = zeros(size(W)); W8(16:17) = .5; %for 800 hPa temp
for  nt=1:12, 
    for nn=1:NN;
        for n=1:N;
            TV(nt,nn,n) = squeeze(T(nt,:,nn,n)) * W; %vertically averaged temperature
            T9(nt,nn,n) = squeeze(T(nt,:,nn,n)) * W9; %900 hPa temperature
            T8(nt,nn,n) = squeeze(T(nt,:,nn,n)) * W8; %800 hPa temperature
        end
    end
end



%make the polar and equatorward averaged temperatures
TVZ =squeeze(mean(shiftdim(TV,2))); %zonal mean vertically integrated temp
TSZ =squeeze(mean(shiftdim(TS,2))); %ditto surface temp.
cross = 30; %crossover latitude
don=find(lat>=30); dot = find(lat > -cross & lat < cross); dos = find(lat <= -cross);

GTV=trapz(sin(lat(:)*pi/180), squeeze(mean(TVZ))')/2;
GTS=trapz(sin(lat(:)*pi/180), squeeze(mean(TSZ))')/2;
for mn=1:12;
    TVN (mn) = trapz(sin(lat(don)*pi/180), squeeze(TVZ(mn,don))')/trapz(sin(lat(don)*pi/180), ones(size(squeeze(TVZ(mn,don))'))) - GTV;
    TVT (mn) = trapz(sin(lat(dot)*pi/180), squeeze(TVZ(mn,dot))')/trapz(sin(lat(dot)*pi/180), ones(size(squeeze(TVZ(mn,dot))'))) - GTV;
    TVS (mn) = trapz(sin(lat(dos)*pi/180), squeeze(TVZ(mn,dos))')/trapz(sin(lat(dos)*pi/180), ones(size(squeeze(TVZ(mn,dos))'))) - GTV;
    
    TSN (mn) = trapz(sin(lat(don)*pi/180), squeeze(TSZ(mn,don))')/trapz(sin(lat(don)*pi/180), ones(size(squeeze(TSZ(mn,don))'))) - GTS;
    TST (mn) = trapz(sin(lat(dot)*pi/180), squeeze(TSZ(mn,dot))')/trapz(sin(lat(dot)*pi/180), ones(size(squeeze(TSZ(mn,dot))'))) - GTS;
    TSS (mn) = trapz(sin(lat(dos)*pi/180), squeeze(TSZ(mn,dos))')/trapz(sin(lat(dos)*pi/180), ones(size(squeeze(TSZ(mn,dos))'))) - GTS;
end

save('3_box_temps_tropmat', 'TVN', 'TVT', 'TVS', 'TSN', 'TST', 'TSS', 'GTV', 'GTS')

