gmaze_pv/subfct/diagVOLU.m

% [V,Cm,E,Vt,CC] = diagVOLU(FLAG,C1,C2,CLASS,LON,LAT,DPT,DV,[Ca(Z,Y,X),Cb(Z,Y,X),...])
%
% DESCRIPTION:
% Compute the volume of water for a particular CLASS of potential
% temperature or density.
% Also compute mean values of additional 3D fields (such as Ca, Cb ...) along
% the CLASS of the analysed field.
%
% The volume is accounted as:
%   CLASS(i) <= FIELD < CLASS(i+1)
%
% INPUTS: 
% FLAG    : Can either be: 0, 1 or 2
%           0: Compute volume of potential density classes
%              from C1=THETA and C2=SALT
%           1: Compute volume of potential density classes
%              from C1=SIGMA_THETA
%           2: Compute volume of temperature classes
%              from C1=THETA
% C1,C2   : Depends on option FLAG:
%           - FLAG = 0 : 
%                        C1 : Temperature (^oC)
%                        C2 : Salinity (PSU)
%           - FLAG = 1 : 
%                        C1 : Potential density (kg/m3) 
%                        C2 : Not used
%           - FLAG = 2 : 
%                        C1 : Temperature (^oC)
%                        C2 : Not used
% ClASS   : Range to explore (eg: [20:.1:30] for potential density)
% LON,LAT,DPT : axis (DPT < 0)
% dV      : Matrix of grid volume elements (m3) centered in (lon,lat,dpt) 
% Ca,Cb,...: Any additional 3D fields (unlimited)  
%
%
% OUTPUTS:
% V       : Volume of each CLASS (m3)
% Cm      : Mean value of the classified field (allow to check errors)
% E       : Each time a grid point is counted, a 1 is added to this 3D matrix
%           Allow to check double count of a point or unexplored areas
% Vt      : Is the total volume explored (Vt)
% CC      : Contains the mean value of additional fields Ca, Cb ....
%
% NOTES:
% - Fields are on the format: C(DPT,LAT,LON)
% - The potential density is computed with the equation of state routine from
%   the MITgcm called densjmd95.m 
%   (see: http://mitgcm.org/cgi-bin/viewcvs.cgi/MITgcm_contrib/gmaze_pv/subfct/densjmd95.m)
% - if dV is filled with NaN, dV is computed by the function
%
%
% AUTHOR: 
% Guillaume Maze / MIT
% 
% HISTORY:
% - Created: 06/29/2007
%

% 

function varargout = diagVOLU(FLAG,C1,C2,CLASS,LON,LAT,DPT,DV,varargin)


% 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPROC
% Variables:
ndpt = size(C1,1);
nlat = size(C1,2);
nlon = size(C1,3);
CLASS = sort(CLASS(:));
[Z b c] = meshgrid(DPT,LON,LAT);clear b c, Z = permute(Z,[2 3 1]);

% Determine fields from which we'll take class contours:
switch FLAG
  
 case 0 % Need to compute SIGMA THETA
  THETA = C1;
  SALT  = C2;
  ST = densjmd95(SALT,THETA,0.09998*9.81*abs(Z)) - 1000; 
  CROP = ST;
  
 case 1
  ST = C1; % Potential density
  CROP = ST;
  
 case 2
  THETA = C1; % Potential temperature
  CROP = THETA;
end
  
% Volume elements:
if length(find(isnan(DV)==1)) == ndpt*nlat*nlon
  if exist('subfct_getdV','file')
    DV = subfct_getdV(DPT,LAT,LON);
  else
    DV  = local_getdV(DPT,LAT,LON);
  end
end

% Need to compute volume integral over these 3D fields
nIN = nargin-8;
if nIN >= 1
  doEXTRA = 1;
else
  doEXTRA = 0;
end

% 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% VOLUME INTEGRATION
explored = zeros(ndpt,nlat,nlon);
% Volume integral:
for iC = 1 : length(CLASS)-1
  mask   = zeros(ndpt,nlat,nlon);
  mask(find( (CLASS(iC) <= CROP) & (CROP < CLASS(iC+1)) )) = 1;
  explored = explored + mask;
  VOL(iC) = nansum(nansum(nansum(DV.*mask,1),2),3);
  
  if VOL(iC) ~= 0
     CAR(iC) = nansum(nansum(nansum(CROP.*DV.*mask,1),2),3)./VOL(iC);
     if doEXTRA
       for ii = 1 : nIN
           C = varargin{ii};
           CAREXTRA(ii,iC) = nansum(nansum(nansum(C.*DV.*mask,1),2),3)./VOL(iC);
       end %for ii      
     end %if doEXTRA
  else
     CAR(iC) = NaN;
     if doEXTRA
       for ii = 1 : nIN
           CAREXTRA(ii,iC) = NaN;
       end %for ii      
     end %if doEXTRA
  end  
end %for iC

% In order to compute the total volume of the domain:
CROP(find(isnan(CROP)==0)) = 1;  
CROP(find(isnan(CROP)==1)) = 0;
Vt = nansum(nansum(nansum(DV.*CROP,1),2),3);

% 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% OUTPUTS
switch nargout
 case 1
  varargout(1) = {VOL};
 case 2
  varargout(1) = {VOL};
  varargout(2) = {CAR};
 case 3
  varargout(1) = {VOL};
  varargout(2) = {CAR};
  varargout(3) = {explored}; 
 case 4
  varargout(1) = {VOL};
  varargout(2) = {CAR};
  varargout(3) = {explored}; 
  varargout(4) = {Vt};
 case 5
  varargout(1) = {VOL};
  varargout(2) = {CAR};
  varargout(3) = {explored}; 
  varargout(4) = {Vt};
  varargout(5) = {CAREXTRA}; 
end %switch


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes the 3D dV volume elements.
% Copy of the subfct_getDV function from gmaze_pv package
function DV = local_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	% [V,Cm,E,Vt,CC] = diagVOLU(FLAG,C1,C2,CLASS,LON,LAT,DPT,DV,[Ca(Z,Y,X),Cb(Z,Y,X),...])
2			%
3			% DESCRIPTION:
4			% Compute the volume of water for a particular CLASS of potential
5			% temperature or density.
6			% Also compute mean values of additional 3D fields (such as Ca, Cb ...) along
7			% the CLASS of the analysed field.
8			%
9			% The volume is accounted as:
10			% CLASS(i) <= FIELD < CLASS(i+1)
11			%
12			% INPUTS:
13			% FLAG : Can either be: 0, 1 or 2
14			% 0: Compute volume of potential density classes
15			% from C1=THETA and C2=SALT
16			% 1: Compute volume of potential density classes
17			% from C1=SIGMA_THETA
18			% 2: Compute volume of temperature classes
19			% from C1=THETA
20			% C1,C2 : Depends on option FLAG:
21			% - FLAG = 0 :
22			% C1 : Temperature (^oC)
23			% C2 : Salinity (PSU)
24			% - FLAG = 1 :
25			% C1 : Potential density (kg/m3)
26			% C2 : Not used
27			% - FLAG = 2 :
28			% C1 : Temperature (^oC)
29			% C2 : Not used
30			% ClASS : Range to explore (eg: [20:.1:30] for potential density)
31			% LON,LAT,DPT : axis (DPT < 0)
32			% dV : Matrix of grid volume elements (m3) centered in (lon,lat,dpt)
33			% Ca,Cb,...: Any additional 3D fields (unlimited)
34			%
35			%
36			% OUTPUTS:
37			% V : Volume of each CLASS (m3)
38			% Cm : Mean value of the classified field (allow to check errors)
39			% E : Each time a grid point is counted, a 1 is added to this 3D matrix
40			% Allow to check double count of a point or unexplored areas
41			% Vt : Is the total volume explored (Vt)
42			% CC : Contains the mean value of additional fields Ca, Cb ....
43			%
44			% NOTES:
45			% - Fields are on the format: C(DPT,LAT,LON)
46			% - The potential density is computed with the equation of state routine from
47			% the MITgcm called densjmd95.m
48			% (see: http://mitgcm.org/cgi-bin/viewcvs.cgi/MITgcm_contrib/gmaze_pv/subfct/densjmd95.m)
49			% - if dV is filled with NaN, dV is computed by the function
50			%
51			%
52			% AUTHOR:
53			% Guillaume Maze / MIT
54			%
55			% HISTORY:
56			% - Created: 06/29/2007
57			%
58
59			%
60
61			function varargout = diagVOLU(FLAG,C1,C2,CLASS,LON,LAT,DPT,DV,varargin)
62
63
64			% 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPROC
65			% Variables:
66			ndpt = size(C1,1);
67			nlat = size(C1,2);
68			nlon = size(C1,3);
69			CLASS = sort(CLASS(:));
70			[Z b c] = meshgrid(DPT,LON,LAT);clear b c, Z = permute(Z,[2 3 1]);
71
72			% Determine fields from which we'll take class contours:
73			switch FLAG
74
75			case 0 % Need to compute SIGMA THETA
76			THETA = C1;
77			SALT = C2;
78			ST = densjmd95(SALT,THETA,0.099989.81abs(Z)) - 1000;
79			CROP = ST;
80
81			case 1
82			ST = C1; % Potential density
83			CROP = ST;
84
85			case 2
86			THETA = C1; % Potential temperature
87			CROP = THETA;
88			end
89
90			% Volume elements:
91			if length(find(isnan(DV)==1)) == ndptnlatnlon
92			if exist('subfct_getdV','file')
93			DV = subfct_getdV(DPT,LAT,LON);
94			else
95			DV = local_getdV(DPT,LAT,LON);
96			end
97			end
98
99			% Need to compute volume integral over these 3D fields
100			nIN = nargin-8;
101			if nIN >= 1
102			doEXTRA = 1;
103			else
104			doEXTRA = 0;
105			end
106
107			% 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% VOLUME INTEGRATION
108			explored = zeros(ndpt,nlat,nlon);
109			% Volume integral:
110			for iC = 1 : length(CLASS)-1
111			mask = zeros(ndpt,nlat,nlon);
112			mask(find( (CLASS(iC) <= CROP) & (CROP < CLASS(iC+1)) )) = 1;
113			explored = explored + mask;
114			VOL(iC) = nansum(nansum(nansum(DV.*mask,1),2),3);
115
116			if VOL(iC) ~= 0
117			CAR(iC) = nansum(nansum(nansum(CROP.DV.mask,1),2),3)./VOL(iC);
118			if doEXTRA
119			for ii = 1 : nIN
120			C = varargin{ii};
121			CAREXTRA(ii,iC) = nansum(nansum(nansum(C.DV.mask,1),2),3)./VOL(iC);
122			end %for ii
123			end %if doEXTRA
124			else
125			CAR(iC) = NaN;
126			if doEXTRA
127			for ii = 1 : nIN
128			CAREXTRA(ii,iC) = NaN;
129			end %for ii
130			end %if doEXTRA
131			end
132			end %for iC
133
134			% In order to compute the total volume of the domain:
135			CROP(find(isnan(CROP)==0)) = 1;
136			CROP(find(isnan(CROP)==1)) = 0;
137			Vt = nansum(nansum(nansum(DV.*CROP,1),2),3);
138
139			% 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% OUTPUTS
140			switch nargout
141			case 1
142			varargout(1) = {VOL};
143			case 2
144			varargout(1) = {VOL};
145			varargout(2) = {CAR};
146			case 3
147			varargout(1) = {VOL};
148			varargout(2) = {CAR};
149			varargout(3) = {explored};
150			case 4
151			varargout(1) = {VOL};
152			varargout(2) = {CAR};
153			varargout(3) = {explored};
154			varargout(4) = {Vt};
155			case 5
156			varargout(1) = {VOL};
157			varargout(2) = {CAR};
158			varargout(3) = {explored};
159			varargout(4) = {Vt};
160			varargout(5) = {CAREXTRA};
161			end %switch
162
163
164
165
166
167
168
169			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
170			% This function computes the 3D dV volume elements.
171			% Copy of the subfct_getDV function from gmaze_pv package
172			function DV = local_getdV(Z,Y,X)
173
174			nz = length(Z);
175			ny = length(Y);
176			nx = length(X);
177
178			DV = zeros(nz,ny,nx);
179
180			% Vertical elements:
181			for iz = 1 : nz % Toward the deep ocean (because DPT<0)
182			% Vertical grid length centered at Z(iy)
183			if iz == 1
184			dz = abs(Z(1)) + abs(sum(diff(Z(iz:iz+1))/2));
185			elseif iz == nz % We don't know the real ocean depth
186			dz = abs(sum(diff(Z(iz-1:iz))/2));
187			else
188			dz = abs(sum(diff(Z(iz-1:iz+1))/2));
189			end
190			DZ(iz) = dz;
191			end
192
193			% Surface and Volume elements:
194			for ix = 1 : nx
195			for iy = 1 : ny
196			% Zonal grid length centered in X(ix),Y(iY)
197			if ix == 1
198			dx = abs(m_lldist([X(ix) X(ix+1)],[1 1]*Y(iy)))/2;
199			elseif ix == nx
200			dx = abs(m_lldist([X(ix-1) X(ix)],[1 1]*Y(iy)))/2;
201			else
202			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;
203			end
204
205			% Meridional grid length centered in X(ix),Y(iY)
206			if iy == 1
207			dy = abs(m_lldist([1 1]*X(ix),[Y(iy) Y(iy+1)]))/2;
208			elseif iy == ny
209			dy = abs(m_lldist([1 1]*X(ix),[Y(iy-1) Y(iy)]))/2;
210			else
211			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;
212			end
213
214			% Surface element:
215			DA = dxdy.ones(1,nz);
216
217			% Volume element:
218			DV(:,iy,ix) = DZ.*DA;
219			end %for iy
220			end %for ix
221