% 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/seas/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));

%define the mass weigthing vector (normalized)
trop_start=7;
W=(Phalf(2:end) - Phalf(1:end-1)) / (max(Phalf) - min(Phalf));
W_trop=zeros(size(W)); W_trop(7:end) = (Phalf(trop_start+1:end) - Phalf(trop_start:end-1)) / (max(Phalf(trop_start:end)) - min(Phalf(trop_start:end)));
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
            TV_trop(nt,nn,n) = squeeze(T(nt,:,nn,n)) * W_trop; 
            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


save('TV_trop', 'TS', 'TV', 'TV_trop', 'T9', 'T8', 'lat', 'lon')
