%%Script to compute Global Values
%This script will compute the globally integrated annual averaged buoyant
%and stress pv flux, as well as, a rough value for the ocean surface area. 

%%Variables:
%conversionFactor: 2e16, conversion factor from pvflux to W/m^2
%latitude: the latitude
%northAtlanticArea: area for the North Atlantic
%latitudeDomain: latitude range of the domain
%longitudeDomain: longitude range of the domain
%pvBuoyancyGlobalAverage: time averaged buoyant pv flux
%pvStressGlobalAverage: time averaged stress pv flux
%pvBuoyancyGlobalTime: globally integrated annual averaged buoyant pv flux
%pvStressGlobalTime: globally integrated annual averaged stress pv flux
%globalArea: Roughly the ocean surface area


%%Uses the latitude to calculate the area for the North Atlantic
conversionFactor=4.e3*1.e3/2.e-4/1.e-4*1.e2;
dy=6370*pi/180*1.e3;
latitude=-89.5:89.5;
dx=dy*cosd(latitude);
northAtlanticArea=dx*dy;

%%Uses the inputted data to calculate the time averaged buoyant and stress
%%pv flux
latitudeDomain=1:180;
longitudeDomain=1:360;
pvBuoyancyGlobalAverage=zeros(length(latitudeDomain),length(longitudeDomain));
pvStressGlobalAverage=zeros(length(latitudeDomain),length(longitudeDomain));

for k=1988:2007;
fid=fopen(['qnet/pvbf2_' int2str(k)],'r','b');
pvBuoyancy=fread(fid,'real*8');pvbf=reshape(pvbf,365,180,360);
pvBuoyancyGlobalAverage=pvBuoyancyGlobalAverage+squeeze(nansum(pvBuoyancy,1));
temp=nansum(squeeze(nansum(pvBuoyancy,1)),2);
pvBuoyancyGlobalTime(k-1987,1)=nansum(temp(1:90).*northAtlanticArea(1:90)');
pvBuoyancyGlobalTime(k-1987,2)=nansum(temp(91:180).*northAtlanticArea(91:180)');
fid=fopen(['ts/pvst_' int2str(k)],'r','b');   
pvStress=fread(fid,'real*8');pvStress=reshape(pvStress,365,180,360);
pvStressGlobalAverage=pvStressGlobalAverage+squeeze(nansum(pvStress,1));
temp=nansum(squeeze(nansum(pvStress,1)),2);
pvStressGlobalTime(k-1987,1)=nansum(temp(1:90).*northAtlanticArea(1:90)');
pvStressGlobalTime(k-1987,2)=nansum(temp(91:180).*northAtlanticArea(91:180)');
k,end;

%%Uses the conversion factor to convert the values from pv flux to W/m^2
pvStressGlobalAverage=pvStressGlobalAverage/7300*conversionFactor;
pvBuoyancyGlobalAverage=pvBuoyancyGlobalAverage/7300*conversionFactor;
pvBuoyancyGlobalTime=pvBuoyancyGlobalTime/365*conversionFactor;
pvStressGlobalTime=pvStressGlobalTime/365*conversionFactor;

%%Calculates the ocean surface area
globalArea=360*nansum(northAtlanticArea)*.7;



