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	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