quicwind/skew.m

   1 function g=skew(f,h,s)
   2 % d=skew(f,h,s)
   3 % transform vector field on staggered skewed mesh so that flux balance
   4 % in the cells is just difference of g on the sides
   5 % arguments:
   6 %    f{1}, f{2}, f{3} vector components on staggered mesh (midpoints of sides)
   7 %    h(1:3) mesh step
   8 %    s{1},s{2}  terrain gradient at mesh cell midpoints
   9 % output:
  10 %    u = df{1}/x1 + df{2}/dx2 + df{3}/dx3
  11 g{1}=f{1}*(h(2)*h(3));
  12 g{2}=f{2}*(h(1)*h(3));
  13 g{3}=f{3}*(h(1)*h(2));
  14 if exist('s','var'),
  15     % average u and v to midpoint of lower side of the cells
  16     % except at
  17     uavg = (u(1:end-1,:,1:end-2)+u(2:end,:,1:end-2)
  18             +u(1:end-1,:,1:end-2)+u(2:end,:,1:end-2))/4;
  19     dvdz = (v(:,1:end-1,3:end)+v(:,2:end,3:end)...
  20           -v(:,1:end-1,1:end-2)-v(:,2:end,1:end-2))/(4*h(3));
  21     sx = s{1};
  22     sy = s{2};
  23     for k=2:size(d,3)-1
  24         d(:,:,k) = d(:,:,k) + dudz .* sx + dvdz .* sy;
  25     end
  26 end