function m = calcUVBave(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' ];

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 ]);
    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