high_res_cube/eddy_flux/calcUVBave_ttcp.m

function m = calcUVBave_ttcp(Spat, Tpat, Upat, Vpat, ilist, nlay, ne)

% Function m = calcTBave(Spat, Tpat, Upat, Vpat, ilist, nlay, ne)
% 
%
% INPUTS
% Tpat     string containing T-data file pattern
% Upat     string containing U-data file pattern
% Vpat     string containing V-data file pattern
% ilist    list of iteration values
% nlay     number of z layers
% ne       number of values alone each edge of each face
%
% OUTPUTS
% m        output array of dimension ns*nps

if nargin ~= 7
  disp('Error: wrong number of arguments')
  exit
end

deltatT = 1200;                % model time step size (s)
tavefreq = 259200;             % averaging period (s)
startDate = datenum(1992,1,1); % model integration starting date

oldyear = -1;
newyear = oldyear;

nps = 6 * ne*ne;
tx = 85;
ty = 85;
nt = ne*ne*6/(tx*ty);
cx = ne;
cy = ne;
nz = 1;

delR   = [ 
    10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ...
    10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ...
    31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ...
    93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ...
    139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ...
    341.50,364.50,387.50,410.50,433.50,456.50 ];
R = cumsum(delR) - 0.5*delR;

vnam = [ 'U_'; 'V_'; 'T_'; 'S_'; 'B_'; 'B1'; ...
         'UT'; 'VT'; 'US'; 'VS'; 'UB'; 'VB' ];

lcstr = [ '! ( rm -f %s ; ttcp -r > %s ) > /dev/null 2>&1 &' ];
rcstr = [ '!ssh 10.12.2.200 ''( dd bs=6242400 skip=%d ' ...
          ' count=1 if=/home/edhill/%s' ...
          ' | ttcp -t 10.12.2.41 ) > /dev/null 2>&1''' ];

disp(['lcstr = "' lcstr '"']);
disp(['rcstr = "' rcstr '"']);

