| 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.09998*9.81*abs(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)) == ndpt*nlat*nlon |
| 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 = dx*dy.*ones(1,nz); |
| 216 |
|
| 217 |
% Volume element: |
| 218 |
DV(:,iy,ix) = DZ.*DA; |
| 219 |
end %for iy |
| 220 |
end %for ix |
| 221 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph