gmaze_pv/subfct/subfct_getdV.m

% DV = subfct_getdV(DEPTH,LATITUDE,LONGITUDE)
% Compute 3D volume elements matrix from geographical
% axis Z(<0,downward), Y and X

function DV = subfct_getdV(Z,Y,X)

nz = length(Z);
ny = length(Y);
nx = length(X);

DV = zeros(nz,ny,nx);

% Vertical elements:
for iz = 1 : nz % Toward the deep ocean (because DPT<0)
        % Vertical grid length centered at Z(iy)
        if iz == 1
          dz = abs(Z(1)) + abs(sum(diff(Z(iz:iz+1))/2));
        elseif iz == nz % We don't know the real ocean depth
          dz = abs(sum(diff(Z(iz-1:iz))/2));
        else
          dz = abs(sum(diff(Z(iz-1:iz+1))/2));
        end
        DZ(iz) = dz;
end

% Surface and Volume elements:
for ix = 1 : nx
  for iy = 1 : ny
      % Zonal grid length centered in X(ix),Y(iY)
      if ix == 1
         dx = abs(m_lldist([X(ix) X(ix+1)],[1 1]*Y(iy)))/2;
      elseif ix == nx 
         dx = abs(m_lldist([X(ix-1) X(ix)],[1 1]*Y(iy)))/2;
      else
         dx = abs(m_lldist([X(ix-1) X(ix)],[1 1]*Y(iy)))/2+abs(m_lldist([X(ix) X(ix+1)],[1 1]*Y(iy)))/2;
      end       
 
      % Meridional grid length centered in X(ix),Y(iY)
      if iy == 1
        dy = abs(m_lldist([1 1]*X(ix),[Y(iy) Y(iy+1)]))/2;
      elseif iy == ny
        dy = abs(m_lldist([1 1]*X(ix),[Y(iy-1) Y(iy)]))/2;
      else      
        dy = abs(m_lldist([1 1]*X(ix),[Y(iy-1) Y(iy)]))/2+abs(m_lldist([1 1]*X(ix),[Y(iy) Y(iy+1)]))/2;
      end

      % Surface element:
      DA = dx*dy.*ones(1,nz);
      
      % Volume element:
      DV(:,iy,ix) = DZ.*DA;
  end %for iy
end %for ix

1	% DV = subfct_getdV(DEPTH,LATITUDE,LONGITUDE)
2	% Compute 3D volume elements matrix from geographical
3	% axis Z(<0,downward), Y and X
4
5	function DV = subfct_getdV(Z,Y,X)
6
7	nz = length(Z);
8	ny = length(Y);
9	nx = length(X);
10
11	DV = zeros(nz,ny,nx);
12
13	% Vertical elements:
14	for iz = 1 : nz % Toward the deep ocean (because DPT<0)
15	% Vertical grid length centered at Z(iy)
16	if iz == 1
17	dz = abs(Z(1)) + abs(sum(diff(Z(iz:iz+1))/2));
18	elseif iz == nz % We don't know the real ocean depth
19	dz = abs(sum(diff(Z(iz-1:iz))/2));
20	else
21	dz = abs(sum(diff(Z(iz-1:iz+1))/2));
22	end
23	DZ(iz) = dz;
24	end
25
26	% Surface and Volume elements:
27	for ix = 1 : nx
28	for iy = 1 : ny
29	% Zonal grid length centered in X(ix),Y(iY)
30	if ix == 1
31	dx = abs(m_lldist([X(ix) X(ix+1)],[1 1]*Y(iy)))/2;
32	elseif ix == nx
33	dx = abs(m_lldist([X(ix-1) X(ix)],[1 1]*Y(iy)))/2;
34	else
35	dx = abs(m_lldist([X(ix-1) X(ix)],[1 1]Y(iy)))/2+abs(m_lldist([X(ix) X(ix+1)],[1 1]Y(iy)))/2;
36	end
37
38	% Meridional grid length centered in X(ix),Y(iY)
39	if iy == 1
40	dy = abs(m_lldist([1 1]*X(ix),[Y(iy) Y(iy+1)]))/2;
41	elseif iy == ny
42	dy = abs(m_lldist([1 1]*X(ix),[Y(iy-1) Y(iy)]))/2;
43	else
44	dy = abs(m_lldist([1 1]X(ix),[Y(iy-1) Y(iy)]))/2+abs(m_lldist([1 1]X(ix),[Y(iy) Y(iy+1)]))/2;
45	end
46
47	% Surface element:
48	DA = dxdy.ones(1,nz);
49
50	% Volume element:
51	DV(:,iy,ix) = DZ.*DA;
52	end %for iy
53	end %for ix
54