il = 1;
for il = 1:nlay
  
  a = sprintf('Slice: %g',il);
  %  disp(a);
  
  %  Zero the accumulators
  nm = 0;
  for io = 1:size(vnam,1)
    comm = sprintf('%s_acc = zeros(ne,ne,6);', vnam(io,:));
    eval(comm)
  end
  
  it = 1;
  for it = 1:length(ilist)
    
    iter = ilist(it);
    str = datestr(startDate + (iter*deltatT-tavefreq/2)/60/60/24);
    dv = datevec(startDate + (iter*deltatT-tavefreq/2)/60/60/24);
    newyear = dv(1);
    
    sfname = sprintf(Spat, ilist(it));
    tfname = sprintf(Tpat, ilist(it));
    ufname = sprintf(Upat, ilist(it));
    vfname = sprintf(Vpat, ilist(it));
    disp([ a '  ' str '  ' tfname ]);
    skip = il - 1;
    eval( sprintf(lcstr, sfname, sfname) );
    eval( sprintf(rcstr, skip,   sfname) );
    eval( sprintf(lcstr, tfname, tfname) );
    eval( sprintf(rcstr, skip,   tfname) );
    eval( sprintf(lcstr, ufname, ufname) );
    eval( sprintf(rcstr, skip,   ufname) );
    eval( sprintf(lcstr, vfname, vfname) );
    eval( sprintf(rcstr, skip,   vfname) );
    sfid = fopen(sfname, 'r', 'ieee-be');
    tfid = fopen(tfname, 'r', 'ieee-be');
    ufid = fopen(ufname, 'r', 'ieee-be');
    vfid = fopen(vfname, 'r', 'ieee-be');
    offset = (il - 1)*nps*4;
    
    %  Convert the horrible JPL mdsio layout to six sane faces
    %fseek(sfid, offset, 'bof');
    scube = unmangleJPL1( reshape(fread(sfid, nps, 'real*4'),tx*nt,ty), ...
                          ne, tx );
    
    %fseek(tfid, offset, 'bof');
    tcube = unmangleJPL1( reshape(fread(tfid, nps, 'real*4'),tx*nt,ty), ...
                          ne, tx );

    %fseek(ufid, offset, 'bof');
    ucube = unmangleJPL1( reshape(fread(ufid, nps, 'real*4'),tx*nt,ty), ...
                          ne, tx );
    
    %fseek(vfid, offset, 'bof');
    vcube = unmangleJPL1( reshape(fread(vfid, nps, 'real*4'),tx*nt,ty), ...
                          ne, tx );

    fclose(sfid);
    fclose(tfid);
    fclose(ufid);
    fclose(vfid);
    
    %  Compute the density
    bcube   = zeros(size(scube));
    bcubem1 = zeros(size(scube));
    for icf = 1:6
      % pressure_for_eos.F :
      %   locPres(i,j) = -rhoConst*rC(k)*gravity*maskC(i,j,k,bi,bj)
      %   rhoConst = 9.997999999999999E+02  /* Ref density ( kg/m^3 ) */
      %   gravity  = 9.810000000000000E+00  /* Grav accel ( m/s^2 ) */
      p = 0.09998 * 9.81 * R(il);
      bcube(:,:,icf) = ...
          eh3densjmd95( squeeze(scube(:,:,icf)), ...
                        squeeze(tcube(:,:,icf)), p ) * -9.81 / 1000;
      if il > 1
        p = 0.09998 * 9.81 * R(il-1);
        bcubem1(:,:,icf) = ...
            eh3densjmd95( scube(:,:,icf), ...
                          tcube(:,:,icf), p ) * -9.81 / 1000;
      end
    end
    tmask   = double(tcube ~= 0.0);
    scube   = scube .* tmask;
    bcube   = bcube .* tmask;
    bcubem1 = bcubem1 .* tmask;

    %  Fill the (ne+1)*(ne+1) grid with T,S,B values
    ni = ne;    nip1 = ni + 1;
    nj = ne;    njp1 = nj + 1;
    tcubep1 = zeros( [ nip1 njp1 6 ] );
    scubep1 = zeros( [ nip1 njp1 6 ] );
    bcubep1 = zeros( [ nip1 njp1 6 ] );
    for icf = 1:6
      tcubep1(2:nip1,2:njp1,icf) = tcube(:,:,icf);
      scubep1(2:nip1,2:njp1,icf) = scube(:,:,icf);
      bcubep1(2:nip1,2:njp1,icf) = bcube(:,:,icf);
    end

    %  Do the upwind-edge T exchanges
    tcubep1(1,2:njp1,1) = tcube(ni:-1:1,nj,5);   % - 
    tcubep1(2:nip1,1,1) = tcube(1:ni,nj,6);      % -
    tcubep1(1,2:njp1,2) = tcube(ni,1:nj,1);      % -
    tcubep1(2:nip1,1,2) = tcube(ni,nj:-1:1,6);   % -
    tcubep1(1,2:njp1,3) = tcube(ni:-1:1,nj,1);   % -
    tcubep1(2:nip1,1,3) = tcube(1:ni,nj,2);      % -
    tcubep1(1,2:njp1,4) = tcube(ni,1:nj,3);      % -
    tcubep1(2:nip1,1,4) = tcube(ni,nj:-1:1,2);   % -
    tcubep1(1,2:njp1,5) = tcube(ni:-1:1,nj,3);   % -
    tcubep1(2:nip1,1,5) = tcube(1:ni,nj,4);      % -
    tcubep1(1,2:njp1,6) = tcube(ni,1:nj,5);      % -
    tcubep1(2:nip1,1,6) = tcube(ni,nj:-1:1,4);   % -

    %  Do the upwind-edge S exchanges
    scubep1(1,2:njp1,1) = scube(ni:-1:1,nj,5);   % - 
    scubep1(2:nip1,1,1) = scube(1:ni,nj,6);      % -
    scubep1(1,2:njp1,2) = scube(ni,1:nj,1);      % -
    scubep1(2:nip1,1,2) = scube(ni,nj:-1:1,6);   % -
    scubep1(1,2:njp1,3) = scube(ni:-1:1,nj,1);   % -
    scubep1(2:nip1,1,3) = scube(1:ni,nj,2);      % -
    scubep1(1,2:njp1,4) = scube(ni,1:nj,3);      % -
    scubep1(2:nip1,1,4) = scube(ni,nj:-1:1,2);   % -
    scubep1(1,2:njp1,5) = scube(ni:-1:1,nj,3);   % -
    scubep1(2:nip1,1,5) = scube(1:ni,nj,4);      % -
    scubep1(1,2:njp1,6) = scube(ni,1:nj,5);      % -
    scubep1(2:nip1,1,6) = scube(ni,nj:-1:1,4);   % -

    %  Do the upwind-edge B exchanges
    bcubep1(1,2:njp1,1) = bcube(ni:-1:1,nj,5);   % - 
    bcubep1(2:nip1,1,1) = bcube(1:ni,nj,6);      % -
    bcubep1(1,2:njp1,2) = bcube(ni,1:nj,1);      % -
    bcubep1(2:nip1,1,2) = bcube(ni,nj:-1:1,6);   % -
    bcubep1(1,2:njp1,3) = bcube(ni:-1:1,nj,1);   % -
    bcubep1(2:nip1,1,3) = bcube(1:ni,nj,2);      % -
    bcubep1(1,2:njp1,4) = bcube(ni,1:nj,3);      % -
    bcubep1(2:nip1,1,4) = bcube(ni,nj:-1:1,2);   % -
    bcubep1(1,2:njp1,5) = bcube(ni:-1:1,nj,3);   % -
    bcubep1(2:nip1,1,5) = bcube(1:ni,nj,4);      % -
    bcubep1(1,2:njp1,6) = bcube(ni,1:nj,5);      % -
    bcubep1(2:nip1,1,6) = bcube(ni,nj:-1:1,4);   % -

    %  Get T values on the U,V grid points
    masku = 1.0 - abs(diff(double(tcubep1(2:nip1,:,:) == 0.0),1,2));
    maskv = 1.0 - abs(diff(double(tcubep1(:,2:nip1,:) == 0.0),1,1));
    diffu = 0.5*diff(tcubep1(2:nip1,:,:),1,2);
    diffv = 0.5*diff(tcubep1(:,2:nip1,:),1,1);
    tonu = tcube(:,:,:) + masku.*diffu;
    tonv = tcube(:,:,:) + maskv.*diffv;

    %  Get S values on the U,V grid points
    diffu = 0.5*diff(scubep1(2:nip1,:,:),1,2);
    diffv = 0.5*diff(scubep1(:,2:nip1,:),1,1);
    sonu = scube(:,:,:) + masku.*diffu;
    sonv = scube(:,:,:) + maskv.*diffv;

    %  Get B values on the U grid points
    diffu = 0.5*diff(bcubep1(2:nip1,:,:),1,2);
    diffv = 0.5*diff(bcubep1(:,2:nip1,:),1,1);
    bonu = bcube(:,:,:) + masku.*diffu;
    bonv = bcube(:,:,:) + maskv.*diffv;
    
    if ((newyear ~= oldyear) && (nm > 0)) || (it == length(ilist))

      %  Write if its a new year or the last of all entries
      disp(['  ==>  Writing files for ==> ' sprintf('%d',oldyear)]);
      for io = 1:size(vnam,1)
        comm = sprintf( ...
            '%s_ido = fopen(''%s_tave_%d'', ''a'', ''ieee-be'');', ...
            vnam(io,:), vnam(io,:), oldyear );
        eval(comm)
        comm = sprintf( '%s_acc = %s_acc ./ nm;',vnam(io,:), ...
                        vnam(io,:));
        eval(comm)
        comm = sprintf('fwrite(%s_ido, %s_acc, ''real*4'');', ...
                       vnam(io,:), vnam(io,:) );
        eval(comm)
        comm = sprintf( 'fclose(%s_ido);', vnam(io,:) );
        eval(comm)
        
        %  Re-zero the accumulators
        comm = sprintf('%s_acc = zeros(ne,ne,6);', vnam(io,:));
        eval(comm)
      end

      nm = 0;
      
    end
    
    %  Accumulate partial sums
    U__acc = U__acc + ucube;
    V__acc = V__acc + vcube;
    T__acc = T__acc + tcube;
    S__acc = S__acc + scube;
    B__acc = B__acc + bcube;
    B1_acc = B1_acc + bcubem1;
    
    UT_acc = UT_acc + tonu.*ucube;
    VT_acc = VT_acc + tonv.*vcube;
    US_acc = US_acc + sonu.*ucube;
    VS_acc = VS_acc + sonv.*vcube;
    UB_acc = UB_acc + bonu.*ucube;
    VB_acc = VB_acc + bonv.*vcube;

    nm = nm + 1;
    
    oldyear = newyear;
    
  end
