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 ../'])
