%% Script to Execute Climatology Calculations
%This script inputs the buoyancy and stress potential vorticity (PV) flux 
%data sets, as well as, the sst data for the North Atlantic and uses those 
%sets to compute several climatology values (daily, annual, winter,
%summer).

%%Variables
%latitudeDomain: latitude range of the domain
%longitudeDomain: longitude range of the domain
%pvBuoyancy: buoyancy pv flux data for the North Atlantic, 20 year data
%   from 1988-2007
%pvStress: stress pv flux data for the North Atlantic, 20 year data from
%   1988-2007
%tst: sst data for the North Atlantic, 20 year data from 1988-2007
%pvBuoyancyMask: masked, Eighteen Degree Water (EDW) buoyancy pv flux data 
%for North Atlantic, 20 year data set
%pvStressMask: masked, EDW stress pv flux data for the North Atlantic, 20
%   year data set
%pvBuoyancyClimatologyDaily: daily buoyancy flux climatology of the edw for
%   the North Atlantic, array is (365,88,91)
%pvStressClimatologyDaily:
%pvBuoyancyClimatologyAnnual: annual bouyancy flux climatology of the edw
%for the North Atlantic, array is (88,91)
%pvStressClimatologyAnnual: annual stress flux climatology of the edw
%for the North Atlantic, array is (88,91)
%pvBuoyancyClimatologyDJF: December-January-February (DJF) buoyancy pv flux 
%   climatology of the EDW for the North Atlantic
%pvBuoyancyClimatologyJJA: June-July-August (JJA) buoyancy pv flux 
%   climatology of the EDW for the North Atlantic
%pvStressClimatologyDJF: DJF stress pv flux climatology of the EDW for the 
%   North Atlantic
%pvStressClimatologyJJA: JJA stress pv flux climatology of the EDW for the 
%   North Atlantic

%%Initial Values
%Sets the initial values for latitude and longitude, while also creating
%the initial vectors for PV due to buoyancy and stress
latitudeDomain=93:180;
longitudeDomain=270:360;
pvBuoyancy=zeros(20,365,length(latitudeDomain),length(longitudeDomain));
pvStress=zeros(20,365,length(latitudeDomain),length(longitudeDomain));
tst=zeros(20,365,length(latitudeDomain),length(longitudeDomain));

%%Data Insertion
%Inputs the appropriate 20 year data for the PV due to buoyancy and stress
for k=1988:2007;
fid=fopen(['qnet/pvbf2_' int2str(k)],'r','b');
pvbf=fread(fid,'real*8');pvbf=reshape(pvbf,365,180,360);
pvBuoyancy(k-1987,:,:,:)=pvbf(:,latitudeDomain,longitudeDomain);
fid=fopen(['ts/pvst_' int2str(k)],'r','b');   
pvws=fread(fid,'real*8');pvws=reshape(pvws,365,180,360);
pvStress(k-1987,:,:,:)=pvws(:,latitudeDomain,longitudeDomain);
filename=['ts/ts_oaflux_' int2str(k) '.nc'];
nc=netcdf(filename);
ts=nc{'tmpsf'}(:);
tst(k-1987,:,:,:)=ts(1:365,latitudeDomain,longitudeDomain)*.01;
k,end;

%%Data Masking
%Masks the data for PV flux in the North Atlantic due to buoyancy and stress
latitudeDomain=93:180;
longitudeDomain=270:360;
pvbfmask=zeros(20,365,length(latitudeDomain),length(longitudeDomain));
pvwsmask=zeros(20,365,length(latitudeDomain),length(longitudeDomain));
for k=1:20;
for d=1:365;
mask=ones(length(latitudeDomain),length(longitudeDomain));
th=squeeze(tst(k,d,:,:));
mask(th>19)=0;mask(th<17)=0;
pvBuoyancyMask(k,d,:,:)=(squeeze(pvBuoyancy(k,d,:,:)).*mask);
pvStressMask(k,d,:,:)=(squeeze(pvStress(k,d,:,:)).*mask);
end;k,end;

%%Daily Climatologies
%Calculates the daily buoyancy and stress PV flux climatologies of the EDW 
%for the North Atlantic (array is 365, 88, 91)
pvBuoyancyClimatologyDaily=squeeze(sum(pvBuoyancyMask,1))/20*2e16;
pvStressClimatologyDaily=squeeze(sum(pvStressMask,1))/20*2e16;

%Annual Climatologies
%Calculates the annual bouyancy and stress PV flux climatology of the EDW 
%for the North Atlantic (array is 88, 91)
pvBuoyancyClimatologyAnnual=squeeze(sum(pvBuoyancyClimatologyDaily,1))/365;
pvStressClimatologyAnnual=squeeze(sum(pvStressClimatologyDaily,1))/365;

%%Winter Climatologies
%Calculates the DJF buoyancy and stress PV flux climatologies of the EDW 
%for the North Atlantic
pvBuoyancyClimatologyDJF=squeeze(sum(pvBuoyancyClimatologyDaily(335:365,:,:),1));
pvBuoyancyClimatologyDJF=pvBuoyancyClimatologyDJF+squeeze(sum(pvBuoyancyClimatologyDaily(1:59,:,:),1));
pvBuoyancyClimatologyDJF=pvBuoyancyClimatologyDJF/(length(335:365)+length(1:59));

pvStressClimatologyDJF=squeeze(sum(pvStressClimatologyDaily(335:365,:,:),1));
pvStressClimatologyDJF=pvStressClimatologyDJF+squeeze(sum(pvStressClimatologyDaily(1:59,:,:),1));
pvStressClimatologyDJF=pvStressClimatologyDJF/(length(335:365)+length(1:59));

%%Summer Climatologies
%Calculates the JJA buoyancy and stress PV flux climatologies of the EDW 
%for the North Atlantic
pvBuoyancyClimatologyJJA=squeeze(sum(pvBuoyancyClimatologyDaily(152:243,:,:),1));
pvBuoyancyClimatologyJJA=pvBuoyancyClimatologyJJA/(length(152:243));

pvStressClimatologyJJA=squeeze(sum(pvStressClimatologyDaily(152:243,:,:),1));
pvStressClimatologyJJA=pvStressClimatologyJJA/(length(152:243));

