function scoord_roms(theta_s, theta_b, Tcline, N, kgrid, column, index, varargin);
%
% originally scoord2.m but ...
% I changed a few of the commands below to better handle netcdf file opening, etc.
% Also changed a few plot commands and added a small plot of entire domain to better
% see where your cross section is being plotted.
% Some command line options need to be specified with their values - see example below.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%----jcwarner June 2002----%%%%%%%
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 1998 Rutgers University.                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
% function scoord_roms(theta_s, theta_b, Tcline,N, kgrid,column, index, varargin)%
% variable args in = ncfile, h                                              %
%                                                                           %
% This routine computes the depths of RHO- or W-points for a grid section   %
% along columns (ETA-axis) or rows (XI-axis).                               %
%                                                                           %
% Required Input:                                                           %
%                                                                           %
%    theta_s   S-coordinate surface control parameter (scalar):             %
%                [0 < theta_s < 20].                                        %
%    theta_b   S-coordinate bottom control parameter (scalar):              %
%                [0 < theta_b < 1].                                         %
%    Tcline    Width (m) of surface or bottom boundary layer in which       %
%              higher vertical resolution is required during streching      %
%              (scalar).                                                    %
%    N         Number of vertical levels (scalar).                          %
%    kgrid     Depth grid type logical switch:                              %
%                kgrid = 0   ->  depths of RHO-points.                      %
%                kgrid = 1   ->  depths of W-points.                        %
%    column    Grid direction logical switch:                               %
%                column = 1  ->  column section.                            %
%                column = 0  ->  row section.                               %
%    index     Column or row to compute (scalar):                           %
%                if column = 1, then   1 <= index <= Lp                     %
%                if column = 0, then   1 <= index <= Mp                     %
%                                                                           %
%    Variable Arguments In (varargin):                                      %
%                                                                           %
%    ncfile    Grid NetCDF filename (string).                               %
%    h         Array of depth values.                                       %
%                                                                           %
%   examples:                                                               %
%   scoord_roms(1,1,10,20,1,0,3,'ncfile','condie_his.nc');
%   this example would read the ncfile to get the variable "h" 
%
%   scoord_roms(1,1,10,20,1,0,3,'h',depth);
%   this example uses array "depth" in local workspace
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%
%  Scan varargin list and get 'h' from local workspace or from netcdf file 
%
x_label='x';
y_label='y';

if strcmp(varargin{1},'h')
  h=varargin{2};
  x_rho=[1:size(h,2)];
  x_rho=repmat(x_rho,size(h,2),1);
  y_rho=[1:size(h,1)]';
  y_rho=repmat(y_rho,1,size(h,1));
end
if(exist('h')~=1)
      if strcmp(varargin{1},'ncfile')
        ncfile = [varargin{2}];
      end
      nc = netcdf(ncfile,'nowrite');
      h = nc{'h'}(:);
      if isempty(nc{'x_rho'})
        x_rho = nc{'lon_rho'}(:,:);
        y_rho = nc{'lat_rho'}(:,:);
        if(exist('x_label')~=1),x_label='Longitude (deg)'; end
        if(exist('y_label')~=1),y_label='Latitude (deg)'; end
      else
        x_rho = nc{'x_rho'}(:,:)/1000;
        y_rho = nc{'y_rho'}(:,:)/1000;
        if(exist('x_label')~=1),x_label='Distance (km)'; end
        if(exist('y_label')~=1),y_label='Distance (km)'; end
      end
end
if (exist('x_rho')~=1)
  x_rho=size(h,2);
  y_rho=size(h,1);
  x_label='distance';
  y_label='distance';
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%----------------------------------------------------------------------------
%  Read in bathymetry and their positions.
%----------------------------------------------------------------------------


%----------------------------------------------------------------------------
%  Set several parameters.
%----------------------------------------------------------------------------

c1=1.0;
c2=2.0;
p5=0.5;
Np=N+1;
ds=1.0/N;
hmin=min(min(h));
hmax=max(max(h));
hc=min(hmin,Tcline);
[Mp Lp]=size(h);

%----------------------------------------------------------------------------
% Test input to see if it's in an acceptable form.
%----------------------------------------------------------------------------

if (column),

  if (index < 1 | index > Lp),
    disp(' ');
    disp([setstr(7),'*** Error:  SCOORD - illegal column index.',setstr(7)]);
    disp([setstr(7),'                     valid range:  1 <= index <= ',...
         num2str(Lp),setstr(7)]);
    disp(' ');
    return
  else,
    disp(' ');
    disp([' SCOORD - computing grid section depths along column: ',...
         num2str(index)]);
    disp(' ');
  end,

else,

  if (index < 1 | index > Mp),
    disp(' ');
    disp([setstr(7),'*** Error:  SCOORD - illegal row index.',setstr(7)]);
    disp([setstr(7),'                     valid range:  1 <= index <= ',...
         num2str(Mp),setstr(7)]);
    disp(' ');
    return
  else,
    disp(' ');
    disp([' SCOORD - computing grid section depths along row: ',...
         num2str(index)]);
    disp(' ');
  end,

end,

% Report S-coordinate parameters.

report = 1;
if report
disp(['          theta_s = ',num2str(theta_s)]);
disp(['          theta_b = ',num2str(theta_b)]);
disp(['          Tcline  = ',num2str(Tcline)]);
disp(['          hc      = ',num2str(hc)]);
disp(['          hmin    = ',num2str(hmin)]);
disp(['          hmax    = ',num2str(hmax)]);
disp(' ');
end

%----------------------------------------------------------------------------
% Define S-Curves at vertical RHO- or W-points (-1 < sc < 0).
%----------------------------------------------------------------------------

if (kgrid == 1),
  Nlev=Np;
  lev=0:N;
  sc=-c1+lev.*ds;
else
  Nlev=N;
  lev=1:N;
  sc=-c1+(lev-p5).*ds;
end,
Ptheta=sinh(theta_s.*sc)./sinh(theta_s);
Rtheta=tanh(theta_s.*(sc+p5))./(c2*tanh(p5*theta_s))-p5;
Cs=(c1-theta_b).*Ptheta+theta_b.*Rtheta;
Cd_r=(c1-theta_b).*cosh(theta_s.*sc)./sinh(theta_s)+ ...
     theta_b./tanh(p5*theta_s)./(c2.*(cosh(theta_s.*(sc+p5))).^2);
Cd_r=Cd_r.*theta_s;

%============================================================================
% Compute depths at requested grid section.  Assume zero free-surface.
%============================================================================

zeta=zeros(size(h));

%----------------------------------------------------------------------------
% Column section: section along ETA-axis.
%----------------------------------------------------------------------------

if (column),

  z=zeros(Mp,Nlev);
  for k=1:Nlev,
    z(:,k)=zeta(:,index).*(c1+sc(k))+hc.*sc(k)+ ...
             (h(:,index)-hc).*Cs(k);
  end,

%----------------------------------------------------------------------------
% Row section: section along XI-axis.
%----------------------------------------------------------------------------

else,

  z=zeros(Lp,Nlev);
  for k=1:Nlev,
    z(:,k)=zeta(index,:)'.*(c1+sc(k))+hc.*sc(k)+ ...
             (h(index,:)'-hc).*Cs(k);
  end,

end,

%============================================================================
% Plot grid section.
%============================================================================
figure;
%
% Plot second axes of entire grid with line at cross-section location
%

ax1=axes('position',[0.1 0.1 0.35 0.8]);
if (column),
    
    eta=y_rho(:,index)';
    eta2=[eta(1) eta eta(Mp)];
    hs=-h(:,index)';
    zmin=min(hs);
    hs=[zmin hs zmin];
    hold off;

    fill(eta2',hs,[0.6 0.7 0.6]);
    hold on;
    plot(eta',z,'k');

    %plot vertical gridlines
    for i=1:length(eta')
        plot(repmat(eta(i)',1,length(z(i,:))),z(i,:),'k-');
    end

%    set(gca,'xlim',[-Inf Inf],'ylim',[zmin 0]);
    if (kgrid == 0),
      title(['Grid (RHO-points) Section at  XI = ',num2str(index)]);
    else
      title(['Grid (W-points) Section at  XI = ',num2str(index)]);
    end,
    xlabel(x_label);
    ylabel('depth  (m)');
%    set(gcf,'Units','Normalized',...
%            'Position',[0.2 0.1 0.6 0.8],...
%            'PaperUnits','Normalized',...
%            'PaperPosition',[0.2 0.1 0.6 0.8]);
 
  else,

    xi=x_rho(index,:);
    xi2=[xi(1) xi xi(Lp)];
    hs=-h(index,:);
    zmin=min(hs);
    hs=[zmin hs zmin];
    hold off;

    fill(xi2,hs,[0.6 0.7 0.6]);

    hold on;
    plot(xi,z,'k');
    %plot vertical gridlines
    for i=1:length(xi)
        plot(repmat(xi(i),1,length(z(i,:))),z(i,:),'k-');
    end
    set(gca,'xlim',[-Inf Inf],'ylim',[zmin 0]);
    if (kgrid == 0),
      title(['Grid (RHO-points) Section at  ETA = ',num2str(index)]);
    else
      title(['Grid (W-points) Section at  ETA = ',num2str(index)]);
    end
    xlabel(x_label);
    ylabel('depth  (m)');
%    set(gcf,'Units','Normalized',...
%            'Position',[0.2 0.1 0.6 0.8],...
%            'PaperUnits','Normalized',...
%            'PaperPosition',[0.2 0.1 0.6 0.8]);

end

%pcolor of depth

ax2=axes('position',[0.55 0.1 0.35 0.8]);
%imagesc(h)

if (1)
%mask_rho = nc{'mask_rho'}(:,:);
%if isempty(mask_rho)
%  mask_rho=ones(size(h));
%end
%zz=find(mask_rho==0);
%mask_rho(zz)=nan;
%pcolor(x_rho,y_rho,h.*mask_rho)
pcolor(h)
%shading flat
zz=colorbar;
set(get(zz,'ylabel'),'string','Bathymetry');
hold on
%if (column)
%  plot([index index], [0 size(h,1)],'linewidth',[2],'color','k')
%else
%  plot([0 size(h,2)],[index index],'linewidth',[2],'color','k')
%end
if (column)
%  plot([x_rho(1,index) x_rho(1,index) ], [y_rho(1,index) y_rho(end,index)],'linewidth',[2],'color','k')
%  plot(x_rho(:,index) , y_rho(:,index) ,'linewidth',[2],'color','k')
else
%  plot([x_rho(index,1) x_rho(index,end) ], [y_rho(index,1) y_rho(index,1)],'linewidth',[2],'color','k')
%  plot(x_rho(index,:) ,y_rho(index,:) ,'linewidth',[2],'color','k')
end

xlabel(x_label);
ylabel(x_label);

end