end

1	edhill	1.1	function m = calcUVBave_ttcp(Spat, Tpat, Upat, Vpat, ilist, nlay, ne)
2
3			% Function m = calcTBave(Spat, Tpat, Upat, Vpat, ilist, nlay, ne)
4			%
5			%
6			% INPUTS
7			% Tpat string containing T-data file pattern
8			% Upat string containing U-data file pattern
9			% Vpat string containing V-data file pattern
10			% ilist list of iteration values
11			% nlay number of z layers
12			% ne number of values alone each edge of each face
13			%
14			% OUTPUTS
15			% m output array of dimension ns*nps
16
17			if nargin ~= 7
18			disp('Error: wrong number of arguments')
19			exit
20			end
21
22			deltatT = 1200; % model time step size (s)
23			tavefreq = 259200; % averaging period (s)
24			startDate = datenum(1992,1,1); % model integration starting date
25
26			oldyear = -1;
27			newyear = oldyear;
28
29			nps = 6 * ne*ne;
30			tx = 85;
31			ty = 85;
32			nt = nene6/(tx*ty);
33			cx = ne;
34			cy = ne;
35			nz = 1;
36
37			delR = [
38			10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ...
39			10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ...
40			31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ...
41			93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ...
42			139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ...
43			341.50,364.50,387.50,410.50,433.50,456.50 ];
44			R = cumsum(delR) - 0.5*delR;
45
46			vnam = [ 'U_'; 'V_'; 'T_'; 'S_'; 'B_'; 'B1'; ...
47			'UT'; 'VT'; 'US'; 'VS'; 'UB'; 'VB' ];
48
49			lcstr = [ '! ( rm -f %s ; ttcp -r > %s ) > /dev/null 2>&1 &' ];
50			rcstr = [ '!ssh 10.12.2.200 ''( dd bs=6242400 skip=%d ' ...
51			' count=1 if=/home/edhill/%s' ...
52			' \| ttcp -t 10.12.2.41 ) > /dev/null 2>&1''' ];
53
54			disp(['lcstr = "' lcstr '"']);
55			disp(['rcstr = "' rcstr '"']);
56
57			il = 1;
58			for il = 1:nlay
59
60			a = sprintf('Slice: %g',il);
61			% disp(a);
62
63			% Zero the accumulators
64			nm = 0;
65			for io = 1:size(vnam,1)
66			comm = sprintf('%s_acc = zeros(ne,ne,6);', vnam(io,:));
67			eval(comm)
68			end
69
70			it = 1;
71			for it = 1:length(ilist)
72
73			iter = ilist(it);
74			str = datestr(startDate + (iter*deltatT-tavefreq/2)/60/60/24);
75			dv = datevec(startDate + (iter*deltatT-tavefreq/2)/60/60/24);
76			newyear = dv(1);
77
78			sfname = sprintf(Spat, ilist(it));
79			tfname = sprintf(Tpat, ilist(it));
80			ufname = sprintf(Upat, ilist(it));
81			vfname = sprintf(Vpat, ilist(it));
82			disp([ a ' ' str ' ' tfname ]);
83			skip = il - 1;
84			eval( sprintf(lcstr, sfname, sfname) );
85			eval( sprintf(rcstr, skip, sfname) );
86			eval( sprintf(lcstr, tfname, tfname) );
87			eval( sprintf(rcstr, skip, tfname) );
88			eval( sprintf(lcstr, ufname, ufname) );
89			eval( sprintf(rcstr, skip, ufname) );
90			eval( sprintf(lcstr, vfname, vfname) );
91			eval( sprintf(rcstr, skip, vfname) );
92			sfid = fopen(sfname, 'r', 'ieee-be');
93			tfid = fopen(tfname, 'r', 'ieee-be');
94			ufid = fopen(ufname, 'r', 'ieee-be');
95			vfid = fopen(vfname, 'r', 'ieee-be');
96			offset = (il - 1)nps4;
97
98			% Convert the horrible JPL mdsio layout to six sane faces
99			%fseek(sfid, offset, 'bof');
100			scube = unmangleJPL1( reshape(fread(sfid, nps, 'real4'),txnt,ty), ...
101			ne, tx );
102
103			%fseek(tfid, offset, 'bof');
104			tcube = unmangleJPL1( reshape(fread(tfid, nps, 'real4'),txnt,ty), ...
105			ne, tx );
106
107			%fseek(ufid, offset, 'bof');
108			ucube = unmangleJPL1( reshape(fread(ufid, nps, 'real4'),txnt,ty), ...
109			ne, tx );
110
111			%fseek(vfid, offset, 'bof');
112			vcube = unmangleJPL1( reshape(fread(vfid, nps, 'real4'),txnt,ty), ...
113			ne, tx );
114
115			fclose(sfid);
116			fclose(tfid);
117			fclose(ufid);
118			fclose(vfid);
119
120			% Compute the density
121			bcube = zeros(size(scube));
122			bcubem1 = zeros(size(scube));
123			for icf = 1:6
124			% pressure_for_eos.F :
125			% locPres(i,j) = -rhoConstrC(k)gravity*maskC(i,j,k,bi,bj)
126			% rhoConst = 9.997999999999999E+02 /* Ref density ( kg/m^3 ) */
127			% gravity = 9.810000000000000E+00 /* Grav accel ( m/s^2 ) */
128			p = 0.09998 * 9.81 * R(il);
129			bcube(:,:,icf) = ...
130			eh3densjmd95( squeeze(scube(:,:,icf)), ...
131			squeeze(tcube(:,:,icf)), p ) * -9.81 / 1000;
132			if il > 1
133			p = 0.09998 * 9.81 * R(il-1);
134			bcubem1(:,:,icf) = ...
135			eh3densjmd95( scube(:,:,icf), ...
136			tcube(:,:,icf), p ) * -9.81 / 1000;
137			end
138			end
139			tmask = double(tcube ~= 0.0);
140			scube = scube .* tmask;
141			bcube = bcube .* tmask;
142			bcubem1 = bcubem1 .* tmask;
143
144			% Fill the (ne+1)*(ne+1) grid with T,S,B values
145			ni = ne; nip1 = ni + 1;
146			nj = ne; njp1 = nj + 1;
147			tcubep1 = zeros( [ nip1 njp1 6 ] );
148			scubep1 = zeros( [ nip1 njp1 6 ] );
149			bcubep1 = zeros( [ nip1 njp1 6 ] );
150			for icf = 1:6
151			tcubep1(2:nip1,2:njp1,icf) = tcube(:,:,icf);
152			scubep1(2:nip1,2:njp1,icf) = scube(:,:,icf);
153			bcubep1(2:nip1,2:njp1,icf) = bcube(:,:,icf);
154			end
155
156			% Do the upwind-edge T exchanges
157			tcubep1(1,2:njp1,1) = tcube(ni:-1:1,nj,5); % -
158			tcubep1(2:nip1,1,1) = tcube(1:ni,nj,6); % -
159			tcubep1(1,2:njp1,2) = tcube(ni,1:nj,1); % -
160			tcubep1(2:nip1,1,2) = tcube(ni,nj:-1:1,6); % -
161			tcubep1(1,2:njp1,3) = tcube(ni:-1:1,nj,1); % -
162			tcubep1(2:nip1,1,3) = tcube(1:ni,nj,2); % -
163			tcubep1(1,2:njp1,4) = tcube(ni,1:nj,3); % -
164			tcubep1(2:nip1,1,4) = tcube(ni,nj:-1:1,2); % -
165			tcubep1(1,2:njp1,5) = tcube(ni:-1:1,nj,3); % -
166			tcubep1(2:nip1,1,5) = tcube(1:ni,nj,4); % -
167			tcubep1(1,2:njp1,6) = tcube(ni,1:nj,5); % -
168			tcubep1(2:nip1,1,6) = tcube(ni,nj:-1:1,4); % -
169
170			% Do the upwind-edge S exchanges
171			scubep1(1,2:njp1,1) = scube(ni:-1:1,nj,5); % -
172			scubep1(2:nip1,1,1) = scube(1:ni,nj,6); % -
173			scubep1(1,2:njp1,2) = scube(ni,1:nj,1); % -
174			scubep1(2:nip1,1,2) = scube(ni,nj:-1:1,6); % -
175			scubep1(1,2:njp1,3) = scube(ni:-1:1,nj,1); % -
176			scubep1(2:nip1,1,3) = scube(1:ni,nj,2); % -
177			scubep1(1,2:njp1,4) = scube(ni,1:nj,3); % -
178			scubep1(2:nip1,1,4) = scube(ni,nj:-1:1,2); % -
179			scubep1(1,2:njp1,5) = scube(ni:-1:1,nj,3); % -
180			scubep1(2:nip1,1,5) = scube(1:ni,nj,4); % -
181			scubep1(1,2:njp1,6) = scube(ni,1:nj,5); % -
182			scubep1(2:nip1,1,6) = scube(ni,nj:-1:1,4); % -
183
184			% Do the upwind-edge B exchanges
185			bcubep1(1,2:njp1,1) = bcube(ni:-1:1,nj,5); % -
186			bcubep1(2:nip1,1,1) = bcube(1:ni,nj,6); % -
187			bcubep1(1,2:njp1,2) = bcube(ni,1:nj,1); % -
188			bcubep1(2:nip1,1,2) = bcube(ni,nj:-1:1,6); % -
189			bcubep1(1,2:njp1,3) = bcube(ni:-1:1,nj,1); % -
190			bcubep1(2:nip1,1,3) = bcube(1:ni,nj,2); % -
191			bcubep1(1,2:njp1,4) = bcube(ni,1:nj,3); % -
192			bcubep1(2:nip1,1,4) = bcube(ni,nj:-1:1,2); % -
193			bcubep1(1,2:njp1,5) = bcube(ni:-1:1,nj,3); % -
194			bcubep1(2:nip1,1,5) = bcube(1:ni,nj,4); % -
195			bcubep1(1,2:njp1,6) = bcube(ni,1:nj,5); % -
196			bcubep1(2:nip1,1,6) = bcube(ni,nj:-1:1,4); % -
197
198			% Get T values on the U,V grid points
199			masku = 1.0 - abs(diff(double(tcubep1(2:nip1,:,:) == 0.0),1,2));
200			maskv = 1.0 - abs(diff(double(tcubep1(:,2:nip1,:) == 0.0),1,1));
201			diffu = 0.5*diff(tcubep1(2:nip1,:,:),1,2);
202			diffv = 0.5*diff(tcubep1(:,2:nip1,:),1,1);
203			tonu = tcube(:,:,:) + masku.*diffu;
204			tonv = tcube(:,:,:) + maskv.*diffv;
205
206			% Get S values on the U,V grid points
207			diffu = 0.5*diff(scubep1(2:nip1,:,:),1,2);
208			diffv = 0.5*diff(scubep1(:,2:nip1,:),1,1);
209			sonu = scube(:,:,:) + masku.*diffu;
210			sonv = scube(:,:,:) + maskv.*diffv;
211
212			% Get B values on the U grid points
213			diffu = 0.5*diff(bcubep1(2:nip1,:,:),1,2);
214			diffv = 0.5*diff(bcubep1(:,2:nip1,:),1,1);
215			bonu = bcube(:,:,:) + masku.*diffu;
216			bonv = bcube(:,:,:) + maskv.*diffv;
217
218			if ((newyear ~= oldyear) && (nm > 0)) \|\| (it == length(ilist))
219
220			% Write if its a new year or the last of all entries
221			disp([' ==> Writing files for ==> ' sprintf('%d',oldyear)]);
222			for io = 1:size(vnam,1)
223			comm = sprintf( ...
224			'%s_ido = fopen(''%s_tave_%d'', ''a'', ''ieee-be'');', ...
225			vnam(io,:), vnam(io,:), oldyear );
226			eval(comm)
227			comm = sprintf( '%s_acc = %s_acc ./ nm;',vnam(io,:), ...
228			vnam(io,:));
229			eval(comm)
230			comm = sprintf('fwrite(%s_ido, %s_acc, ''real*4'');', ...
231			vnam(io,:), vnam(io,:) );
232			eval(comm)
233			comm = sprintf( 'fclose(%s_ido);', vnam(io,:) );
234			eval(comm)
235
236			% Re-zero the accumulators
237			comm = sprintf('%s_acc = zeros(ne,ne,6);', vnam(io,:));
238			eval(comm)
239			end
240
241			nm = 0;
242
243			end
244
245			% Accumulate partial sums
246			U__acc = U__acc + ucube;
247			V__acc = V__acc + vcube;
248			T__acc = T__acc + tcube;
249			S__acc = S__acc + scube;
250			B__acc = B__acc + bcube;
251			B1_acc = B1_acc + bcubem1;
252
253			UT_acc = UT_acc + tonu.*ucube;
254			VT_acc = VT_acc + tonv.*vcube;
255			US_acc = US_acc + sonu.*ucube;
256			VS_acc = VS_acc + sonv.*vcube;
257			UB_acc = UB_acc + bonu.*ucube;
258			VB_acc = VB_acc + bonv.*vcube;
259
260			nm = nm + 1;
261
262			oldyear = newyear;
263
264			end
265			end
266