%compute the vertically integrated mse for each month
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=1984;
endyear=1990;
for yr=startyear:endyear;
dir = '/home/disk/eos10/dargan/aquafv/output/seas12mml/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);
plev = nc{'pfull'}(:); NL=length(plev);
pb=[0 mean([plev(2:end) plev(1:end-1)]') 1000]; %boundaries of each pressure level
w = (pb(2:end)-pb(1:end-1))/(max(pb)-min(pb))'; %mass weighting vector for vertical integrals

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  TOA %%%%%%%%%%%%%%%%%%%%%%%%%%
data(1:12, 1:NL,1:NN) = squeeze(mean(shiftdim(nc{'mse'}(:,:,:),3))); %zonal mean mse
data2(1:12, 1:NL,1:NN) = squeeze(mean(shiftdim(nc{'dse'}(:,:,:),3))); %zonal mean dse

 for mn=1:12
     d=squeeze(data(mn,:,:));
     d2=squeeze(data2(mn,:,:));
     MSE(yr-startyear+1,mn,1:NN) = (w * d) * (10^5/9.81); %vertically integrated MSE in J/m^2;
     DSE(yr-startyear+1,mn,1:NN) = (w * d2) * (10^5/9.81); %vertically integrated MSE in J/m^2;
     %GMSE(yr-startyear+1,mn)=2*pi*(a^2)*trapz(sin(lat*pi/180), squeeze(MSE(mn,:))); %global mean MSE in J
     %GDSE(yr-startyear+1,mn)=2*pi*(a^2)*trapz(sin(lat*pi/180), squeeze(DSE(mn,:))); %global mean MSE in J
 end

 
 
%  figure(3); clf;
%  plot(lat, MSE);
%  legend('1','2','3','4','5','6','7','8','9','10','11','12')

end
MSE=squeeze(mean(GSE));
DSE=squeeze(mean(DSE));
 

MSE_TEN=(MSE([2:12 1],:) - MSE([12 1:11],:))/ ( 365/12 * 24 * 3600); %central difference
DSE_TEN=(DSE([2:12 1],:) - DSE([12 1:11],:))/ ( 365/12 * 24 * 3600); %central difference

save('MSE_TEN', 'lat','MSE_TEN', 'DSE_TEN')