Back to home page

MITgcm

 
 

    


Warning, /verification/fizhi-cs-aqualev20/scripts/APEprocess1.m is written in an unsupported language. File is not indexed.

view on githubraw file Latest commit 7ae8fb32 on 2006-04-03 20:55:29 UTC
7ae8fb32b5 Andr*0001 function [APEprocess1] = APEprocess1(dirsuffix)
                0002 datadir=['mnc_out_',dirsuffix];
                0003 eval(['cd ' datadir])
                0004 % Load grid fields
                0005 griddata = rdmnc('grid.*');
                0006  XC = griddata.XC;
                0007  YC = griddata.YC;
                0008  XG = griddata.XG; XG = XG(1:end-1,1:end-1);
                0009  YG = griddata.YG; YG = YG(1:end-1,1:end-1);
                0010  dxC = griddata.dxC;
                0011  dyC = griddata.dyC;
                0012  dxG = griddata.dxG;
                0013  dyG = griddata.dyG;
                0014  RAC = griddata.rA;
                0015  HFacC = griddata.HFacC;
                0016  HFacW = griddata.HFacW;
                0017  HFacS = griddata.HFacS;
                0018  ZC = griddata.RC;
                0019  ZF = griddata.RF;
                0020  drC = griddata.drC;
                0021  drF = griddata.drF;
                0022 
                0023 % Load ML fields
                0024 mldata = rdmnc('MLfields.*');
                0025 ntimes = size(mldata.T,1);
                0026 temp = mldata.THETA(:,:,:,1);
                0027 nx = size(temp,1);
                0028 ny = size(temp,2);
                0029 nPg = nx*ny;
                0030 nr = size(temp,3);
                0031 % Reference geopotential profile for 17 levels:
                0032 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];
                0033 % load COS & SIN of rotation angle:
                0034 fid=fopen('/u/molod/proj_cs32_2uEvN.bin','r','b'); uvEN=fread(fid,nx*ny*2,'real*8'); fclose(fid);
                0035 uvEN=reshape(uvEN,[nPg 2]);
                0036 AngleCS=uvEN(:,1);
                0037 AngleSN=uvEN(:,2);
                0038 % Some initialization of arrays
                0039 uav=zeros(nx,ny,nr);
                0040 vav=zeros(nx,ny,nr);
                0041 wav=zeros(nx,ny,nr);
                0042 thetaav=zeros(nx,ny,nr);
                0043 saltav=zeros(nx,ny,nr);
                0044 phiav=zeros(nx,ny,nr);
                0045 hfactot=zeros(nx,ny,nr);
                0046 uzav=zeros(64,nr);
                0047 vzav=zeros(64,nr);
                0048 wzav=zeros(64,nr);
                0049 thetazav=zeros(64,nr);
                0050 saltzav=zeros(64,nr);
                0051 phizav=zeros(64,nr);
                0052 usqz=zeros(64,nr);
                0053 vsqz=zeros(64,nr);
                0054 thetasqz=zeros(64,nr);
                0055 phisqz=zeros(64,nr);
                0056 omegasqz=zeros(64,nr);
                0057 qsqz=zeros(64,nr);
                0058 uvz=zeros(64,nr);
                0059 uomegaz=zeros(64,nr);
                0060 vomegaz=zeros(64,nr);
                0061 vthetaz=zeros(64,nr);
                0062 omegathetaz=zeros(64,nr);
                0063 vqz=zeros(64,nr);
                0064 omegaqz=zeros(64,nr);
                0065 vphiz=zeros(64,nr);
                0066 hfacztot=zeros(64,nr);
                0067 for it=1:ntimes
                0068 % 1) compute zonal mean of scalars (theta,q,phi,omega) using calc_ZonAv_CS:
                0069 temp = mldata.THETA(:,:,:,it);
                0070 hfac4calcz = ones(nx,ny,nr);
                0071 hfac4calcz(find(temp(:,:,:)==0))=0.;
                0072 hfactot = hfactot + hfac4calcz;
                0073 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',hfac4calcz);
                0074 thetaz=tempZ(:,:,1);
                0075 hfacz = ones(64,nr);
                0076 hfacz(:,1) = dum1(:,1) ./ dum1(:,17);
                0077 hfacztot = hfacztot + hfacz;
                0078 temp = mldata.SALT(:,:,:,it);
                0079 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',hfac4calcz);
                0080 saltz=tempZ(:,:,1);
                0081 temp = mldata.PHIHYD(:,:,:,it);
                0082 % Add reference geopotential to anomaly
                0083 for ilev = 1:nr
                0084 ind2d = find(hfac4calcz(:,:,ilev)~=0);
                0085 atemp = squeeze(temp(:,:,ilev));
                0086 atemp(ind2d)=atemp(ind2d) + phiref(ilev);
                0087 temp(:,:,ilev) = atemp;
                0088 end
                0089 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',hfac4calcz);
                0090 phiz=tempZ(:,:,1);
                0091 temp = mldata.WVEL(:,:,:,it);
                0092 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',hfac4calcz);
                0093 wz=tempZ(:,:,1);
                0094 % 2) rotate u and v then compute zonal mean using calc_ZonAv_CS:
                0095 u = mldata.UVEL(:,:,:,it);
                0096 u = u(1:nx,1:ny,:);
                0097 v = mldata.VVEL(:,:,:,it);
                0098 v = v(1:nx,1:ny,:);
                0099 [uE,vN] = rotate_uv2uvEN(u,v,AngleCS,AngleSN);
                0100 temp = uE;
                0101 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',hfac4calcz);
                0102 uz=tempZ(:,:,1);
                0103 temp = vN;
                0104 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',hfac4calcz);
                0105 vz=tempZ(:,:,1);
                0106 % 3)compute the products and accumulate into a time average:
                0107 ind00=find(hfacz~=0);
                0108 usqz(ind00)=usqz(ind00) + uz(ind00).*uz(ind00).*hfacz(ind00);
                0109 vsqz(ind00)=vsqz(ind00) + vz(ind00).*vz(ind00).*hfacz(ind00);
                0110 thetasqz(ind00)=thetasqz(ind00) + thetaz(ind00).*thetaz(ind00).*hfacz(ind00);
                0111 phisqz(ind00)=phisqz(ind00) + phiz(ind00).*phiz(ind00).*hfacz(ind00);
                0112 omegasqz(ind00)=omegasqz(ind00) + wz(ind00).*wz(ind00).*hfacz(ind00);
                0113 qsqz(ind00)=qsqz(ind00) + saltz(ind00).*saltz(ind00).*hfacz(ind00);
                0114 uvz(ind00)=uvz(ind00) + uz(ind00).*vz(ind00).*hfacz(ind00);
                0115 uomegaz(ind00)=uomegaz(ind00) + uz(ind00).*wz(ind00).*hfacz(ind00);
                0116 vomegaz(ind00)=vomegaz(ind00) + vz(ind00).*wz(ind00).*hfacz(ind00);
                0117 vthetaz(ind00)=vthetaz(ind00) + vz(ind00).*thetaz(ind00).*hfacz(ind00);
                0118 omegathetaz(ind00)=omegathetaz(ind00) + wz(ind00).*thetaz(ind00).*hfacz(ind00);
                0119 vqz(ind00)=vqz(ind00) + vz(ind00).*saltz(ind00).*hfacz(ind00);
                0120 omegaqz(ind00)=omegaqz(ind00) + wz(ind00).*saltz(ind00).*hfacz(ind00);
                0121 vphiz(ind00)=vphiz(ind00) + vz(ind00).*phiz(ind00).*hfacz(ind00);
                0122 %
                0123 % Compute time average of zonal means:
                0124 uzav(ind00)=uzav(ind00) + uz(ind00).*hfacz(ind00);
                0125 vzav(ind00)=vzav(ind00) + vz(ind00).*hfacz(ind00);
                0126 wzav(ind00)=wzav(ind00) + wz(ind00).*hfacz(ind00);
                0127 thetazav(ind00)=thetazav(ind00) + thetaz(ind00).*hfacz(ind00);
                0128 saltzav(ind00)=saltzav(ind00) + saltz(ind00).*hfacz(ind00);
                0129 phizav(ind00)=phizav(ind00) + phiz(ind00).*hfacz(ind00);
                0130 %
                0131 % 4)Compute time average of u, v, Theta, q, Omega, Phi (full fields)
                0132 ind0=find(hfac4calcz~=0);
                0133 uav(ind0)=uav(ind0) + u(ind0);
                0134 vav(ind0)=vav(ind0) + v(ind0);
                0135 ttemp=squeeze(mldata.WVEL(:,:,:,it));
                0136 wav(ind0)=wav(ind0) + ttemp(ind0);
                0137 ttemp=squeeze(mldata.THETA(:,:,:,it));
                0138 thetaav(ind0)=thetaav(ind0) + ttemp(ind0);
                0139 ttemp=squeeze(mldata.SALT(:,:,:,it));
                0140 saltav(ind0)=saltav(ind0) + ttemp(ind0);
                0141 ttemp=squeeze(mldata.PHIHYD(:,:,:,it));
                0142 phiav(ind0)=phiav(ind0) + ttemp(ind0);
                0143 %
                0144 end
                0145 % Divide the sums by the appropriate counter to get time average
                0146 ind1=find(hfacztot~=0);
                0147 usqz(ind1)=usqz(ind1) ./ hfacztot(ind1);
                0148 vsqz(ind1)=vsqz(ind1) ./ hfacztot(ind1);
                0149 thetasqz(ind1)=thetasqz(ind1) ./ hfacztot(ind1);
                0150 phisqz(ind1)=phisqz(ind1) ./ hfacztot(ind1);
                0151 omegasqz(ind1)=omegasqz(ind1) ./ hfacztot(ind1);
                0152 qsqz(ind1)=qsqz(ind1) ./ hfacztot(ind1);
                0153 uvz(ind1)=uvz(ind1) ./ hfacztot(ind1);
                0154 uomegaz(ind1)=uomegaz(ind1) ./ hfacztot(ind1);
                0155 vomegaz(ind1)=vomegaz(ind1) ./ hfacztot(ind1);
                0156 vthetaz(ind1)=vthetaz(ind1) ./ hfacztot(ind1);
                0157 omegathetaz(ind1)=omegathetaz(ind1) ./ hfacztot(ind1);
                0158 vqz(ind1)=vqz(ind1) ./ hfacztot(ind1);
                0159 omegaqz(ind1)=omegaqz(ind1) ./ hfacztot(ind1);
                0160 vphiz(ind1)=vphiz(ind1) ./ hfacztot(ind1);
                0161 uzav(ind1)=uzav(ind1) ./ hfacztot(ind1);
                0162 vzav(ind1)=vzav(ind1) ./ hfacztot(ind1);
                0163 wzav(ind1)=wzav(ind1) ./ hfacztot(ind1);
                0164 thetazav(ind1)=thetazav(ind1) ./ hfacztot(ind1);
                0165 saltzav(ind1)=saltzav(ind1) ./ hfacztot(ind1);
                0166 phizav(ind1)=phizav(ind1) ./ hfacztot(ind1);
                0167 %
                0168 ind2=find(hfactot~=0);
                0169 uav(ind2)=uav(ind2) ./ hfactot(ind2);
                0170 vav(ind2)=vav(ind2) ./ hfactot(ind2);
                0171 wav(ind2)=wav(ind2) ./ hfactot(ind2);
                0172 thetaav(ind2)=thetaav(ind2) ./ hfactot(ind2);
                0173 saltav(ind2)=saltav(ind2) ./ hfactot(ind2);
                0174 phiav(ind2)=phiav(ind2) ./ hfactot(ind2);
                0175 %
                0176 % Add reference geopotential to anomaly
                0177 for ilev = 1:nr
                0178 atemp = squeeze(phiav(:,:,ilev));
                0179 ind2d2 = find(hfactot(:,:,ilev)~=0);
                0180 atemp(ind2d2)=atemp(ind2d2) + phiref(ilev);
                0181 phiav(:,:,ilev) = atemp;
                0182 end
                0183 %
                0184 % Release memory of all time-dependant ML quantities
                0185 clear mldata
                0186 %
                0187 % Some steps that have to wait until all 10-day periods are processed
                0188 % 5) Compute squared time averages:
                0189 % 6) Compute zonal mean of squared time averaged:
                0190 %
                0191 % Write out (time averaged) zonal means of first moments
                0192 %
                0193 fid = fopen('uzon.interim','w','b'); fwrite(fid,uzav,'real*8'); fclose(fid);
                0194 fid = fopen('vzon.interim','w','b'); fwrite(fid,vzav,'real*8'); fclose(fid);
                0195 fid = fopen('wzon.interim','w','b'); fwrite(fid,wzav,'real*8'); fclose(fid);
                0196 fid = fopen('thetazon.interim','w','b'); fwrite(fid,thetazav,'real*8'); fclose(fid);
                0197 fid = fopen('saltzon.interim','w','b'); fwrite(fid,saltzav,'real*8'); fclose(fid);
                0198 fid = fopen('phizon.interim','w','b'); fwrite(fid,phizav,'real*8'); fclose(fid);
                0199 weightoutz = hfacztot ./ ntimes;
                0200 fid = fopen('hfacztot1.interim','w','b'); fwrite(fid,weightoutz,'real*8'); fclose(fid);
                0201 %
                0202 % Write out (time averaged) products of zonal means
                0203 %
                0204 fid = fopen('usqz.interim','w','b'); fwrite(fid,usqz,'real*8'); fclose(fid);
                0205 fid = fopen('vsqz.interim','w','b'); fwrite(fid,vsqz,'real*8'); fclose(fid);
                0206 fid = fopen('thetasqz.interim','w','b'); fwrite(fid,thetasqz,'real*8'); fclose(fid);
                0207 fid = fopen('phisqz.interim','w','b'); fwrite(fid,phisqz,'real*8'); fclose(fid);
                0208 fid = fopen('omegasqz.interim','w','b'); fwrite(fid,omegasqz,'real*8'); fclose(fid);
                0209 fid = fopen('qsqz.interim','w','b'); fwrite(fid,qsqz,'real*8'); fclose(fid);
                0210 fid = fopen('uvz.interim','w','b'); fwrite(fid,uvz,'real*8'); fclose(fid);
                0211 fid = fopen('uomegaz.interim','w','b'); fwrite(fid,uomegaz,'real*8'); fclose(fid);
                0212 fid = fopen('vomegaz.interim','w','b'); fwrite(fid,vomegaz,'real*8'); fclose(fid);
                0213 fid = fopen('vthetaz.interim','w','b'); fwrite(fid,vthetaz,'real*8'); fclose(fid);
                0214 fid = fopen('omegathetaz.interim','w','b'); fwrite(fid,omegathetaz,'real*8'); fclose(fid);
                0215 fid = fopen('vqz.interim','w','b'); fwrite(fid,vqz,'real*8'); fclose(fid);
                0216 fid = fopen('omegaqz.interim','w','b'); fwrite(fid,omegaqz,'real*8'); fclose(fid);
                0217 fid = fopen('vphiz.interim','w','b'); fwrite(fid,vphiz,'real*8'); fclose(fid);
                0218 %
                0219 % Write out time averages of first moments
                0220 %
                0221 fid = fopen('uave.interim','w','b'); fwrite(fid,uav,'real*8'); fclose(fid);
                0222 fid = fopen('vave.interim','w','b'); fwrite(fid,vav,'real*8'); fclose(fid);
                0223 fid = fopen('wave.interim','w','b'); fwrite(fid,wav,'real*8'); fclose(fid);
                0224 fid = fopen('thetaave.interim','w','b'); fwrite(fid,thetaav,'real*8'); fclose(fid);
                0225 fid = fopen('saltave.interim','w','b'); fwrite(fid,saltav,'real*8'); fclose(fid);
                0226 fid = fopen('phiave.interim','w','b'); fwrite(fid,phiav,'real*8'); fclose(fid);
                0227 weightout = hfactot ./ ntimes;
                0228 fid = fopen('hfactot1.interim','w','b'); fwrite(fid,weightout,'real*8'); fclose(fid);
                0229 %
                0230 %Load MF fields
                0231 mfdata = rdmnc('MFfields.*');
                0232 ntimes = size(mfdata.T,1);
                0233 % allocate some space for time averages
                0234 uvelsq = zeros(nx,ny,nr);
                0235 vvelsq = zeros(nx,ny,nr);
                0236 thetasq = zeros(nx,ny,nr);
                0237 wvelsq = zeros(nx,ny,nr);
                0238 phihydsq = zeros(nx,ny,nr);
                0239 saltsq = zeros(nx,ny,nr);
                0240 uvvelc = zeros(nx,ny,nr);
                0241 wuvel = zeros(nx,ny,nr);
                0242 wvvel = zeros(nx,ny,nr);
                0243 uvelth = zeros(nx,ny,nr);
                0244 vvelth = zeros(nx,ny,nr);
                0245 wvelth = zeros(nx,ny,nr);
                0246 uvelslt = zeros(nx,ny,nr);
                0247 vvelslt = zeros(nx,ny,nr);
                0248 wvelslt = zeros(nx,ny,nr);
                0249 uvelphi = zeros(nx,ny,nr);
                0250 vvelphi = zeros(nx,ny,nr);
                0251 hfacztot=zeros(64,nr);
                0252 %
                0253 %1)Compute time averages of quantities - cannot use mean function because of undefs
                0254 %      also - variances can't be negative (u**2,v**2, w**2, q**2, t**2 and phi**2)
                0255 for it=1:ntimes
                0256 temp00 = mfdata.THETASQ(1:nx,1:ny,:,it);
                0257 hfac4calcz = ones(nx,ny,nr);
                0258 hfac4calcz(find(temp00(:,:,:)==0))=0.;
                0259 [ind5]=find(hfac4calcz~=0);
                0260 temp0 = mfdata.UVELSQ(1:nx,1:ny,:,it);
                0261 temp0(find(temp0<0.))=0.;
                0262 uvelsq(ind5) = uvelsq(ind5) + temp0(ind5);
                0263 temp0 = mfdata.VVELSQ(1:nx,1:ny,:,it);
                0264 temp0(find(temp0<0.))=0.;
                0265 vvelsq(ind5) = vvelsq(ind5) + temp0(ind5);
                0266 temp0 = mfdata.THETASQ(1:nx,1:ny,:,it);
                0267 temp0(find(temp0<0.))=0.;
                0268 thetasq(ind5) = thetasq(ind5) + temp0(ind5);
                0269 temp0 = mfdata.WVELSQ(1:nx,1:ny,:,it);
                0270 temp0(find(temp0<0.))=0.;
                0271 wvelsq(ind5) = wvelsq(ind5) + temp0(ind5);
                0272 temp0 = mfdata.PHIHYDSQ(1:nx,1:ny,:,it);
                0273 temp0(find(temp0<0.))=0.;
                0274 phihydsq(ind5) = phihydsq(ind5) + temp0(ind5);
                0275 temp0 = mfdata.SALTSQ(1:nx,1:ny,:,it);
                0276 temp0(find(temp0<0.))=0.;
                0277 saltsq(ind5) = saltsq(ind5) + temp0(ind5);
                0278 temp0 = mfdata.UV_VEL_C(1:nx,1:ny,:,it);
                0279 uvvelc(ind5) = uvvelc(ind5) + temp0(ind5);
                0280 temp0 = mfdata.WU_VEL(1:nx,1:ny,:,it);
                0281 wuvel(ind5) = wuvel(ind5) + temp0(ind5);
                0282 temp0 = mfdata.WV_VEL(1:nx,1:ny,:,it);
                0283 wvvel(ind5) = wvvel(ind5) + temp0(ind5);
                0284 temp0 = mfdata.UVELTH(1:nx,1:ny,:,it);
                0285 uvelth(ind5) = uvelth(ind5) + temp0(ind5);
                0286 temp0 = mfdata.VVELTH(1:nx,1:ny,:,it);
                0287 vvelth(ind5) = vvelth(ind5) + temp0(ind5);
                0288 temp0 = mfdata.WVELTH(1:nx,1:ny,:,it);
                0289 wvelth(ind5) = wvelth(ind5) + temp0(ind5);
                0290 temp0 = mfdata.UVELSLT(1:nx,1:ny,:,it);
                0291 uvelslt(ind5) = uvelslt(ind5) + temp0(ind5);
                0292 temp0 = mfdata.VVELSLT(1:nx,1:ny,:,it);
                0293 vvelslt(ind5) = vvelslt(ind5) + temp0(ind5);
                0294 temp0 = mfdata.WVELSLT(1:nx,1:ny,:,it);
                0295 wvelslt(ind5) = wvelslt(ind5) + temp0(ind5);
                0296 temp0 = mfdata.UVELPHI(1:nx,1:ny,:,it);
                0297 uvelphi(ind5) = uvelphi(ind5) + temp0(ind5);
                0298 temp0 = mfdata.VVELPHI(1:nx,1:ny,:,it);
                0299 vvelphi(ind5) = vvelphi(ind5) + temp0(ind5);
                0300 end
                0301 %Release memory for ML time-dependant quantities
                0302 clear mfdata
                0303 % Divide the sums by the appropriate counter to get time average
                0304 %
                0305 ind3=find(hfactot~=0);
                0306 uvelsq(ind3) = uvelsq(ind3) ./ hfactot(ind3);
                0307 vvelsq(ind3) = vvelsq(ind3) ./ hfactot(ind3);
                0308 thetasq(ind3) = thetasq(ind3) ./ hfactot(ind3);
                0309 wvelsq(ind3) = wvelsq(ind3) ./ hfactot(ind3);
                0310 phihydsq(ind3) = phihydsq(ind3) ./ hfactot(ind3);
                0311 saltsq(ind3) = saltsq(ind3) ./ hfactot(ind3);
                0312 uvvelc(ind3) = uvvelc(ind3) ./ hfactot(ind3);
                0313 wuvel(ind3) = wuvel(ind3) ./ hfactot(ind3);
                0314 wvvel(ind3) = wvvel(ind3) ./ hfactot(ind3);
                0315 uvelth(ind3) = uvelth(ind3) ./ hfactot(ind3);
                0316 vvelth(ind3) = vvelth(ind3) ./ hfactot(ind3);
                0317 wvelth(ind3) = wvelth(ind3) ./ hfactot(ind3);
                0318 uvelslt(ind3) = uvelslt(ind3) ./ hfactot(ind3);
                0319 vvelslt(ind3) = vvelslt(ind3) ./ hfactot(ind3);
                0320 wvelslt(ind3) = wvelslt(ind3) ./ hfactot(ind3);
                0321 uvelphi(ind3) = uvelphi(ind3) ./ hfactot(ind3);
                0322 vvelphi(ind3) = vvelphi(ind3) ./ hfactot(ind3);
                0323 %2)Use calc2ndmom to compute UV, U squared and V squared
                0324 [UVtot,UVtrans,U2tot,U2trans,V2tot,V2trans,errFlag]=calc2ndmom(uav,vav,uvelsq,vvelsq,uvvelc,AngleCS,AngleSN);
                0325 %3)Move to centers and rotate moments related to u and v transports:
                0326 [wuE,wvN] = rotate_uv2uvEN(wuvel,wvvel,AngleCS,AngleSN);
                0327 [uEth,vNth] = rotate_uv2uvEN(uvelth,vvelth,AngleCS,AngleSN);
                0328 [uEslt,vNslt] = rotate_uv2uvEN(uvelslt,vvelslt,AngleCS,AngleSN);
                0329 [uEphi,vNphi] = rotate_uv2uvEN(uvelphi,vvelphi,AngleCS,AngleSN);
                0330 %4) Compute zonal means:
                0331 temp = UVtot;
                0332 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout);
                0333 UVtotz=tempZ(:,:,1);
                0334 temp = U2tot;
                0335 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout);
                0336 U2totz=tempZ(:,:,1);
                0337 temp = V2tot;
                0338 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout);
                0339 V2totz=tempZ(:,:,1);
                0340 temp = wvelsq;
                0341 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout);
                0342 wvelsqz=tempZ(:,:,1);
                0343 temp = thetasq;
                0344 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout);
                0345 potempsqz=tempZ(:,:,1);
                0346 temp = phihydsq;
                0347 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout);
                0348 phihydsqz=tempZ(:,:,1);
                0349 temp = saltsq;
                0350 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout);
                0351 saltsqz=tempZ(:,:,1);
                0352 temp = wuE;
                0353 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout);
                0354 wuEz=tempZ(:,:,1);
                0355 temp = wvN;
                0356 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout);
                0357 wvNz=tempZ(:,:,1);
                0358 temp = vNth;
                0359 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout);
                0360 vNthz=tempZ(:,:,1);
                0361 temp = wvelth;
                0362 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout);
                0363 wvelthz=tempZ(:,:,1);
                0364 temp = vNslt;
                0365 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout);
                0366 vNsltz=tempZ(:,:,1);
                0367 temp = wvelslt;
                0368 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout);
                0369 wvelsltz=tempZ(:,:,1);
                0370 temp = vNphi;
                0371 [tempZ,dum1,dum2] = calc_ZonAv_CS(temp,1,1,0,XC,YC,XG,YG,RAC,'./',weightout);
                0372 vNphiz=tempZ(:,:,1);
                0373 %
                0374 weightoutz2 = ones(64,nr);
                0375 weightoutz2(:,1) = dum1(:,1) ./ dum1(:,17);
                0376 %
                0377 % Add reference geopotential to anomaly
                0378 for ilev = 1:nr
                0379 ind1d = find(weightoutz(:,ilev)~=0);
                0380 phihydsqz(ind1d,ilev)=phihydsqz(ind1d,ilev) - phiref(ilev).*phiref(ilev) + 2.*phiref(ilev).*phizav(ind1d,ilev);
                0381 vNphiz(ind1d,ilev)=vNphiz(ind1d,ilev) + phiref(ilev).*vzav(ind1d,ilev);
                0382 end
                0383 %
                0384 % Write out time averages of zonal mean second moments
                0385 fid = fopen('uvtotz.interim','w','b'); fwrite(fid,UVtotz,'real*8'); fclose(fid);
                0386 fid = fopen('u2totz.interim','w','b'); fwrite(fid,U2totz,'real*8'); fclose(fid);
                0387 fid = fopen('v2totz.interim','w','b'); fwrite(fid,V2totz,'real*8'); fclose(fid);
                0388 fid = fopen('wvelsqz.interim','w','b'); fwrite(fid,wvelsqz,'real*8'); fclose(fid);
                0389 fid = fopen('potempsqz.interim','w','b'); fwrite(fid,potempsqz,'real*8'); fclose(fid);
                0390 fid = fopen('phihydsqz.interim','w','b'); fwrite(fid,phihydsqz,'real*8'); fclose(fid);
                0391 fid = fopen('saltsqz.interim','w','b'); fwrite(fid,saltsqz,'real*8'); fclose(fid);
                0392 fid = fopen('wuEz.interim','w','b'); fwrite(fid,wuEz,'real*8'); fclose(fid);
                0393 fid = fopen('wvNz.interim','w','b'); fwrite(fid,wvNz,'real*8'); fclose(fid);
                0394 fid = fopen('vNthz.interim','w','b'); fwrite(fid,vNthz,'real*8'); fclose(fid);
                0395 fid = fopen('wvelthz.interim','w','b'); fwrite(fid,wvelthz,'real*8'); fclose(fid);
                0396 fid = fopen('vNsltz.interim','w','b'); fwrite(fid,vNsltz,'real*8'); fclose(fid);
                0397 fid = fopen('wvelsltz.interim','w','b'); fwrite(fid,wvelsltz,'real*8'); fclose(fid);
                0398 fid = fopen('vNphiz.interim','w','b'); fwrite(fid,vNphiz,'real*8'); fclose(fid);
                0399 fid = fopen('hfacztot2.interim','w','b'); fwrite(fid,weightoutz2,'real*8'); fclose(fid);
                0400 %
                0401 % Last step has to wait until all 10-day peices are processed: Compute terms
                0402 %
                0403 eval(['cd ../'])