function [z]=depths(fname, gname, igrid, idims, tindex);
%
% DEPTHS:  Compute the ROMS depths associated with a 3D variable
%          Modified 20180806 for the presence of ice shelves.
% [z]=depths(fname, gname, igrid, idims, tindex)
%
% This function computes the depths at the requested staggered C-grid.
% If the time record (tindex) is not provided, a zero free-surface is
% assumed and the unperturbed depths are returned.
%
% On Input:
%
%    fname       NetCDF data file name (character string)
%    gname       NetCDF grid file name (character string)
%    igrid       Staggered grid C-type (integer):
%                  igrid=1  => density points
%                  igrid=2  => streamfunction points
%                  igrid=3  => u-velocity points
%                  igrid=4  => v-velocity points
%                  igrid=5  => w-velocity points
%    idims       Depths dimension order switch (integer):
%                  idims=0  => (i,j,k)  column-major order (Fortran)
%                  idims=1  => (j,i,k)  row-major order (C-language)
%    tindex      Time index (integer)
%
% On Output:
%
%    z           Depths (3D array; meters, negative)
%

% svn $Id: depths.m 754 2015-01-07 23:23:40Z arango $
%=========================================================================%
%  Copyright (c) 2002-2015 The ROMS/TOMS Group                            %
%    Licensed under a MIT/X style license                                 %
%    See License_ROMS.txt                           Hernan G. Arango      %
%=========================================================================%

% Check arguments and set default values.

if (nargin < 3)
  igrid=1;
end,

if (nargin < 4),
  idims=0;
end,

if (nargin < 5),
  tindex=0;
end,

% Initialize vertical tranformation and stretching function to
% original formulation values.

%too dangerous; psl20201222 Vtransform=1;
%too dangerous; psl20201222 Vstretching=1;
got.hc=false;

%--------------------------------------------------------------------------
% Read in S-coordinate parameters.
%--------------------------------------------------------------------------

S=ncinfo( fname );
% S=nc_vnames(fname);
nvars=length(S.Variables);

for n=1:nvars
  name=char(S.Variables(n).Name);
  switch name
    case 'Vtransform'
      Vtransform=ncread(fname,'Vtransform');
    case 'Vstretching'
      Vstretching=ncread(fname,'Vstretching');
    case 'hc'
      got.hc=true;
    case 'sc_r'
      name_r='sc_r';
    case 's_rho'
      name_r='s_rho';
    case 'sc_w'
      name_w='sc_w';
    case 's_w'
      name_w='s_w';
    case 'zice'
      zice = ncread( fname, 'zice' );
  end
end

sc_r=ncread(fname,name_r);
Cs_r=ncread(fname,'Cs_r');

sc_w=ncread(fname,name_w);
Cs_w=ncread(fname,'Cs_w');

N=length(sc_r);
Np=N+1;

if (length(sc_w) == N),
  sc_w=[-1 sc_w'];
  Cs_w=[-1 Cs_w'];
end

%------------------------------------------------------------------------
% Get bottom topography.
%------------------------------------------------------------------------

h=ncread(gname,'h');
[Lp Mp]=size(h);
L=Lp-1;
M=Mp-1;

  if ( exist( 'zice' ) ~= 1 )
    zice = zeros( Lp, Mp );
  end

switch ( igrid ),
  case 1
    if (idims), h=h'; zice=zice'; end,
%   if (idims), h=h'; end,
  case 2
    hp=0.25.*(h(1:L,1:M)+h(2:Lp,1:M)+h(1:L,2:Mp)+h(2:Lp,2:Mp));
    zicep=0.25.*(zice(1:L,1:M)+zice(2:Lp,1:M)+zice(1:L,2:Mp)+zice(2:Lp,2:Mp));
    if (idims), hp=hp'; zicep=zicep'; end,
%   if (idims), hp=hp'; end,
  case 3
    hu=0.5.*(h(1:L,1:Mp)+h(2:Lp,1:Mp));
    ziceu=0.5.*(zice(1:L,1:Mp)+zice(2:Lp,1:Mp));
    if (idims), hu=hu'; ziceu=ziceu'; end,
%   if (idims), hu=hu'; end,
  case 4
    hv=0.5.*(h(1:Lp,1:M)+h(1:Lp,2:Mp));
    zicev=0.5.*(zice(1:Lp,1:M)+zice(1:Lp,2:Mp));
    if (idims), hv=hv'; zicev=zicev'; end,
%   if (idims), hv=hv'; end,
  case 5
    if (idims), h=h'; zice=zice'; end,
%   if (idims), h=h'; end,
end

% Set critical depth parameter.

if (got.hc),
  hc=ncread(fname,'hc');
else
  hc=min(min(h));
end

%------------------------------------------------------------------------
% Get free-surface
%------------------------------------------------------------------------

if (tindex == 0),
  zeta=zeros([Lp Mp]);
else
% psl:
  if ( any( size( h ) == 1 ) ) % If a station    file...
    zeta=ncread(fname,'zeta', [1,    tindex], [inf,      1], [1,    1] );
  else                         % If a his or avg file...
    zeta=ncread(fname,'zeta', [1, 1, tindex], [inf, inf, 1], [1, 1, 1] );
  end
  zeta( find( isnan( zeta ) ) ) = 0.;
% zeta=ncread(fname,'zeta',tindex);
end

switch ( igrid ),
  case 1
    if (idims), zeta=zeta'; end,
  case 2
    zetap=0.25.*(zeta(1:L,1:M )+zeta(2:Lp,1:M )+                      ...
                 zeta(1:L,2:Mp)+zeta(2:Lp,2:Mp));
    if (idims), zetap=zetap'; end,
  case 3
    zetau=0.5.*(zeta(1:L,1:Mp)+zeta(2:Lp,1:Mp));
    if (idims), zetau=zetau'; end,
  case 4
    zetav=0.5.*(zeta(1:Lp,1:M)+zeta(1:Lp,2:Mp));
    if (idims), zetav=zetav'; end,
  case 5
    if (idims), zeta=zeta'; end,
end

%------------------------------------------------------------------------
% Compute depths.
%------------------------------------------------------------------------

if (Vtransform == 1),
  switch ( igrid ),
    case 1
      for k=1:N,
	z0=(sc_r(k)-Cs_r(k))*hc + Cs_r(k).*(h+zice);
        z(:,:,k)=z0 + zeta.*(1.0 + z0./(h+zice)) + zice;
%       z0=(sc_r(k)-Cs_r(k))*hc + Cs_r(k).*h;
%       z(:,:,k)=z0 + zeta.*(1.0 + z0./h);
      end
    case 2
      for k=1:N,
        z0=(sc_r(k)-Cs_r(k))*hc + Cs_r(k).*(hp+zicep);
        z(:,:,k)=z0 + zetap.*(1.0 + z0./(hp+zicep)) + zicep;
%       z0=(sc_r(k)-Cs_r(k))*hc + Cs_r(k).*hp;
%       z(:,:,k)=z0 + zetap.*(1.0 + z0./hp);
      end
    case 3
      for k=1:N,
        z0=(sc_r(k)-Cs_r(k))*hc + Cs_r(k).*(hu+ziceu);
        z(:,:,k)=z0 + zetau.*(1.0 + z0./(hu+ziceu)) + ziceu;
%       z0=(sc_r(k)-Cs_r(k))*hc + Cs_r(k).*hu;
%       z(:,:,k)=z0 + zetau.*(1.0 + z0./hu);
      end
    case 4
      for k=1:N,
        z0=(sc_r(k)-Cs_r(k))*hc + Cs_r(k).*(hv+zicev);
        z(:,:,k)=z0 + zetav.*(1.0 + z0./(hv+zicev)) + zicev;
%       z0=(sc_r(k)-Cs_r(k))*hc + Cs_r(k).*hv;
%       z(:,:,k)=z0 + zetav.*(1.0 + z0./hv);
      end
    case 5
      z(:,:,1)=-h;
      for k=2:Np,
        z0=(sc_w(k)-Cs_w(k))*hc + Cs_w(k).*(h+zice);
        z(:,:,k)=z0 + zeta.*(1.0 + z0./(h+zice)) + zice;
%       z0=(sc_w(k)-Cs_w(k))*hc + Cs_w(k).*h;
%       z(:,:,k)=z0 + zeta.*(1.0 + z0./h);
      end
  end
elseif (Vtransform == 2),
  switch ( igrid ),
    case 1
      for k=1:N,
        z0=(hc.*sc_r(k)+Cs_r(k).*(h+zice))./(hc+h+zice);
        z(:,:,k)=zeta+(zeta+h+zice).*z0;
%       z0=(hc.*sc_r(k)+Cs_r(k).*h)./(hc+h);
%       z(:,:,k)=zeta+(zeta+h).*z0;
      end,
    case 2
      for k=1:N,
        z0=(hc.*sc_r(k)+Cs_r(k).*(hp+zicep))./(hc+hp+zicep);
        z(:,:,k)=zetap+(zetap+hp+zicep).*z0;
%       z0=(hc.*sc_r(k)+Cs_r(k).*hp)./(hc+hp);
%       z(:,:,k)=zetap+(zetap+hp).*z0;
      end,
    case 3
      for k=1:N,
        z0=(hc.*sc_r(k)+Cs_r(k).*(hu+ziceu))./(hc+hu+ziceu);
        z(:,:,k)=zetau+(zetau+hu+ziceu).*z0;
%       z0=(hc.*sc_r(k)+Cs_r(k).*hu)./(hc+hu);
%       z(:,:,k)=zetau+(zetau+hu).*z0;
      end,
    case 4
      for k=1:N,
        z0=(hc.*sc_r(k)+Cs_r(k).*(hv+zicev))./(hc+hv+zicev);
        z(:,:,k)=zetav+(zetav+hv+zicev).*z0;
%       z0=(hc.*sc_r(k)+Cs_r(k).*hv)./(hc+hv);
%       z(:,:,k)=zetav+(zetav+hv).*z0;
      end,
    case 5
      for k=1:Np,
        z0=(hc.*sc_w(k)+Cs_w(k).*(h+zice))./(hc+h+zice);
        z(:,:,k)=zeta+(zeta+h+zice).*z0;
%       z0=(hc.*sc_w(k)+Cs_w(k).*h)./(hc+h);
%       z(:,:,k)=zeta+(zeta+h).*z0;
      end
  end
end

return
