%%Climatological PV Flux Calculations
%This script calculates the climatological buoyant & stress PV flux based
%on NMLD mixed layer depth climatology

%%Variables
%latitudeDomain: latitude range of the domain
%longitudeDomain: longitude range of the domain
%conversionFactor: 2e16, conversion factor from PV flux to W/m^2
%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 1998-2007
%pvBuoyancyIntegrated: 20 years daily intefrated PV flux into the EDW due
%   to buoyancy
%pvStressIntegrated: 20 years daily intefrated PV flux into the EDW due to
%	stress
%pvBuoyancyMasked: masked, EDW buoyancy PV flux data for the North Atlantic, 20
%   year data
%pvStressMasked: masked, EDW stress PV flux data for the North Atlantic, 20
%   year data
%pvBuoyancyMixedLayer: climatological buoyant PV flux based on NMLD mixed layer depth
%   climatology
%pvStressMixedLayer: climatological stress PV flux based on NMLD mixed layer depth
%   climatology

%%Sets the given values for the latitude, longitude, & conversion factor
latitudeDomain=93:180;
longitudeDomain=270:360;
conversionFactor=4.e3*1.e3/2.e-4/1.e-4*1.e2;

%Creates the buoyancy & stress PV flux vectors
pvBuoyancy=zeros(20,365,length(latitudeDomain),length(longitudeDomain));
pvStress=zeros(20,365,length(latitudeDomain),length(longitudeDomain));
tst=zeros(20,365,length(latitudeDomain),length(longitudeDomain));

%%Fills in the data for the PV due to buoyancy & stress PV flux vectors
for k=1988:2007;
fid=fopen(['qnet/pvBuoyancyMixedLayer_' 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;

%%Creates the vectors for the integrated PV flux due to buoyancy & stress
latitudeDomain=93:180;
longitudeDomain=270:360;
dy=6370*pi/180*1.e3;
latitude=-89.5:89.5;
dx=dy*cosd(lat(latitudeDomain));
northAtlanticarea=dx*dy;
pvBuoyancyIntegrated=[1:365]*0;
pvStressIntegrated=pvBuoyancyIntegrated;

%%Inputs the data into the PV integrated vectors due to buoyancy & stress
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;
pvBuoyancyMasked(:,:)=(squeeze(pvBuoyancy(k,d,:,:)).*mask);
pvStressMasked(:,:)=(squeeze(pvStress(k,d,:,:)).*mask);
pvBuoyancyIntegrated(d)=pvBuoyancyIntegrated(d)+nansum(squeeze(nansum(pvBuoyancyMasked,2).*northAtlanticArea'));
pvStressIntegrated(d)=pvStressIntegrated(d)+nansum(squeeze(nansum(pvStressMasked,2).*northAtlanticArea'));
end;k,end;

%%Calculates the climatological PV flux based on NMLD mixed layer depth
%%climatology
pvBuoyancyMixedLayer=pvBuoyancyIntegrated/20*2.e16;
pvStressMixedLayer=pvStressIntegrated/20*2.e16;



