function [dz,z]=rocky_grd(h,alpha,dz1,n);
% function [dz,z]=test_grd(h,alpha,dz1,n);
%     fuunction for making log spaced boundary layer that merges smoothly 
%                               into uniform vertical spacing
%                  h is depth
%                  alpha is factor of increease in dz (greater than 1,
%                                                              less than 2)
%                  dz1 is the spacing of the bottom grid in meters
%                  n is the number of grid cells
%                  dz is output delta z (length n)
%                  z is output height above bottom (length(n)
%                  note that these variables will have to be converted to
%                  conform to roms           
%   note that the bottom grid is not exactly as specified but it is close.

if nargin<4;   %  these are reasonable default parameters
  alpha=1.5;   %  factor of increase in dz for each vertical level in boundary layer
  dz1=.1;      %  dz of bottom grid cell (10 cm) 
  n=20;        %  number of grid cells
end;

if nargin==0;
  h=10;        %  depth  (meters) 
end;

[Im,Jm]=size(h);

k=1:n;         %  index array of sigma levels

dzbl0=dz1*alpha.^(k-1);  % initial boundary layer dz array 

for j=1:Jm,
  for i=1:Im,
    nbl=log( h(i,j)/(n-2)/dz1 )/log(alpha);   % sigma level of transition

    if (nbl < .1),    %  this takes care of cases in which the depth
      nbl=.1;         %  is too shallow for there to be a boundary layer grid
    end,   

    fac=(tanh(5.*(1-k./nbl))+1)./2;  %  hyperbolic tangent weighting
                                     %  function for boundary layer

    dzbl=dzbl0.*fac ;                %  boundary layer grid 
    zbl=sum(dzbl);                   %  height contribution of boundary layer
    dz_outer=(1-fac)*(h-zbl)/sum(1-fac);  %  outer grid

    dz=dzbl+dz_outer;                %  combined grid of dz
    z=cumsum(dz);                   %  heights
  end,
end,

plot(dz,-z,'.');
