quicwind/grad3z.m

   1 function g=grad3z(f,h,b)
   2 % u=grad3z(f,h)
   3 % compute gradient in 3d assuming zero boundary conditions
   4 % u=grad3z(f,h,b)
   5 % in 3rd coordinate use reflection at bottom instead of zero
   6 % input:
   7 %    f       3d array
   8 %    h(1:3)  stepsize
   9 %    b       use zero boundary condition on output at the bottom
  10
  11 if ~exist('b','var')
  12     b = false;
  13 end
  14
  15 % wrap in zeros
  16 n = size(f);
  17 % fz = zeros(n+2);
  18 fz(2:n(1)+1,2:n(2)+1,2:n(3)+1)=f;
  19 if b,
  20     fz(:,:,1) = fz(:,:,2);
  21 end
  22
  23 % derivatives
  24 g{1}=(fz(2:end,2:end-1,2:end-1)-fz(1:end-1,2:end-1,2:end-1))/h(1);
  25 g{2}=(fz(2:end-1,2:end,2:end-1)-fz(2:end-1,1:end-1,2:end-1))/h(2);
  26 g{3}=(fz(2:end-1,2:end-1,2:end)-fz(2:end-1,2:end-1,1:end-1))/h(3);
  27 end
  28