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	gmaze	1.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