| 1 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 2 |
%% SUB-FUNCTIONS %% |
| 3 |
%% USED BY: VOLBET2ISO %% |
| 4 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 5 |
% USE DIRECTLY AT YOUR OWN RISK ! % |
| 6 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 7 |
|
| 8 |
|
| 9 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 10 |
% Master sub-function: |
| 11 |
function varargout = subfct_getvol(CHP,Z,Y,X,LIMITS) |
| 12 |
|
| 13 |
% Limits: |
| 14 |
%disp(strcat('Limits: ',num2str(LIMITS))); |
| 15 |
O = LIMITS(1); |
| 16 |
MZ = LIMITS(2); |
| 17 |
MY = sort( LIMITS(3:4) ); |
| 18 |
MX = sort( LIMITS(5:6) ); |
| 19 |
|
| 20 |
% Compute the volume: |
| 21 |
[V Vmat dV] = getvol(Z,Y,X,O,MZ,MY,MX,CHP); |
| 22 |
|
| 23 |
% Outputs: |
| 24 |
switch nargout |
| 25 |
case 1 |
| 26 |
varargout(1) = {V}; |
| 27 |
case 2 |
| 28 |
varargout(1) = {V}; |
| 29 |
varargout(2) = {Vmat}; |
| 30 |
case 3 |
| 31 |
varargout(1) = {V}; |
| 32 |
varargout(2) = {Vmat}; |
| 33 |
varargout(3) = {dV}; |
| 34 |
end %switch nargout |
| 35 |
|
| 36 |
|
| 37 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 38 |
% This function computes the volume limited southward by |
| 39 |
% MY(1), northward by iso-O (or MY(2) if iso-O reaches it), |
| 40 |
% eastward by MX(1), westward by MX(2) and downward by iso-O |
| 41 |
% (or MZ if iso-O reaches it). |
| 42 |
% |
| 43 |
% Updates: |
| 44 |
% 20060616: bug fixed in CHPsouth, CHPnorth nanmean |
| 45 |
% 20060615: add a default meridional gradient (negative) |
| 46 |
% test on x,y,z limits |
| 47 |
% |
| 48 |
function varargout = getvol(Z,Y,X,O,MZ,MY,MX,CHP) |
| 49 |
|
| 50 |
%% Dim: |
| 51 |
nz = length(Z); |
| 52 |
ny = length(Y); |
| 53 |
nx = length(X); |
| 54 |
%disp(num2str([nz ny nx])); |
| 55 |
|
| 56 |
%% Indices: |
| 57 |
izmax = min( find( Z>MZ ) ); if isempty(izmax),izmax=nz;end; |
| 58 |
iymin = max( find( Y<MY(1) ) ); if isempty(iymin),iymin=1; end; |
| 59 |
iymax = min( find( Y>MY(2) ) ); if isempty(iymax),iymax=ny;end; |
| 60 |
ixmin = max( find( X<MX(1) ) ); if isempty(ixmin),ixmin=1; end; |
| 61 |
ixmax = min( find( X>MX(2) ) ); if isempty(ixmax),ixmax=nx;end; |
| 62 |
%disp(num2str([1 izmax iymin iymax ixmin ixmax])); |
| 63 |
|
| 64 |
%% 1- determine the 3D matrix of volume elements defined by |
| 65 |
% the grid: |
| 66 |
dV = getdV(Z,Y,X); |
| 67 |
|
| 68 |
%% 2- compute the 3D volume matrix where 1 means dV must be |
| 69 |
% counted and 0 must not: |
| 70 |
|
| 71 |
V = ones(nz,ny,nx); % initialy keep all points |
| 72 |
|
| 73 |
% Exclude northward iso-O limits: |
| 74 |
% NB: here the test depends on the meridional gradient of CHP |
| 75 |
% if CHP increase (resp. decreases) with LAT then we must |
| 76 |
% keep lower (resp. higher) values than O limit |
| 77 |
% a: determine test type: |
| 78 |
N = iymax - iymin + 1; % Number of Y points in the domain |
| 79 |
CHPsouth = nanmean(nanmean(squeeze(CHP(1,iymin:iymin+fix(N/2),ixmin:ixmax)),1),2); |
| 80 |
CHPnorth = nanmean(nanmean(squeeze(CHP(1,iymin+fix(N/2):iymax,ixmin:ixmax)),1),2); |
| 81 |
SNgrad = (CHPnorth - CHPsouth)./abs(CHPnorth - CHPsouth); |
| 82 |
if isnan(SNgrad), SNgrad=-1; end % Assume negative gradient by default |
| 83 |
%disp(strcat('Northward gradient sign is:',num2str(SNgrad))); |
| 84 |
switch SNgrad |
| 85 |
case 1, testype = 'le'; % Less than or equal |
| 86 |
case -1, testype = 'ge'; % Greater than or equal |
| 87 |
end %switch |
| 88 |
% b: exclude points |
| 89 |
for iz=1:izmax |
| 90 |
chpZ = squeeze(CHP(iz,:,:)); |
| 91 |
V(iz,:,:)=double(feval(testype,chpZ,O)); |
| 92 |
end %for iz |
| 93 |
|
| 94 |
% Exclude southward limit: |
| 95 |
V(:,1:iymin,:) = zeros(nz,iymin,nx); |
| 96 |
|
| 97 |
% Exclude northward limit: |
| 98 |
V(:,iymax:ny,:) = zeros(nz,(ny-iymax)+1,nx); |
| 99 |
|
| 100 |
% Exclude westward limit: |
| 101 |
V(:,:,1:ixmin) = zeros(nz,ny,ixmin); |
| 102 |
|
| 103 |
% Exclude eastward limit: |
| 104 |
V(:,:,ixmax:nx) = zeros(nz,ny,(nx-ixmax)+1); |
| 105 |
|
| 106 |
% Exclude downward limit: |
| 107 |
V(izmax:nz,:,:) = zeros((nz-izmax)+1,ny,nx); |
| 108 |
|
| 109 |
|
| 110 |
%% 3- Then compute the volume by summing dV elements |
| 111 |
% for non 0 V points |
| 112 |
Vkeep = V.*dV; |
| 113 |
Vkeep = sum(sum(sum(Vkeep))); |
| 114 |
|
| 115 |
%% 4- Outputs: |
| 116 |
switch nargout |
| 117 |
case 1 |
| 118 |
varargout(1) = {Vkeep}; % Volume single value |
| 119 |
case 2 |
| 120 |
varargout(1) = {Vkeep}; |
| 121 |
varargout(2) = {V}; % Logical V matrix with included/excluded points |
| 122 |
case 3 |
| 123 |
varargout(1) = {Vkeep}; |
| 124 |
varargout(2) = {V}; |
| 125 |
varargout(3) = {dV}; % Volume elements matrix |
| 126 |
end %switch nargout |
| 127 |
|
| 128 |
|
| 129 |
|
| 130 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 131 |
% This function computes the 3D dV volume elements. |
| 132 |
function dv = getdV(Z,Y,X); |
| 133 |
|
| 134 |
nz = length(Z); |
| 135 |
ny = length(Y); |
| 136 |
nx = length(X); |
| 137 |
|
| 138 |
%%% Compute the DY: |
| 139 |
% Assuming Y is independant of ix: |
| 140 |
d = m_lldist([1 1]*X(1),Y); |
| 141 |
dy = [d(1)/2 (d(2:length(d))+d(1:length(d)-1))/2 d(length(d))/2]; |
| 142 |
dy = meshgrid(dy,X)'; |
| 143 |
|
| 144 |
%%% Compute the DX: |
| 145 |
clear d |
| 146 |
for iy = 1 : ny |
| 147 |
d(:,iy) = m_lldist(X,Y([iy iy])); |
| 148 |
end |
| 149 |
dx = [d(1,:)/2 ; ( d(2:size(d,1),:) + d(1:size(d,1)-1,:) )./2 ; d(size(d,1),:)/2]; |
| 150 |
dx = dx'; |
| 151 |
|
| 152 |
%% Compute the horizontal DS surface element: |
| 153 |
ds = dx.*dy; |
| 154 |
|
| 155 |
%% Compute the DZ: |
| 156 |
d = diff(Z); |
| 157 |
dz = [d(1)/2 ( d(2:length(d)) + d(1:length(d)-1) )./2 d(length(d))/2]; |
| 158 |
|
| 159 |
%% Then compute the DV volume elements: |
| 160 |
dv = ones(nz,ny,nx)*NaN; |
| 161 |
for iz=1:nz |
| 162 |
dv(iz,:,:) = dz(iz).*ds; |
| 163 |
end |
| 164 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph