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