function [t,s] = shelfice_hydrography(ix,iy,is,t0,s0);
  [nx ny nz nt] = size(t0);
  niter = 25;
  for kt = 1:nt
    disp(sprintf('horizontal and vertical extrapolation for time slab %u',kt))
    ttm1 = NaN*ones(nx,ny);
    ssm1 = NaN*ones(nx,ny);
    for k=1:nz
      t(:,:,k,kt) = xyexpand(t0(:,:,k,kt),25);
      s(:,:,k,kt) = xyexpand(s0(:,:,k,kt),25);
      % vertical extrapolation
      tt = t(:,:,k,kt);
      ss = s(:,:,k,kt);
      km1 = max(k-1,1);
      for ks = 1:length(is)
	% nearest extrapolation from above
	if isnan(tt(is(ks)))
	  tt(is(ks)) = ttm1(is(ks));
	end
	if isnan(ss(is(ks)))
	  ss(is(ks)) = ssm1(is(ks));
	end
      end
      ttm1 = tt;
      ssm1 = ss;
      tt0 = t(:,:,k,kt);
      ss0 = s(:,:,k,kt);
      % apply mask
      tt0(is) = tt(is);
      ss0(is) = ss(is);
      t(:,:,k,kt) = tt0;
      s(:,:,k,kt) = ss0;
    end
  end

  return