function [dbdy] = calc_dbdy(b, dut,dvt, lluy,llvy)

sz = size(b);
ne = sz(1);

dbdy = zeros(size(b));

for k = 1:6
  for i = 1:ne
    for j = 1:ne

      num = 0;
      sum = 0.0;
      if (i > 1)
        num = num + 1;
        sum = sum + ...
              lluy(i,j,k)*(b(i,j,k) - b(i-1,j,k))/dut(i,j,k);
      end
      if (i < ne)
        num = num + 1;
        sum = sum + ...
              lluy(i,j,k)*(b(i+1,j,k) - b(i,j,k))/dut(i,j,k);
      end
      if (j > 1)
        num = num + 1;
        sum = sum + ...
              llvy(i,j,k)*(b(i,j,k) - b(i,j-1,k))/dvt(i,j,k);
      end
      if (j < 1)
        num = num + 1;
        sum = sum + ...
              llvy(i,j,k)*(b(i,j+1,k) - b(i,j,k))/dvt(i,j,k);
      end

      if num > 0
        dbdy(i,j,k) = sum / num;
      end
      
    end
  end

end