quicwind/div3.m

   1 function d=div3(f,h,s)
   2 % d=div3(f,h)
   3 % compute divergence of vector field
   4 % arguments:
   5 %    f{1}, f{2}, f{3} vector components on staggered mesh (midpoints of sides)
   6 %    h(1:3) mesh step
   7 %    s{1},s{2}  terrain gradient at mesh cell midpoints
   8 % output:
   9 %    u = df{1}/x1 + df{2}/dx2 + df{3}/dx3
  10 if ~exist('h','var')
  11     h=[1,1,1];
  12 end
  13 u=f{1};
  14 v=f{2};
  15 w=f{3};
  16 dudx = (u(2:end,:,:)-u(1:end-1,:,:))/h(1);
  17 dvdy = (v(:,2:end,:)-v(:,1:end-1,:))/h(2);
  18 dwdz = (w(:,:,2:end)-w(:,:,1:end-1))/h(3);
  19 d = dudx + dvdy + dwdz;
  20 if exist('s','var')
  21     error('slant not supported');
  22     %
  23     dudz = (u(1:end-1,:,3:end)+u(2:end,:,3:end)...
  24           -u(1:end-1,:,1:end-2)-u(2:end,:,1:end-2))/(4*h(3));
  25     dvdz = (v(:,1:end-1,3:end)+v(:,2:end,3:end)...
  26           -v(:,1:end-1,1:end-2)-v(:,2:end,1:end-2))/(4*h(3));
  27     sx = s{1};
  28     sy = s{2};
  29     for k=2:size(d,3)-1
  30         d(:,:,k) = d(:,:,k) + dudz .* sx + dvdz .* sy;
  31     end
  32 end
  33