%  The data is on a grid with 474 rows and 576 columns.  
%  The bottom left corner of the grid is 53.1358333 degrees lat and -3.99875 degrees long.  
%  The grid increments are 1/600 for lat and 1/400 for long.  
%  The format of the data is written out in is 20f9.2 (-999.99 is land).


datadir = fullfile('..', 'raw');
dir(datadir)
filename = fullfile(datadir, 'liv.bathymetry.dat')
[fid1,message]=fopen(filename,'r','n')% load POLCOMS bathymetry


%LB
 nlat = 474;
 nlon = 576;
 lonw= -3.99875; lone= -2.56125;
 lats= 53.1358333; latn= 53.92416663;


dx = (lone-lonw) / (nlon-1); %1/400
dy = (latn-lats) / (nlat-1); %1/600
xx= [lonw:dx:lone] ;
yy= [lats:dy:latn];

 
bath = zeros(1,nlon*nlat);

int = fix(nlon*nlat / 50) ;
for i=1:int
    i1 = (i-1)*50 + 1;
    i2 = i*50;
    bath(i1:i2) = fscanf(fid1,'%e',50);
end
i1 = i2 + 1;
i2 = ((198*224)-44350)+i1-1;
a = fix(i2-i1) + 1;
bath(i1:i2) = fscanf(fid1,'%e',a);
status = fclose(fid1)

bath1 = zeros(nlon,nlat);
for j=1:nlat
    i1 = (j-1)*nlon + 1;
    i2 = j*nlon;
    bath1(:,j) = bath(i1:i2);
end

data=flipud(-bath1');

%I = find(bath1==0); bath1(I) = NaN;
I = find(bath1==-9999.99); bath1(I) = NaN;
%I = find(bath1==-999); bath1(I) = NaN;

(fileparts(mfilename('fullpath')))
outputfile = 'liverpool_topobathy.nc';
nc_create_empty(outputfile);
nc_dump (outputfile);
nc_add_dimension(outputfile, 'latitude', nlat)
nc_add_dimension(outputfile, 'longitude', nlon)
nc_dump (outputfile);
s.Name      = 'lat';
s.Nctype    = nc_double;
s.Dimension = {'latitude'};
s.Attribute = struct('Name', {'standard_name', 'units', 'axis'}, 'Value', {'latitude', 'degree_north', 'X'});
nc_addvar(outputfile, s);

s.Name      = 'lon';
s.Nctype    = nc_double;
s.Dimension = {'longitude'};
s.Attribute = struct('Name', {'standard_name', 'units', 'axis'}, 'Value', {'longitude', 'degree_east', 'Y'});
nc_addvar(outputfile, s);
nc_dump (outputfile);
s.Name      = 'altitude';
s.Nctype    = nc_double;
s.Dimension = {'longitude', 'latitude'};
s.Attribute = struct('Name', {'standard_name', 'units', 'comment', 'coordinates', '_FillValue', 'positive'}, ...
    'Value', {'surface_altitude', 'm', 'altitude above MSL (probably)', 'lat lon', NaN, 'down'});
nc_addvar(outputfile, s);
nc_dump (outputfile);
nc_varput(outputfile, 'altitude', bath1)
nc_varput(outputfile, 'lat', yy)
nc_varput(outputfile, 'lon', xx)
nc_dump (outputfile);
x = nc_varget(outputfile, 'lon');
y = nc_varget(outputfile, 'lat');
z = nc_varget(outputfile, 'altitude');
