% Create a simple rectangular grid for KS fjord

% Set-up Lambert projection 
mstruct = defaultm('lambertstd');
mstruct.mapparallels = 71.42;
mstruct.origin = [71.42 -52.6]; % center of grid
mstruct = defaultm(mstruct);

% Grid dimensions
Lp = 3*2^7+2;
Mp = 3*2^5+2;

% Build rectangular grid from central point
[I,J] = meshgrid(1:Lp,1:Mp);
I = I-Lp/2;
J = J-Mp/2;

% Trial and error: Guess at suitable domain size in projected coordinates
% to get dx,dy you want. An alternative would be to experiment with 
% projinv calls to relate actual distance in latitude to dy units.
XL = 0.014;  % domain size in Lambert - not sure how this relates to geo
EL = 0.0035;

% grid spacing in Lambert coordinates
dx = XL/Lp;
dy = EL/Mp;

% convert i,j to dx,dy units 
X = dx*I;
Y = dy*J;

% option to roate the grid from east-north
theta = 5*pi/180;
if theta ~= 0
  % rotate the grid
  R = [cos(theta) -sin(theta); sin(theta) cos(theta)];
  XY = [X(:) Y(:)]';
  XY2 = R*XY;
  X = reshape(XY2(1,:),size(X));
  Y = reshape(XY2(2,:),size(Y));
end

% Convert projected coordinate back to lon/lat
% [rlat,rlon] = minvtran(mstruct,X,Y); 
[rlat,rlon] = projinv(mstruct,X,Y); % gives slightly different results 

% Coastline for plotting - this loads GSHHS files that can be downoaded
% from various sources. Here, I have kept them with the myroms.org 
% distribution of m_map (for historical reasons). 
clear S
if ~exist('S','var')
  % find where I keep the GSHHS coastline files with m_map toolbox
  dd = dir(which('m_grid'));
  % resolution options coarse, low, intermediate, high, fine
  % use                c       l    i             h     f 
  file = fullfile(dd.folder,'private','gshhs_i.b');
  %                                         ^^ resolution
  indexfile = gshhs(file, 'createindex');
  S = gshhs(file,range(rlat),range(rlon));
end

% Plot
clf
plot(rlon,rlat,'r.')
a = axis;
hold on
for k=1:size(S)
  han = plot(S(k).Lon,S(k).Lat,'k-');
  han.LineWidth = 2;
end
hold off
axis(a)

% write the roms grid file
% set_grid is part or myroms.org matlab tools
set_grid(rlon, rlat, 'ks_test.nc','linear','i')

%%
pm = ncread('ks_test.nc','pm');
pn = ncread('ks_test.nc','pn');
disp(['range of 1/pm (i.e. delta-x) is ' num2str(min(1./(pm(:)))) ...
  ' to ' num2str(max(1./(pm(:)))) ' meters'])
disp(['range of 1/pn (i.e. delta-y) is ' num2str(min(1./(pn(:)))) ...
  ' to ' num2str(max(1./(pn(:)))) ' meters'])