eddy_flux/c22/calc_dbdy.m

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


1	edhill	1.1	function [dbdy] = calc_dbdy(b, dut,dvt, lluy,llvy)
2
3			sz = size(b);
4			ne = sz(1);
5
6			dbdy = zeros(size(b));
7
8			for k = 1:6
9			for i = 1:ne
10			for j = 1:ne
11
12			num = 0;
13			sum = 0.0;
14			if (i > 1)
15			num = num + 1;
16			sum = sum + ...
17			lluy(i,j,k)*(b(i,j,k) - b(i-1,j,k))/dut(i,j,k);
18			end
19			if (i < ne)
20			num = num + 1;
21			sum = sum + ...
22			lluy(i,j,k)*(b(i+1,j,k) - b(i,j,k))/dut(i,j,k);
23			end
24			if (j > 1)
25			num = num + 1;
26			sum = sum + ...
27			llvy(i,j,k)*(b(i,j,k) - b(i,j-1,k))/dvt(i,j,k);
28			end
29			if (j < 1)
30			num = num + 1;
31			sum = sum + ...
32			llvy(i,j,k)*(b(i,j+1,k) - b(i,j,k))/dvt(i,j,k);
33			end
34
35			if num > 0
36			dbdy(i,j,k) = sum / num;
37			end
38
39			end
40			end
41
42			end
43
44
45