%Make UVic ESCM grid for matrix extraction
% K Kvale November 2014
clear all;

sphericalGrid=1;
gridType='spherical';
deltaT=108000; 

ncid=netcdf.open('UVIC/UVic_template/data/G_grid.nc');

varid=netcdf.inqVarID(ncid,'longitude');
x=netcdf.getVar(ncid,varid,'double'); % longitude of grid box centres (degrees x)
x=x(2:end-1);
nx=length(x);
varid=netcdf.inqVarID(ncid,'latitude');
y=netcdf.getVar(ncid,varid,'double'); % latitude of grid box centres (degrees y)
y=y(2:end-1);
ny=length(y);
varid=netcdf.inqVarID(ncid,'depth');
z=netcdf.getVar(ncid,varid,'double'); % depth in meters
nz=length(z);

varid=netcdf.inqVarID(ncid,'G_dxt');
dphi=netcdf.getVar(ncid,varid,'double');
dphi=dphi(2:end-1);
dphi=repmat(dphi,[1 ny]);

varid=netcdf.inqVarID(ncid,'G_dyt');
dth=netcdf.getVar(ncid,varid,'double');
dth=dth(2:end-1);
dth=repmat(dth',[nx 1]);

varid=netcdf.inqVarID(ncid,'G_dzt');
dznom=netcdf.getVar(ncid,varid,'double'); % box thickness in meters

netcdf.close(ncid)

dz=repmat(dznom,[1 nx ny]);
dz=permute(dz,[2 3 1]);

% Temporary secondary data file load
ncid=netcdf.open('tavg.00001.01.11.nc');
  
% figure out dry/wet mask
varid=netcdf.inqVarID(ncid,'G_kmt');
kmt=netcdf.getVar(ncid,varid,'double'); % kmt
kmt=repmat(kmt,[1 1 nz]);
K=repmat([1:nz]',[1 nx ny]);
K=permute(K,[2 3 1]);
bathy=1*(K<=kmt);
kd=find(isnan(bathy) | bathy==0); % dry points
kw=find(~isnan(bathy) & bathy~=0); % wet points
bathy(kd)=0;
bathy(kw)=1;

dz=dz.*bathy;

% find index of deepest wet point 
ideep=kmt(:,:,1);

% cell volumes
varid=netcdf.inqVarID(ncid,'G_areaT');
da=netcdf.getVar(ncid,varid,'double'); % cell areas
da=repmat(da,[1 1 nz]).*bathy;
dv=da.*dz;

% grid width in m
varid=netcdf.inqVarID(ncid,'G_dxt');
dx=netcdf.getVar(ncid,varid,'double'); % width t grid
dx=repmat(dx,[1 ny]);
varid=netcdf.inqVarID(ncid,'G_dyt');
dy=netcdf.getVar(ncid,varid,'double'); % width t grid
dy=repmat(dy',[nx 1]);

netcdf.close(ncid)


% 
% % Assume spherical grid
% a=6370e3; % radius in UVIC
% 
% dy=netcdf.getVar(ncid,5,'double');% height of boxes in meters
% y1=abs(y); %get rid of negative latitudes
% y2=y1*pi/180; % convert from degrees to radians
% cir=2*pi*a*cos(y2); %circumference of Earth at each latitude
% dx=cir/100; %just the slice between two longitudes
% dy=repmat(dy,[1 nx]);
% %dy=permute(dy,[2 1]);
% dx=repmat(dx,[1 nx]);
% %dx=permute(dx,[2 1]);

save grid kd kw nx ny nz x y z dv ideep bathy dznom dz dx dy dth dphi ...
	 da sphericalGrid deltaT gridType
