% Load grids
xcs=rdmds(fullfile(inputDirPath,'XC'));
ycs=rdmds(fullfile(inputDirPath,'YC')); 
xcg=rdmds(fullfile(inputDirPath,'XG'));
ycg=rdmds(fullfile(inputDirPath,'YG'));
x=rdmds(fullfile(inputDirPath,'XC'));
y=rdmds(fullfile(inputDirPath,'YC'));

% Load data
t=rdmds(fullfile(inputDirPath,'Ttave'),NaN);
s=rdmds(fullfile(inputDirPath,'Stave'),NaN);

% Plot temperature
figure(1)
shift=-1;
c1=-2;
c2=30;

grph_CS(sq(t(:,:,1)),xcs,ycs,xcg,ycg,c1,c2,shift)
title('Annual average surface temperature (C)');

% Plot salinity
figure(2)
shift=-1;
c1=25;
c2=40;

grph_CS(sq(s(:,:,1)),xcs,ycs,xcg,ycg,c1,c2,shift)
title('Annual average surface salinity ');


%=======================  EH3  =======================

clear all, close all
tx = 85;
ty = 85;
nt = 216;
cx = 510;
cy = 510;
nz = 1;
%  XCS
phi = readbin('XC.data',[tx*nt,ty,nz]);
foo
%%for i = 1:6
%%  surf(te(:,:,i)), view(2), shading interp
%%  pause(1)
%%end
xcs = zeros(6*510,510);
for i = 1:6
  xb = (i-1)*510 + 1;  xe = xb + 510 - 1;
  yb = 1;              ye = 510;
  xcs(xb:xe,yb:ye) = te(:,:,i);
end
xcs = xcs + -360.0*(xcs > 180.0);
clear phi te
%%surf(xcs(:,:)'), view(2), axis equal, shading interp

%  YCS
phi = readbin('YC.data',[tx*nt,ty,nz]);
foo
ycs = zeros(6*510,510);
for i = 1:6
  xb = (i-1)*510 + 1;  xe = xb + 510 - 1;
  yb = 1;              ye = 510;
  ycs(xb:xe,yb:ye) = te(:,:,i);
end
%%ycs = ycs + -360.0*(ycs > 180.0);
%%surf(ycs(:,:)'), view(2), axis equal, shading interp

%  XGS
phi = readbin('XG.data',[tx*nt,ty,nz]);
foo
xgs = zeros(6*510,510);
for i = 1:6
  xb = (i-1)*510 + 1;  xe = xb + 510 - 1;
  yb = 1;              ye = 510;
  xgs(xb:xe,yb:ye) = te(:,:,i);
end
xgs = xgs + -360.0*(xgs > 180.0);
%%surf(xgs(:,:)'), view(2), axis equal, shading interp

%  YGS
phi = readbin('YG.data',[tx*nt,ty,nz]);
foo
ygs = zeros(6*510,510);
for i = 1:6
  xb = (i-1)*510 + 1;  xe = xb + 510 - 1;
  yb = 1;              ye = 510;
  ygs(xb:xe,yb:ye) = te(:,:,i);
end
%%ygs = ygs + -360.0*(ygs > 180.0);
%%surf(ygs(:,:)'), view(2), axis equal, shading interp
clear phi te

[ min(min(xcs)) max(max(xcs)) ; ...
  min(min(ycs)) max(max(ycs)) ; ...
  min(min(xgs)) max(max(xgs)) ; ...
  min(min(ygs)) max(max(ygs)) ; ...
]

%  ave-1995
%
phi = readbin('Tave-1995',[tx*nt,ty,nz]);
foo
%%for i = 1:6
%%  surf(te(:,:,i)), view(2), shading interp
%%  pause(2)
%%end
ta1995 = zeros(6*510,510);
for i = 1:6
  xb = (i-1)*510 + 1;  xe = xb + 510 - 1;
  yb = 1;              ye = 510;
  ta1995(xb:xe,yb:ye) = te(:,:,i);
end
clear phi te
%%surf(ta1995(:,:)'), view(2), axis equal, shading interp



%  Plot using "grph_CS"

shift=-1;
c1=0;
c2=0;
grph_CS(sq(ta1995),xcs,ycs,xgs,ygs,c1,c2,shift)
grph_CS(sq(ta1995),xcs,ycs,xgs,ygs,-2,32,shift)