%% Make a readable matrix out of the .bin file 
clear;
prec='real*4';
ieee='b';
fid=fopen('bathy_2deg_nowall.bin','r',ieee); 
read_bat=fread(fid,prec); 
fclose(fid);
bathy = reshape(read_bat,[180 90]);

nx = 180;
ny = 90;

lon = linspace(0,358,nx);
lat = linspace(-90,88,ny);

%% If you already have som readable output 
%h = rdmds('../Depth/Depth');
%lon_c = rdmds('../Grid/XC');
%lat_c = rdmds('../Grid/YC');
%height_c = rdmds('../Grid/RC');
%height_f = rdmds('../Grid/RF');
%lon_g = rdmds('../Grid/XG');
%lat_g = rdmds('../Grid/YG');

%% Smoothing command

times_smooth = 2;
h_smooth = smooth9(bathy,times_smooth);

[r,c]= find(bathy>=0);
h_smooth2 = h_smooth;
for i = 1:length(r)
    h_smooth2(r(i), c(i)) = 0;
end

[x,y] = meshgrid(lon,lat);
figure;
pcolor(x,y,bathy')
%shading('interp')
colorbar

%figure;
%pcolor(x,y,h_smooth')
%shading('interp')
%colorbar

figure;
pcolor(x,y,h_smooth2')
%shading('interp')
colorbar
% 
% V = [0 0];
% figure;
% contourf(x, y,bathy',V,'color','k')
% figure;
% contourf(x, y,h_smooth',V,'color','k')
% figure;
% contourf(x, y,h_smooth2',V,'color','k')
%% Plot topo at one latitude

n = 83;                                        %# Which latitude to plot 

figure;
h_max = min(bathy(:,n));
x = lon(1,:);                                %#initialize x array
y1 = bathy(:,n)';                                  %#create first curve
y2 = (h_max-500)*ones(1,length(lon(1,:)));   %#create second curve
X=[x,fliplr(x)];                                %#create continuous x value array for plotting
Y=[y1,fliplr(y2)];                              %#create y values for out and then back
fill(X,Y,'k')                                   %#plot filled area
xlim([0 358])
ylim([(h_max-500) 0])
% 
% figure;
% h_max = min(h_smooth(:,n));
% x = lon(1,:);                                %#initialize x array
% y1 = h_smooth(:,n)';                           %#create first curve
% y2 = (h_max-500)*ones(1,length(lon(1,:)));   %#create second curve
% X=[x,fliplr(x)];                                %#create continuous x value array for plotting
% Y=[y1,fliplr(y2)];                              %#create y values for out and then back
% fill(X,Y,'k')                                   %#plot filled area
% xlim([0 358])
% ylim([(h_max-500) 0])

figure;
h_max = min(h_smooth2(:,n));
x = lon(1,:);                                %#initialize x array
y1 = h_smooth2(:,n)';                           %#create first curve
y2 = (h_max-500)*ones(1,length(lon(1,:)));   %#create second curve
X=[x,fliplr(x)];                                %#create continuous x value array for plotting
Y=[y1,fliplr(y2)];                              %#create y values for out and then back
fill(X,Y,'k')                                   %#plot filled area
xlim([0 358])
ylim([(h_max-500) 0])

%% Make a smooth topography file
%fid = fopen('smooth_bathy_pangea_2deg.bin', 'w','b');
%fwrite(fid, h_smooth, prec);
%fclose(fid);

% Same land points as non-smooth topography file
fid = fopen('smooth_bathy_2deg_nowall.bin', 'w','b');
fwrite(fid, h_smooth2, prec);
fclose(fid);

%% To check the written temp-file
clear;
prec='real*4';
ieee='b';
%lon = rdmds('../../Grid82/XG');
%lat = rdmds('../../Grid82/YG');
nx = 180;
ny = 90;

lon = linspace(0,358,nx);
lat = linspace(-90,88,ny);
[x,y] = meshgrid(lon,lat);

fid=fopen('smooth_bathy_2deg_nowall.bin','r',ieee); 
smooth_bathy=fread(fid,prec); 
fclose(fid);
bathy2 = reshape(smooth_bathy,[180 90]);

%plot_T = squeeze(temp2(:,:,10));

%V = [0 0];
figure
%contourf(lon,lat,plot_T)
pcolor(x,y,bathy2')
%view(0,90)
shading('interp')
colorbar