| 1 |
% [V,Cm,E,Vt,CC] = diagVOLUinterp(PII,B,dV,C,DX,DY,DZ,Tc,dT) |
| 2 |
% |
| 3 |
% DESCRIPTION: |
| 4 |
% Compute the volume of water embeded in Tc-dT/2 <= C < Tc+dT/2 |
| 5 |
% by interpolation of grid fraction |
| 6 |
% |
| 7 |
% INPUTS: |
| 8 |
% PII: 0/1 matrix defining the volume |
| 9 |
% B : 0/1 matrix defining the boundary's volume (from getVOLbounds) |
| 10 |
% dV : Volume elements matrix |
| 11 |
% C : Input field used to get PII, of size: ndpt x nlat x nlon |
| 12 |
% DX : zonal grid spacing of size: ndpt x nlat x nlon+1 |
| 13 |
% DY : meridional grid spacing of size: ndpt x nlat+1 x nlon |
| 14 |
% DZ : vertical grid spacing of size: ndpt+1 x nlat x nlon |
| 15 |
% Tc : the iso-C contour |
| 16 |
% dT : the bin used to get PII |
| 17 |
% |
| 18 |
% OUTPUTS: |
| 19 |
% V : Interpolated volume of the layer |
| 20 |
% Vraw : Raw volume, as returned by diagVOLU |
| 21 |
% |
| 22 |
% EG of use: |
| 23 |
% Tc = 18; dT = 2; |
| 24 |
% pii = boxcar(THETA,-10000,lon,lat,dpt,Tc,dT); |
| 25 |
% [BN BS BW BE BT BB]=getVOLbounds(pii); B=BN+BS+BE+BW+BT+BB; B(find(B~=0)) = 1; |
| 26 |
% [Vi Vr]=diagVOLUinterp(pii,B,dV_3D,THETA,DX,DY,DZ,Tc,dT) |
| 27 |
% |
| 28 |
% |
| 29 |
% AUTHOR: |
| 30 |
% Guillaume Maze / MIT |
| 31 |
% |
| 32 |
% HISTORY: |
| 33 |
% - Created: 07/24/2007 |
| 34 |
% |
| 35 |
|
| 36 |
% |
| 37 |
|
| 38 |
function varargout = diagVOLUinterp(pii,B,dV_3D,THETA,DX,DY,DZ,Tc,dT) |
| 39 |
|
| 40 |
ndpt = size(THETA,1); |
| 41 |
nlat = size(THETA,2); |
| 42 |
nlon = size(THETA,3); |
| 43 |
|
| 44 |
% Raw value of the volume: |
| 45 |
Vraw = nansum(nansum(nansum( pii.*dV_3D))); |
| 46 |
|
| 47 |
% Raw value without boundary points: |
| 48 |
Vraw_interior = nansum(nansum(nansum( (pii-B).*dV_3D))); |
| 49 |
|
| 50 |
ii = 0; |
| 51 |
npt = length(find(B==1)); |
| 52 |
dV = 0; |
| 53 |
dVraw = 0; |
| 54 |
warning off |
| 55 |
for iz = 1 : ndpt |
| 56 |
for ix = 1 : nlon |
| 57 |
for iy = 1 : nlat |
| 58 |
if B(iz,iy,ix) == 1 |
| 59 |
ii = ii + 1; |
| 60 |
%disp(sprintf('%d/%d/pii=%d',npt,ii,pii(iz,iy,ix))); |
| 61 |
Tloc = THETA(iz,iy,ix); |
| 62 |
clear T |
| 63 |
|
| 64 |
% Northern value: |
| 65 |
dyn = DY(iy+1,ix); c = 1/2; |
| 66 |
if iy+1 < nlat & isnan(THETA(iz,iy+1,ix)) == 0 |
| 67 |
T = THETA(iz,iy+1,ix); |
| 68 |
alphan = (dT/2 - abs(Tloc-Tc))./( abs(T-Tc) - abs(Tloc-Tc) ); |
| 69 |
if alphan > 1/2 & pii(iz,iy+1,ix) ~= 1 & ~isinf(alphan) |
| 70 |
c = alphan-1/2; |
| 71 |
end |
| 72 |
end |
| 73 |
dyn = dyn*c; |
| 74 |
|
| 75 |
% Southern value: |
| 76 |
dys = DY(iy,ix); c = 1/2; |
| 77 |
if iy-1 > 1 & isnan(THETA(iz,iy-1,ix)) == 0 |
| 78 |
T = THETA(iz,iy-1,ix); |
| 79 |
alphas = (dT/2 - abs(Tloc-Tc))./( abs(T-Tc) - abs(Tloc-Tc) ); |
| 80 |
if alphas > 1/2 & pii(iz,iy-1,ix) ~= 1 & ~isinf(alphas) |
| 81 |
c = alphas-1/2; |
| 82 |
end |
| 83 |
end |
| 84 |
dys = dys*c; |
| 85 |
|
| 86 |
% Eastern value: |
| 87 |
dxe = DX(iy,ix+1); c = 1/2; |
| 88 |
if ix+1 < nlon & isnan(THETA(iz,iy,ix+1)) == 0 |
| 89 |
T = THETA(iz,iy,ix+1); |
| 90 |
alphae = (dT/2 - abs(Tloc-Tc))./( abs(T-Tc) - abs(Tloc-Tc) ); |
| 91 |
if alphae > 1/2 & pii(iz,iy,ix+1) ~= 1 & ~isinf(alphae) |
| 92 |
c = alphae-1/2; |
| 93 |
end |
| 94 |
end |
| 95 |
dxe = dxe*c; |
| 96 |
|
| 97 |
% Western value: |
| 98 |
dxw = DX(iy,ix); c = 1/2; |
| 99 |
if ix-1 > 1 & isnan(THETA(iz,iy,ix-1)) == 0 |
| 100 |
T = THETA(iz,iy,ix-1); |
| 101 |
alphaw = (dT/2 - abs(Tloc-Tc))./( abs(T-Tc) - abs(Tloc-Tc) ); |
| 102 |
if alphaw > 1/2 & pii(iz,iy,ix-1) ~= 1 & ~isinf(alphaw) |
| 103 |
c = alphaw-1/2; |
| 104 |
end |
| 105 |
end |
| 106 |
dxw = dxw*c; |
| 107 |
|
| 108 |
% Top value: |
| 109 |
dzt = DZ(iz); c = 1/2; |
| 110 |
if iz-1 > 1 & isnan(THETA(iz-1,iy,ix)) == 0 |
| 111 |
T = THETA(iz-1,iy,ix); |
| 112 |
alphat = (dT/2 - abs(Tloc-Tc))./( abs(T-Tc) - abs(Tloc-Tc) ); |
| 113 |
if alphat > 1/2 & pii(iz-1,iy,ix) ~= 1 & ~isinf(alphat) |
| 114 |
c = alphat-1/2; |
| 115 |
end |
| 116 |
end |
| 117 |
dzt = dzt*c; |
| 118 |
|
| 119 |
% Bottom value: |
| 120 |
dzb = DZ(iz+1); c = 1/2; |
| 121 |
if iz+1 > 1 & isnan(THETA(iz+1,iy,ix)) == 0 |
| 122 |
T = THETA(iz+1,iy,ix); |
| 123 |
alphab = (dT/2 - abs(Tloc-Tc))./( abs(T-Tc) - abs(Tloc-Tc) ); |
| 124 |
if alphab > 1/2 & pii(iz+1,iy,ix) ~= 1 & ~isinf(alphab) |
| 125 |
c = alphab-1/2; |
| 126 |
end |
| 127 |
end |
| 128 |
dzb = dzb*c; |
| 129 |
|
| 130 |
dV(ii) = (dxw+dxe)*(dys+dyn)*(dzt+dzb); |
| 131 |
dVraw(ii) = dV_3D(iz,iy,ix); |
| 132 |
|
| 133 |
end %if boundary point |
| 134 |
end %for iy |
| 135 |
end %for ix |
| 136 |
end %for iz |
| 137 |
warning on |
| 138 |
Vraw2 = Vraw_interior + sum(dVraw); |
| 139 |
Vinterp = Vraw_interior + sum(dV); |
| 140 |
|
| 141 |
if Vraw2 ~= Vraw & 0 |
| 142 |
warning(sprintf('diagVOLUinterp: Raw volumes doesn''t match ! \n Difference in %%: %g',... |
| 143 |
(Vraw2-Vraw)/Vraw*100)) |
| 144 |
end |
| 145 |
|
| 146 |
% 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% OUTPUTS |
| 147 |
switch nargout |
| 148 |
case 1 |
| 149 |
varargout(1) = {Vinterp}; |
| 150 |
case 2 |
| 151 |
varargout(1) = {Vinterp}; |
| 152 |
varargout(2) = {Vraw2}; |
| 153 |
end %switch |
| 154 |
|
| 155 |
|
| 156 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph