clear all;

%The domains used are   Y: 0N - 60N, grid pts every 1 deg  
%                       X: 20W- 50W, grid pts every 1 deg
%                       Z: 0-3000m , grid pts every 10m

%load annual and seasonal temperature and salinity
%the seasonal fields have the matrix form T(Y,X,Z,t) where t=1,2,3,4 depending on
%the season
%the annual fields have the matrix form T(Y,X,Z) 

load X;load Y;load Z;load tempa;load sala;load temps;load sals; 

%extrapolate temp and sal to get rid of the NaN values assigned to the
%missing data points

warning('off');
for i=1:length(X);for j=1:length(Z);tempa_new(:,i,j)=spline(Y,tempa(:,i,j),Y);end;end;
for i=1:length(X);for j=1:length(Z);sala_new(:,i,j)=spline(Y,sala(:,i,j),Y);end;end;
for i=1:length(X);for j=1:length(Z);for it=1:4;temps_new(:,i,j,it)=spline(Y,temps(:,i,j,it),Y);end;end;end;
for i=1:length(X);for j=1:length(Z);for it=1:4;sals_new(:,i,j,it)=spline(Y,sals(:,i,j,it),Y);end;end;end;

%take the zonal average
for i=1:length(Y);for j=1:length(Z);tempa_xav(i,j)=mean(tempa_new(i,:,j));end;end;
for i=1:length(Y);for j=1:length(Z);sala_xav(i,j)=mean(sala_new(i,:,j));end;end;
for i=1:length(Y);for j=1:length(Z);for it=1:4;temps_xav(i,j,it)=mean(temps_new(i,:,j,it));end;end;end;
for i=1:length(Y);for j=1:length(Z);for it=1:4;sals_xav(i,j,it)=mean(sals_new(i,:,j,it));end;end;end;

%define the winter and summer mean
for i=1:length(Y);temp_xav_w(i,:)=temps_xav(i,:,1);end;
for i=1:length(Y);temp_xav_s(i,:)=temps_xav(i,:,3);end;
for i=1:length(Y);sal_xav_w(i,:)=sals_xav(i,:,1);end;
for i=1:length(Y);sal_xav_s(i,:)=sals_xav(i,:,3);end;

%calculate rho using the equations in Gill's Appendix
rho_w=myfunct(temp_xav_w,sal_xav_w,Z)-1000;
rho_s=myfunct(temp_xav_s,sal_xav_s,Z)-1000;
rho_a=myfunct(tempa_xav,sala_xav,Z)-1000;

%plot T,S,rho for annual, winter and summer
figure;subplot(211);h=contourf(Y,Z,tempa_xav.');clabel(h);set(gca,'YDir','reverse');xlabel('y(deg)');ylabel('z (m)');title('T (^o C) - annual mean');
subplot(223);h=contourf(Y,Z,temp_xav_w.');clabel(h);set(gca,'YDir','reverse');xlabel('y(deg)');ylabel('z (m)');title('T (^o C) - winter');
subplot(224);h=contourf(Y,Z,temp_xav_s.');clabel(h);set(gca,'YDir','reverse');xlabel('y(deg)');ylabel('z (m)');title('T (^o C) - summer');

figure;subplot(211);h=contourf(Y,Z,sala_xav.');clabel(h);set(gca,'YDir','reverse');xlabel('y(deg)');ylabel('z (m)');title('S (p.s.u) - annual mean');
subplot(223);h=contourf(Y,Z,sal_xav_w.');clabel(h);set(gca,'YDir','reverse');xlabel('y(deg)');ylabel('z (m)');title('S (p.s.u) - winter');
subplot(224);h=contourf(Y,Z,sal_xav_s.');clabel(h);set(gca,'YDir','reverse');xlabel('y(deg)');ylabel('z (m)');title('S (p.s.u) - summer');

figure;subplot(211);h=contourf(Y,Z,rho_a.');clabel(h);set(gca,'YDir','reverse');xlabel('y(deg)');ylabel('z (m)');title('\rho (kgr/m^3-1000) - annual mean');
subplot(223);h=contourf(Y,Z,rho_w.');clabel(h);set(gca,'YDir','reverse');xlabel('y(deg)');ylabel('z (m)');title('\rho (kgr/m^3-1000) - winter');
subplot(224);h=contourf(Y,Z,rho_s.');clabel(h);set(gca,'YDir','reverse');xlabel('y(deg)');ylabel('z (m)');title('\rho (kgr/m^3-1000) - summer');

%Divide into four boxes : box 1: 0N  - 50N + 0-500m 
%                         box 2: 50N - 60N + 0-500m 
%                         box 3: 50N - 60N + 500-3000m 
%                         box 4: 0N  - 50N + 500-3000m 
%Calculate the mean T,S for summer, winter and annual for the four boxes


t1w=mean(mean(temp_xav_w(1:50,1:14)))
t2w=mean(mean(temp_xav_w(51:60,1:14)))
t3w=mean(mean(temp_xav_w(51:60,15:28)))
t4w=mean(mean(temp_xav_w(1:50,15:28)))

t1s=mean(mean(temp_xav_s(1:50,1:14)))
t2s=mean(mean(temp_xav_s(51:60,1:14)))
t3s=mean(mean(temp_xav_s(51:60,15:28)))
t4s=mean(mean(temp_xav_s(1:50,15:28)))

t1a=mean(mean(tempa_xav(1:50,1:14)))
t2a=mean(mean(tempa_xav(51:60,1:14)))
t3a=mean(mean(tempa_xav(51:60,15:28)))
t4a=mean(mean(tempa_xav(1:50,15:28)))

s1w=mean(mean(sal_xav_w(1:50,1:14)))
s2w=mean(mean(sal_xav_w(51:60,1:14)))
s3w=mean(mean(sal_xav_w(51:60,15:28)))
s4w=mean(mean(sal_xav_w(1:50,15:28)))

s1s=mean(mean(sal_xav_s(1:50,1:14)))
s2s=mean(mean(sal_xav_s(51:60,1:14)))
s3s=mean(mean(sal_xav_s(51:60,15:28)))
s4s=mean(mean(sal_xav_s(1:50,15:28)))

s1a=mean(mean(sala_xav(1:50,1:14)))
s2a=mean(mean(sala_xav(51:60,1:14)))
s3a=mean(mean(sala_xav(51:60,15:28)))
s4a=mean(mean(sala_xav(1:50,15:28)))

