1 |
function [APEprocess1] = APEprocess1(dirsuffix) |
2 |
datadir=['mnc_out_',dirsuffix]; |
3 |
eval(['cd ' datadir]) |
4 |
% Load grid fields |
5 |
griddata = rdmnc('grid.*'); |
6 |
XC = griddata.XC; |
7 |
YC = griddata.YC; |
8 |
XG = griddata.XG; XG = XG(1:end-1,1:end-1); |
9 |
YG = griddata.YG; YG = YG(1:end-1,1:end-1); |
10 |
dxC = griddata.dxC; |
11 |
dyC = griddata.dyC; |
12 |
dxG = griddata.dxG; |
13 |
dyG = griddata.dyG; |
14 |
RAC = griddata.rA; |
15 |
HFacC = griddata.HFacC; |
16 |
HFacW = griddata.HFacW; |
17 |
HFacS = griddata.HFacS; |
18 |
ZC = griddata.RC; |
19 |
ZF = griddata.RF; |
20 |
drC = griddata.drC; |
21 |
drF = griddata.drF; |
22 |
|
23 |
% Load ML fields |
24 |
mldata = rdmnc('MLfields.*'); |
25 |
ntimes = size(mldata.T,1); |
26 |
temp = mldata.THETA(:,:,:,1); |
27 |
nx = size(temp,1); |
28 |
ny = size(temp,2); |
29 |
nPg = nx*ny; |
30 |
nr = size(temp,3); |
31 |
% Reference geopotential profile for 17 levels: |
32 |
phiref=9.8.*[ -0.739878953 646.302002 1338.38696 2894.67627 4099.63135 5484.93359 7116.62549 9115.08008 10321.4688 11741.3574 13494.4102 15862.4404 17833.5605 19667.8887 22136.1934 23854.2266 26366.9375]; |
33 |
% load COS & SIN of rotation angle: |
34 |
fid=fopen('/u/molod/proj_cs32_2uEvN.bin','r','b'); uvEN=fread(fid,nx*ny*2,'real*8'); fclose(fid); |
35 |
uvEN=reshape(uvEN,[nPg 2]); |
36 |
AngleCS=uvEN(:,1); |
37 |
AngleSN=uvEN(:,2); |
38 |
% Some initialization of arrays |
39 |
uav=zeros(nx,ny,nr); |
40 |
vav=zeros(nx,ny,nr); |
41 |
wav=zeros(nx,ny,nr); |
42 |
thetaav=zeros(nx,ny,nr); |
43 |
saltav=zeros(nx,ny,nr); |
44 |
phiav=zeros(nx,ny,nr); |
45 |
hfactot=zeros(nx,ny,nr); |
46 |
uzav=zeros(64,nr); |
47 |
vzav=zeros(64,nr); |
48 |
wzav=zeros(64,nr); |
49 |
thetazav=zeros(64,nr); |
50 |
saltzav=zeros(64,nr); |
51 |
phizav=zeros(64,nr); |
52 |
usqz=zeros(64,nr); |
53 |
vsqz=zeros(64,nr); |
54 |
thetasqz=zeros(64,nr); |
55 |
phisqz=zeros(64,nr); |
56 |
omegasqz=zeros(64,nr); |
57 |
qsqz=zeros(64,nr); |
58 |
uvz=zeros(64,nr); |
59 |
uomegaz=zeros(64,nr); |
60 |
vomegaz=zeros(64,nr); |
61 |
vthetaz=zeros(64,nr); |
62 |
omegathetaz=zeros(64,nr); |
63 |
vqz=zeros(64,nr); |
64 |
omegaqz=zeros(64,nr); |
65 |
vphiz=zeros(64,nr); |
66 |
hfacztot=zeros(64,nr); |
67 |
for it=1:ntimes |
68 |
% 1) compute zonal mean of scalars (theta,q,phi,omega) using calc_ZonAv_CS: |
69 |
temp = mldata.THETA(:,:,:,it); |
70 |
hfac4calcz = ones(nx,ny,nr); |
71 |
hfac4calcz(find(temp(:,:,:)==0))=0.; |
72 |
hfactot = hfactot + hfac4calcz; |
73 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',hfac4calcz); |
74 |
thetaz=tempZ(:,:,1); |
75 |
hfacz = ones(64,nr); |
76 |
hfacz(:,1) = dum1(:,1) ./ dum1(:,17); |
77 |
hfacztot = hfacztot + hfacz; |
78 |
temp = mldata.SALT(:,:,:,it); |
79 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',hfac4calcz); |
80 |
saltz=tempZ(:,:,1); |
81 |
temp = mldata.PHIHYD(:,:,:,it); |
82 |
% Add reference geopotential to anomaly |
83 |
for ilev = 1:nr |
84 |
ind2d = find(hfac4calcz(:,:,ilev)~=0); |
85 |
atemp = squeeze(temp(:,:,ilev)); |
86 |
atemp(ind2d)=atemp(ind2d) + phiref(ilev); |
87 |
temp(:,:,ilev) = atemp; |
88 |
end |
89 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',hfac4calcz); |
90 |
phiz=tempZ(:,:,1); |
91 |
temp = mldata.WVEL(:,:,:,it); |
92 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',hfac4calcz); |
93 |
wz=tempZ(:,:,1); |
94 |
% 2) rotate u and v then compute zonal mean using calc_ZonAv_CS: |
95 |
u = mldata.UVEL(:,:,:,it); |
96 |
u = u(1:nx,1:ny,:); |
97 |
v = mldata.VVEL(:,:,:,it); |
98 |
v = v(1:nx,1:ny,:); |
99 |
[uE,vN] = rotate_uv2uvEN(u,v,AngleCS,AngleSN); |
100 |
temp = uE; |
101 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',hfac4calcz); |
102 |
uz=tempZ(:,:,1); |
103 |
temp = vN; |
104 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',hfac4calcz); |
105 |
vz=tempZ(:,:,1); |
106 |
% 3)compute the products and accumulate into a time average: |
107 |
ind00=find(hfacz~=0); |
108 |
usqz(ind00)=usqz(ind00) + uz(ind00).*uz(ind00).*hfacz(ind00); |
109 |
vsqz(ind00)=vsqz(ind00) + vz(ind00).*vz(ind00).*hfacz(ind00); |
110 |
thetasqz(ind00)=thetasqz(ind00) + thetaz(ind00).*thetaz(ind00).*hfacz(ind00); |
111 |
phisqz(ind00)=phisqz(ind00) + phiz(ind00).*phiz(ind00).*hfacz(ind00); |
112 |
omegasqz(ind00)=omegasqz(ind00) + wz(ind00).*wz(ind00).*hfacz(ind00); |
113 |
qsqz(ind00)=qsqz(ind00) + saltz(ind00).*saltz(ind00).*hfacz(ind00); |
114 |
uvz(ind00)=uvz(ind00) + uz(ind00).*vz(ind00).*hfacz(ind00); |
115 |
uomegaz(ind00)=uomegaz(ind00) + uz(ind00).*wz(ind00).*hfacz(ind00); |
116 |
vomegaz(ind00)=vomegaz(ind00) + vz(ind00).*wz(ind00).*hfacz(ind00); |
117 |
vthetaz(ind00)=vthetaz(ind00) + vz(ind00).*thetaz(ind00).*hfacz(ind00); |
118 |
omegathetaz(ind00)=omegathetaz(ind00) + wz(ind00).*thetaz(ind00).*hfacz(ind00); |
119 |
vqz(ind00)=vqz(ind00) + vz(ind00).*saltz(ind00).*hfacz(ind00); |
120 |
omegaqz(ind00)=omegaqz(ind00) + wz(ind00).*saltz(ind00).*hfacz(ind00); |
121 |
vphiz(ind00)=vphiz(ind00) + vz(ind00).*phiz(ind00).*hfacz(ind00); |
122 |
% |
123 |
% Compute time average of zonal means: |
124 |
uzav(ind00)=uzav(ind00) + uz(ind00).*hfacz(ind00); |
125 |
vzav(ind00)=vzav(ind00) + vz(ind00).*hfacz(ind00); |
126 |
wzav(ind00)=wzav(ind00) + wz(ind00).*hfacz(ind00); |
127 |
thetazav(ind00)=thetazav(ind00) + thetaz(ind00).*hfacz(ind00); |
128 |
saltzav(ind00)=saltzav(ind00) + saltz(ind00).*hfacz(ind00); |
129 |
phizav(ind00)=phizav(ind00) + phiz(ind00).*hfacz(ind00); |
130 |
% |
131 |
% 4)Compute time average of u, v, Theta, q, Omega, Phi (full fields) |
132 |
ind0=find(hfac4calcz~=0); |
133 |
uav(ind0)=uav(ind0) + u(ind0); |
134 |
vav(ind0)=vav(ind0) + v(ind0); |
135 |
ttemp=squeeze(mldata.WVEL(:,:,:,it)); |
136 |
wav(ind0)=wav(ind0) + ttemp(ind0); |
137 |
ttemp=squeeze(mldata.THETA(:,:,:,it)); |
138 |
thetaav(ind0)=thetaav(ind0) + ttemp(ind0); |
139 |
ttemp=squeeze(mldata.SALT(:,:,:,it)); |
140 |
saltav(ind0)=saltav(ind0) + ttemp(ind0); |
141 |
ttemp=squeeze(mldata.PHIHYD(:,:,:,it)); |
142 |
phiav(ind0)=phiav(ind0) + ttemp(ind0); |
143 |
% |
144 |
end |
145 |
% Divide the sums by the appropriate counter to get time average |
146 |
ind1=find(hfacztot~=0); |
147 |
usqz(ind1)=usqz(ind1) ./ hfacztot(ind1); |
148 |
vsqz(ind1)=vsqz(ind1) ./ hfacztot(ind1); |
149 |
thetasqz(ind1)=thetasqz(ind1) ./ hfacztot(ind1); |
150 |
phisqz(ind1)=phisqz(ind1) ./ hfacztot(ind1); |
151 |
omegasqz(ind1)=omegasqz(ind1) ./ hfacztot(ind1); |
152 |
qsqz(ind1)=qsqz(ind1) ./ hfacztot(ind1); |
153 |
uvz(ind1)=uvz(ind1) ./ hfacztot(ind1); |
154 |
uomegaz(ind1)=uomegaz(ind1) ./ hfacztot(ind1); |
155 |
vomegaz(ind1)=vomegaz(ind1) ./ hfacztot(ind1); |
156 |
vthetaz(ind1)=vthetaz(ind1) ./ hfacztot(ind1); |
157 |
omegathetaz(ind1)=omegathetaz(ind1) ./ hfacztot(ind1); |
158 |
vqz(ind1)=vqz(ind1) ./ hfacztot(ind1); |
159 |
omegaqz(ind1)=omegaqz(ind1) ./ hfacztot(ind1); |
160 |
vphiz(ind1)=vphiz(ind1) ./ hfacztot(ind1); |
161 |
uzav(ind1)=uzav(ind1) ./ hfacztot(ind1); |
162 |
vzav(ind1)=vzav(ind1) ./ hfacztot(ind1); |
163 |
wzav(ind1)=wzav(ind1) ./ hfacztot(ind1); |
164 |
thetazav(ind1)=thetazav(ind1) ./ hfacztot(ind1); |
165 |
saltzav(ind1)=saltzav(ind1) ./ hfacztot(ind1); |
166 |
phizav(ind1)=phizav(ind1) ./ hfacztot(ind1); |
167 |
% |
168 |
ind2=find(hfactot~=0); |
169 |
uav(ind2)=uav(ind2) ./ hfactot(ind2); |
170 |
vav(ind2)=vav(ind2) ./ hfactot(ind2); |
171 |
wav(ind2)=wav(ind2) ./ hfactot(ind2); |
172 |
thetaav(ind2)=thetaav(ind2) ./ hfactot(ind2); |
173 |
saltav(ind2)=saltav(ind2) ./ hfactot(ind2); |
174 |
phiav(ind2)=phiav(ind2) ./ hfactot(ind2); |
175 |
% |
176 |
% Add reference geopotential to anomaly |
177 |
for ilev = 1:nr |
178 |
atemp = squeeze(phiav(:,:,ilev)); |
179 |
ind2d2 = find(hfactot(:,:,ilev)~=0); |
180 |
atemp(ind2d2)=atemp(ind2d2) + phiref(ilev); |
181 |
phiav(:,:,ilev) = atemp; |
182 |
end |
183 |
% |
184 |
% Release memory of all time-dependant ML quantities |
185 |
clear mldata |
186 |
% |
187 |
% Some steps that have to wait until all 10-day periods are processed |
188 |
% 5) Compute squared time averages: |
189 |
% 6) Compute zonal mean of squared time averaged: |
190 |
% |
191 |
% Write out (time averaged) zonal means of first moments |
192 |
% |
193 |
fid = fopen('uzon.interim','w','b'); fwrite(fid,uzav,'real*8'); fclose(fid); |
194 |
fid = fopen('vzon.interim','w','b'); fwrite(fid,vzav,'real*8'); fclose(fid); |
195 |
fid = fopen('wzon.interim','w','b'); fwrite(fid,wzav,'real*8'); fclose(fid); |
196 |
fid = fopen('thetazon.interim','w','b'); fwrite(fid,thetazav,'real*8'); fclose(fid); |
197 |
fid = fopen('saltzon.interim','w','b'); fwrite(fid,saltzav,'real*8'); fclose(fid); |
198 |
fid = fopen('phizon.interim','w','b'); fwrite(fid,phizav,'real*8'); fclose(fid); |
199 |
weightoutz = hfacztot ./ ntimes; |
200 |
fid = fopen('hfacztot1.interim','w','b'); fwrite(fid,weightoutz,'real*8'); fclose(fid); |
201 |
% |
202 |
% Write out (time averaged) products of zonal means |
203 |
% |
204 |
fid = fopen('usqz.interim','w','b'); fwrite(fid,usqz,'real*8'); fclose(fid); |
205 |
fid = fopen('vsqz.interim','w','b'); fwrite(fid,vsqz,'real*8'); fclose(fid); |
206 |
fid = fopen('thetasqz.interim','w','b'); fwrite(fid,thetasqz,'real*8'); fclose(fid); |
207 |
fid = fopen('phisqz.interim','w','b'); fwrite(fid,phisqz,'real*8'); fclose(fid); |
208 |
fid = fopen('omegasqz.interim','w','b'); fwrite(fid,omegasqz,'real*8'); fclose(fid); |
209 |
fid = fopen('qsqz.interim','w','b'); fwrite(fid,qsqz,'real*8'); fclose(fid); |
210 |
fid = fopen('uvz.interim','w','b'); fwrite(fid,uvz,'real*8'); fclose(fid); |
211 |
fid = fopen('uomegaz.interim','w','b'); fwrite(fid,uomegaz,'real*8'); fclose(fid); |
212 |
fid = fopen('vomegaz.interim','w','b'); fwrite(fid,vomegaz,'real*8'); fclose(fid); |
213 |
fid = fopen('vthetaz.interim','w','b'); fwrite(fid,vthetaz,'real*8'); fclose(fid); |
214 |
fid = fopen('omegathetaz.interim','w','b'); fwrite(fid,omegathetaz,'real*8'); fclose(fid); |
215 |
fid = fopen('vqz.interim','w','b'); fwrite(fid,vqz,'real*8'); fclose(fid); |
216 |
fid = fopen('omegaqz.interim','w','b'); fwrite(fid,omegaqz,'real*8'); fclose(fid); |
217 |
fid = fopen('vphiz.interim','w','b'); fwrite(fid,vphiz,'real*8'); fclose(fid); |
218 |
% |
219 |
% Write out time averages of first moments |
220 |
% |
221 |
fid = fopen('uave.interim','w','b'); fwrite(fid,uav,'real*8'); fclose(fid); |
222 |
fid = fopen('vave.interim','w','b'); fwrite(fid,vav,'real*8'); fclose(fid); |
223 |
fid = fopen('wave.interim','w','b'); fwrite(fid,wav,'real*8'); fclose(fid); |
224 |
fid = fopen('thetaave.interim','w','b'); fwrite(fid,thetaav,'real*8'); fclose(fid); |
225 |
fid = fopen('saltave.interim','w','b'); fwrite(fid,saltav,'real*8'); fclose(fid); |
226 |
fid = fopen('phiave.interim','w','b'); fwrite(fid,phiav,'real*8'); fclose(fid); |
227 |
weightout = hfactot ./ ntimes; |
228 |
fid = fopen('hfactot1.interim','w','b'); fwrite(fid,weightout,'real*8'); fclose(fid); |
229 |
% |
230 |
%Load MF fields |
231 |
mfdata = rdmnc('MFfields.*'); |
232 |
ntimes = size(mfdata.T,1); |
233 |
% allocate some space for time averages |
234 |
uvelsq = zeros(nx,ny,nr); |
235 |
vvelsq = zeros(nx,ny,nr); |
236 |
thetasq = zeros(nx,ny,nr); |
237 |
wvelsq = zeros(nx,ny,nr); |
238 |
phihydsq = zeros(nx,ny,nr); |
239 |
saltsq = zeros(nx,ny,nr); |
240 |
uvvelc = zeros(nx,ny,nr); |
241 |
wuvel = zeros(nx,ny,nr); |
242 |
wvvel = zeros(nx,ny,nr); |
243 |
uvelth = zeros(nx,ny,nr); |
244 |
vvelth = zeros(nx,ny,nr); |
245 |
wvelth = zeros(nx,ny,nr); |
246 |
uvelslt = zeros(nx,ny,nr); |
247 |
vvelslt = zeros(nx,ny,nr); |
248 |
wvelslt = zeros(nx,ny,nr); |
249 |
uvelphi = zeros(nx,ny,nr); |
250 |
vvelphi = zeros(nx,ny,nr); |
251 |
hfacztot=zeros(64,nr); |
252 |
% |
253 |
%1)Compute time averages of quantities - cannot use mean function because of undefs |
254 |
% also - variances can't be negative (u**2,v**2, w**2, q**2, t**2 and phi**2) |
255 |
for it=1:ntimes |
256 |
temp00 = mfdata.THETASQ(1:nx,1:ny,:,it); |
257 |
hfac4calcz = ones(nx,ny,nr); |
258 |
hfac4calcz(find(temp00(:,:,:)==0))=0.; |
259 |
[ind5]=find(hfac4calcz~=0); |
260 |
temp0 = mfdata.UVELSQ(1:nx,1:ny,:,it); |
261 |
temp0(find(temp0<0.))=0.; |
262 |
uvelsq(ind5) = uvelsq(ind5) + temp0(ind5); |
263 |
temp0 = mfdata.VVELSQ(1:nx,1:ny,:,it); |
264 |
temp0(find(temp0<0.))=0.; |
265 |
vvelsq(ind5) = vvelsq(ind5) + temp0(ind5); |
266 |
temp0 = mfdata.THETASQ(1:nx,1:ny,:,it); |
267 |
temp0(find(temp0<0.))=0.; |
268 |
thetasq(ind5) = thetasq(ind5) + temp0(ind5); |
269 |
temp0 = mfdata.WVELSQ(1:nx,1:ny,:,it); |
270 |
temp0(find(temp0<0.))=0.; |
271 |
wvelsq(ind5) = wvelsq(ind5) + temp0(ind5); |
272 |
temp0 = mfdata.PHIHYDSQ(1:nx,1:ny,:,it); |
273 |
temp0(find(temp0<0.))=0.; |
274 |
phihydsq(ind5) = phihydsq(ind5) + temp0(ind5); |
275 |
temp0 = mfdata.SALTSQ(1:nx,1:ny,:,it); |
276 |
temp0(find(temp0<0.))=0.; |
277 |
saltsq(ind5) = saltsq(ind5) + temp0(ind5); |
278 |
temp0 = mfdata.UV_VEL_C(1:nx,1:ny,:,it); |
279 |
uvvelc(ind5) = uvvelc(ind5) + temp0(ind5); |
280 |
temp0 = mfdata.WU_VEL(1:nx,1:ny,:,it); |
281 |
wuvel(ind5) = wuvel(ind5) + temp0(ind5); |
282 |
temp0 = mfdata.WV_VEL(1:nx,1:ny,:,it); |
283 |
wvvel(ind5) = wvvel(ind5) + temp0(ind5); |
284 |
temp0 = mfdata.UVELTH(1:nx,1:ny,:,it); |
285 |
uvelth(ind5) = uvelth(ind5) + temp0(ind5); |
286 |
temp0 = mfdata.VVELTH(1:nx,1:ny,:,it); |
287 |
vvelth(ind5) = vvelth(ind5) + temp0(ind5); |
288 |
temp0 = mfdata.WVELTH(1:nx,1:ny,:,it); |
289 |
wvelth(ind5) = wvelth(ind5) + temp0(ind5); |
290 |
temp0 = mfdata.UVELSLT(1:nx,1:ny,:,it); |
291 |
uvelslt(ind5) = uvelslt(ind5) + temp0(ind5); |
292 |
temp0 = mfdata.VVELSLT(1:nx,1:ny,:,it); |
293 |
vvelslt(ind5) = vvelslt(ind5) + temp0(ind5); |
294 |
temp0 = mfdata.WVELSLT(1:nx,1:ny,:,it); |
295 |
wvelslt(ind5) = wvelslt(ind5) + temp0(ind5); |
296 |
temp0 = mfdata.UVELPHI(1:nx,1:ny,:,it); |
297 |
uvelphi(ind5) = uvelphi(ind5) + temp0(ind5); |
298 |
temp0 = mfdata.VVELPHI(1:nx,1:ny,:,it); |
299 |
vvelphi(ind5) = vvelphi(ind5) + temp0(ind5); |
300 |
end |
301 |
%Release memory for ML time-dependant quantities |
302 |
clear mfdata |
303 |
% Divide the sums by the appropriate counter to get time average |
304 |
% |
305 |
ind3=find(hfactot~=0); |
306 |
uvelsq(ind3) = uvelsq(ind3) ./ hfactot(ind3); |
307 |
vvelsq(ind3) = vvelsq(ind3) ./ hfactot(ind3); |
308 |
thetasq(ind3) = thetasq(ind3) ./ hfactot(ind3); |
309 |
wvelsq(ind3) = wvelsq(ind3) ./ hfactot(ind3); |
310 |
phihydsq(ind3) = phihydsq(ind3) ./ hfactot(ind3); |
311 |
saltsq(ind3) = saltsq(ind3) ./ hfactot(ind3); |
312 |
uvvelc(ind3) = uvvelc(ind3) ./ hfactot(ind3); |
313 |
wuvel(ind3) = wuvel(ind3) ./ hfactot(ind3); |
314 |
wvvel(ind3) = wvvel(ind3) ./ hfactot(ind3); |
315 |
uvelth(ind3) = uvelth(ind3) ./ hfactot(ind3); |
316 |
vvelth(ind3) = vvelth(ind3) ./ hfactot(ind3); |
317 |
wvelth(ind3) = wvelth(ind3) ./ hfactot(ind3); |
318 |
uvelslt(ind3) = uvelslt(ind3) ./ hfactot(ind3); |
319 |
vvelslt(ind3) = vvelslt(ind3) ./ hfactot(ind3); |
320 |
wvelslt(ind3) = wvelslt(ind3) ./ hfactot(ind3); |
321 |
uvelphi(ind3) = uvelphi(ind3) ./ hfactot(ind3); |
322 |
vvelphi(ind3) = vvelphi(ind3) ./ hfactot(ind3); |
323 |
%2)Use calc2ndmom to compute UV, U squared and V squared |
324 |
[UVtot,UVtrans,U2tot,U2trans,V2tot,V2trans,errFlag]=calc2ndmom(uav,vav,uvelsq,vvelsq,uvvelc,AngleCS,AngleSN); |
325 |
%3)Move to centers and rotate moments related to u and v transports: |
326 |
[wuE,wvN] = rotate_uv2uvEN(wuvel,wvvel,AngleCS,AngleSN); |
327 |
[uEth,vNth] = rotate_uv2uvEN(uvelth,vvelth,AngleCS,AngleSN); |
328 |
[uEslt,vNslt] = rotate_uv2uvEN(uvelslt,vvelslt,AngleCS,AngleSN); |
329 |
[uEphi,vNphi] = rotate_uv2uvEN(uvelphi,vvelphi,AngleCS,AngleSN); |
330 |
%4) Compute zonal means: |
331 |
temp = UVtot; |
332 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout); |
333 |
UVtotz=tempZ(:,:,1); |
334 |
temp = U2tot; |
335 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout); |
336 |
U2totz=tempZ(:,:,1); |
337 |
temp = V2tot; |
338 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout); |
339 |
V2totz=tempZ(:,:,1); |
340 |
temp = wvelsq; |
341 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout); |
342 |
wvelsqz=tempZ(:,:,1); |
343 |
temp = thetasq; |
344 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout); |
345 |
potempsqz=tempZ(:,:,1); |
346 |
temp = phihydsq; |
347 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout); |
348 |
phihydsqz=tempZ(:,:,1); |
349 |
temp = saltsq; |
350 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout); |
351 |
saltsqz=tempZ(:,:,1); |
352 |
temp = wuE; |
353 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout); |
354 |
wuEz=tempZ(:,:,1); |
355 |
temp = wvN; |
356 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout); |
357 |
wvNz=tempZ(:,:,1); |
358 |
temp = vNth; |
359 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout); |
360 |
vNthz=tempZ(:,:,1); |
361 |
temp = wvelth; |
362 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout); |
363 |
wvelthz=tempZ(:,:,1); |
364 |
temp = vNslt; |
365 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout); |
366 |
vNsltz=tempZ(:,:,1); |
367 |
temp = wvelslt; |
368 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout); |
369 |
wvelsltz=tempZ(:,:,1); |
370 |
temp = vNphi; |
371 |
[tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout); |
372 |
vNphiz=tempZ(:,:,1); |
373 |
% |
374 |
weightoutz2 = ones(64,nr); |
375 |
weightoutz2(:,1) = dum1(:,1) ./ dum1(:,17); |
376 |
% |
377 |
% Add reference geopotential to anomaly |
378 |
for ilev = 1:nr |
379 |
ind1d = find(weightoutz(:,ilev)~=0); |
380 |
phihydsqz(ind1d,ilev)=phihydsqz(ind1d,ilev) - phiref(ilev).*phiref(ilev) + 2.*phiref(ilev).*phizav(ind1d,ilev); |
381 |
vNphiz(ind1d,ilev)=vNphiz(ind1d,ilev) + phiref(ilev).*vzav(ind1d,ilev); |
382 |
end |
383 |
% |
384 |
% Write out time averages of zonal mean second moments |
385 |
fid = fopen('uvtotz.interim','w','b'); fwrite(fid,UVtotz,'real*8'); fclose(fid); |
386 |
fid = fopen('u2totz.interim','w','b'); fwrite(fid,U2totz,'real*8'); fclose(fid); |
387 |
fid = fopen('v2totz.interim','w','b'); fwrite(fid,V2totz,'real*8'); fclose(fid); |
388 |
fid = fopen('wvelsqz.interim','w','b'); fwrite(fid,wvelsqz,'real*8'); fclose(fid); |
389 |
fid = fopen('potempsqz.interim','w','b'); fwrite(fid,potempsqz,'real*8'); fclose(fid); |
390 |
fid = fopen('phihydsqz.interim','w','b'); fwrite(fid,phihydsqz,'real*8'); fclose(fid); |
391 |
fid = fopen('saltsqz.interim','w','b'); fwrite(fid,saltsqz,'real*8'); fclose(fid); |
392 |
fid = fopen('wuEz.interim','w','b'); fwrite(fid,wuEz,'real*8'); fclose(fid); |
393 |
fid = fopen('wvNz.interim','w','b'); fwrite(fid,wvNz,'real*8'); fclose(fid); |
394 |
fid = fopen('vNthz.interim','w','b'); fwrite(fid,vNthz,'real*8'); fclose(fid); |
395 |
fid = fopen('wvelthz.interim','w','b'); fwrite(fid,wvelthz,'real*8'); fclose(fid); |
396 |
fid = fopen('vNsltz.interim','w','b'); fwrite(fid,vNsltz,'real*8'); fclose(fid); |
397 |
fid = fopen('wvelsltz.interim','w','b'); fwrite(fid,wvelsltz,'real*8'); fclose(fid); |
398 |
fid = fopen('vNphiz.interim','w','b'); fwrite(fid,vNphiz,'real*8'); fclose(fid); |
399 |
fid = fopen('hfacztot2.interim','w','b'); fwrite(fid,weightoutz2,'real*8'); fclose(fid); |
400 |
% |
401 |
% Last step has to wait until all 10-day peices are processed: Compute terms |
402 |
% |
403 |
eval(['cd ../']) |