mlosch/interp_llc/shelfice_hydrography.m

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
1	function [t,s] = shelfice_hydrography(ix,iy,is,t0,s0);
2	[nx ny nz nt] = size(t0);
3	niter = 25;
4	for kt = 1:nt
5	disp(sprintf('horizontal and vertical extrapolation for time slab %u',kt))
6	ttm1 = NaN*ones(nx,ny);
7	ssm1 = NaN*ones(nx,ny);
8	for k=1:nz
9	t(:,:,k,kt) = xyexpand(t0(:,:,k,kt),25);
10	s(:,:,k,kt) = xyexpand(s0(:,:,k,kt),25);
11	% vertical extrapolation
12	tt = t(:,:,k,kt);
13	ss = s(:,:,k,kt);
14	km1 = max(k-1,1);
15	for ks = 1:length(is)
16	% nearest extrapolation from above
17	if isnan(tt(is(ks)))
18	tt(is(ks)) = ttm1(is(ks));
19	end
20	if isnan(ss(is(ks)))
21	ss(is(ks)) = ssm1(is(ks));
22	end
23	end
24	ttm1 = tt;
25	ssm1 = ss;
26	tt0 = t(:,:,k,kt);
27	ss0 = s(:,:,k,kt);
28	% apply mask
29	tt0(is) = tt(is);
30	ss0(is) = ss(is);
31	t(:,:,k,kt) = tt0;
32	s(:,:,k,kt) = ss0;
33	end
34	end
35
36	return