/[MITgcm]/MITgcm_contrib/gmaze_pv/subfct/diagVOLUinterp.m
ViewVC logotype

Annotation of /MITgcm_contrib/gmaze_pv/subfct/diagVOLUinterp.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (hide annotations) (download)
Wed Sep 19 15:24:41 2007 UTC (17 years, 10 months ago) by gmaze
Branch: MAIN
CVS Tags: HEAD
General Update

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    

  ViewVC Help
Powered by ViewVC 1.1.22