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

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

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


Revision 1.1 - (show 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 % [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