%%Monthly PV Values
%This script calculates the monthly integrated PV values that are due to
%either buoyancy or stress

%%Variables
%latitude: the latitude
%latitudeDomain: latitude range of the domain
%longitudeDomain: longitude range of the domain
%northAtlanticOutcropArea: area of the 17-19 outcrop in the North Atlantic
%   for the 20 year data set
%northAtlanticOutcropAreaAvg: average area of the 17-19 outcrop in the
%   North Atlantic for the 20 year data set
%pvBuoyancyClimatologyMonthly: monthly distribution climatology of buoyancy
%   PV flux (in W/m^2)
%pvStressClimatologyMonthly: monthly distribution climatology of stress
%   PV flux (in W/m^2)
%pvBuoyancyMonthlyIntegrated0090: integrated buoyancy PV flux from 0 to 90W
%pvBuoyancyMonthlyIntegrated0030: integrated buoyancy PV flux from 0 to 30W
%pvBuoyancyMonthlyIntegrated3060: integrated buoyancy PV flux from 30 to 60W
%pvBuoyancyMonthlyIntegrated6090: integrated buoyancy PV flux from 60 to 90W
%pvStressMonthlyIntegrated0090: integrated stress PV flux from 0 to 90W
%pvStressMonthlyIntegrated0030: integrated stress PV flux from 0 to 30W
%pvStressMonthlyIntegrated3060: integrated stress PV flux from 30 to 60W
%pvStressMonthlyIntegrated6090: integrated stress PV flux from 60 to 90W

%%Load the PV flux climatology data 
load pvfluxna_new

%%Initial Values
%Creates the initial monthly climatology and monthly PV integrated data
%vectors
pvBuoyancyClimatologyMonthly=zeros(12,length(latitudeDomain),length(longitudeDomain));
pvStressClimatologyMonthly=zeros(12,length(latitudeDomain),length(longitudeDomain));

northAtlanticOutcropArea=[1:7300]*0;;

pvBuoyancyMonthlyIntegrated0090=(1:12)*0;
pvBuoyancyMonthlyIntegrated0030=pvBuoyancyMonthlyIntegrated0090;
pvBuoyancyMonthlyIntegrated3060=pvBuoyancyMonthlyIntegrated0090;
pvBuoyancyMonthlyIntegrated6090=pvBuoyancyMonthlyIntegrated0090;
pvStressMonthlyIntegrated0090=pvBuoyancyMonthlyIntegrated0090;
pvStressMonthlyIntegrated0030=pvBuoyancyMonthlyIntegrated0090;
pvStressMonthlyIntegrated3060=pvBuoyancyMonthlyIntegrated0090;
pvStressMonthlyIntegrated6090=pvBuoyancyMonthlyIntegrated0090;

%%Calculates the area of the 17-19 outcrop in the North Atlantic
jd=[1 32 60 91 121 152 182 213 244 274 305 335 366];
dy=6370*1.e3*pi/180;
latitude=-89.5:89.5;
dx=dy*cosd(latitude(latitudeDomain));
northAtlanticArea=dx*dy;
	for k=1988:2007;
		for day=1:365;
			mask=ones(length(latitudeDomain),length(longitudeDomain));
			th=squeeze(tst(k-1987,day,:,:));
			mask(th>19)=0;mask(th<17)=0;
			northAtlanticOutcropArea((k-1988)*365+day)=squeeze(sum(sum(mask,2).*northAtlanticArea'));
		end;
	end;
    
%%Calculates the monthly integrated buoyancy and stress PV flux values
	for mon=1:12;
		pvBuoyancyMonthly2d(mon,:,:)=squeeze(nansum(nansum(pvBuoyancyMask(:,jd(mon):jd(mon+1)-1,:,:),1),2))/length(jd(mon):jd(mon+1)-1)/20*2e16;
		pvStressMonthly2d(mon,:,:)=squeeze(nansum(nansum(pvStressMask(:,jd(mon):jd(mon+1)-1,:,:),1),2))/length(jd(mon):jd(mon+1)-1)/20*2e16;
		pvBuoyancyMonthlyIntegrated0090(mon)=sum(squeeze(sum(pvBuoyancyMonthly2d(mon,:,:),3)).*northAtlanticArea,2);
		pvStressMonthlyIntegrated0090(mon)=sum(squeeze(sum(pvStressMonthly2d(mon,:,:),3)).*northAtlanticArea,2);
		pvBuoyancyMonthlyIntegrated0030(mon)=sum(squeeze(sum(pvBuoyancyMonthly2d(mon,:,1:30),3)).*northAtlanticArea,2);
		pvStressMonthlyIntegrated0030(mon)=sum(squeeze(sum(pvStressMonthly2d(mon,:,1:30),3)).*northAtlanticArea,2);
		pvBuoyancyMonthlyIntegrated3060(mon)=sum(squeeze(sum(pvBuoyancyMonthly2d(mon,:,31:60),3)).*northAtlanticArea,2);
		pvStressMonthlyIntegrated3060(mon)=sum(squeeze(sum(pvStressMonthly2d(mon,:,31:60),3)).*northAtlanticArea,2);
		pvBuoyancyMonthlyIntegrated6090(mon)=sum(squeeze(sum(pvBuoyancyMonthly2d(mon,:,61:90),3)).*northAtlanticArea,2);
		pvStressMonthlyIntegrated6090(mon)=sum(squeeze(sum(pvStressMonthly2d(mon,:,61:90),3)).*northAtlanticArea,2);
	mon,end;

	northAtlanticOutcropAreaAnnualAverage=sum(northAtlanticOutcropArea)/length(northAtlanticOutcropArea);
	pvBuoyancyMonthlyIntegrated0090=pvBuoyancyMonthlyIntegrated0090/northAtlanticOutcropAreaAnnualAverage;
	pvStressMonthlyIntegrated0090=pvStressMonthlyIntegrated0090/northAtlanticOutcropAreaAnnualAverage;
	pvBuoyancyMonthlyIntegrated0030=pvBuoyancyMonthlyIntegrated0030/northAtlanticOutcropAreaAnnualAverage;
	pvStressMonthlyIntegrated0030=pvStressMonthlyIntegrated0030/northAtlanticOutcropAreaAnnualAverage;
	pvBuoyancyMonthlyIntegrated3060=pvBuoyancyMonthlyIntegrated3060/northAtlanticOutcropAreaAnnualAverage;
	pvStressMonthlyIntegrated3060=pvStressMonthlyIntegrated3060/northAtlanticOutcropAreaAnnualAverage;
	pvBuoyancyMonthlyIntegrated6090=pvBuoyancyMonthlyIntegrated6090/northAtlanticOutcropAreaAnnualAverage;
	pvBuoyancyMonthlyIntegrated6090=pvBuoyancyMonthlyIntegrated6090/northAtlanticOutcropAreaAnnualAverage;

