%% Read bathymetry.bin file
% F_topo used to set land area to zero in the forcing file regardless of
% the topography used
clear;
prec='real*4';
ieee='b';
nx=180;
ny=90;
fid=fopen('bathy_2deg_nowall.bin','r',ieee); 
F_topo=fread(fid,[nx,ny],prec); 
fclose(fid);

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

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

figure
%contourf(x,y,F_topo')
pcolor(x,y,F_topo')
%view(0,90)
shading('interp')
colorbar
xlabel('Longitude')
ylabel('Latitude')

%F_topo(:,1:find(lat==-84))=0;
%F_topo(:,find(lat==84):end)=0;

%fid = fopen('bathy_2deg_pangea24z.bin', 'w','b');
%fwrite(fid, F_topo, prec);
%fclose(fid);

%% zonally symmetric forcing

%lat=linspace(-82,81,ny);
%lon=linspace(0,359,nx);
ii=find(F_topo>=0);
% % The net heat flux forcing
% q1=200.*cosd(2*lat);
% figure
% plot(lat,q1)
% Q=ones(nx,ny).*repmat(q1,[nx 1]);
% Q(ii)=0;

% The SST clim file
sst1=25*cosd(lat);
figure
plot(lat,sst1)
xlim([-90 88])
SST=ones(nx,ny).*repmat(sst1,[nx 1]);
SST(ii)=0;

% % SW heat flux at surface
% qsw=300.*cosd(1.07*lat);
% figure
% plot(lat,qsw)
% Qsw=ones(nx,ny).*repmat(qsw,[nx 1]);
% Qsw(ii)=0;

%% plot results
% figure
% surf(x,y,Q')
% view(0,90)
% shading('interp')
% xlim([0 360])
% ylim([-90 90])
% colorbar
% 
figure
surf(x,y,SST')
view(0,90)
shading('interp')
xlim([0 360])
ylim([-90 88])
colorbar
xlabel('Longitude')
ylabel('Latitude')
title('\bf Sea Surface Temperature Forcing')
caxis([0 25])

% figure
% surf(x,y,Qsw')
% view(0,90)
% shading('interp')
% xlim([0 360])
% ylim([-90 90])
% colorbar
% xlabel('Longitude')
% ylabel('Latitude')
% title('\bf Shorwave radition forcing at the sea surface')
% %caxis([0 400])


%% Write .bin files
% fid = fopen('Input_files/zonal_qnet.bin', 'w','b');
% fwrite(fid, Q, prec);
% fclose(fid);

 fid = fopen('zonal_sst_2deg_nowall.bin', 'w','b');
 fwrite(fid, SST, prec);
 fclose(fid);

% fid = fopen('Input_files/zonal_qsw.bin', 'w','b');
% fwrite(fid, Qsw, prec);
% fclose(fid);

%% To check the written files
clear;
prec='real*4';
ieee='b';
nx=180;
ny=90;
fid=fopen('bathy_2deg_nowall.bin','r',ieee); 
%fid=fopen('zonal_sst_2deg_nowall.bin','r',ieee); 
F=fread(fid,[nx,ny],prec); 
fclose(fid);

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

figure
pcolor(x,y,F')
%view(0,90)
shading('interp')
colorbar

pcolor(F')
%view(0,90)
shading('interp')
colorbar



%hydrogThetaFile='init_temp.bin'