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

vnam = [ 'T2'; 'S2'; 'B2'; 'U2'; 'V2'; 'UV'; ];

% lcstr = [ '! ( rm -f %s ; ttcp -r > %s ) > /dev/null 2>&1 &' ];
% rcstr = [ '!ssh -x 10.12.2.%d ''( dd bs=6242400 skip=%d ' ...
%           ' count=1 if=/home/edhill/%s' ...
%           ' | ttcp -t 10.12.2.%d ) > /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, ips(2), skip,   sfname, ips(1)) );
    % eval( sprintf(lcstr,         tfname, tfname) );
    % eval( sprintf(rcstr, ips(3), skip,   tfname, ips(1)) );
    % eval( sprintf(lcstr,         ufname, ufname) );
    % eval( sprintf(rcstr, ips(4), skip,   ufname, ips(1)) );
    % eval( sprintf(lcstr,         vfname, vfname) );
    % eval( sprintf(rcstr, ips(5), skip,   vfname, ips(1)) );
    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) U grid and 
    %       the (ne+1)*(ne) V with initial values
    ni = ne;    nip1 = ni + 1;
    nj = ne;    njp1 = nj + 1;
    ucubep1 = ones( [ nip1 nj   6 ] );
    vcubep1 = ones( [ ni   njp1 6 ] );
    for icf = 1:6
      ucubep1(1:ni,1:nj,icf) = ucube(:,:,icf);
      vcubep1(1:ni,1:nj,icf) = vcube(:,:,icf);
    end

    %  Do the "n+1" U exchanges
    ucubep1(nip1,:,1) =  ucube(1,:,2);       % - 
    ucubep1(nip1,:,2) =  vcube(ne:-1:1,1,4); % - 
    ucubep1(nip1,:,3) =  ucube(1,:,4);       % - 
    ucubep1(nip1,:,4) =  vcube(ne:-1:1,1,6); % - 
    ucubep1(nip1,:,5) =  ucube(1,:,6);       % - 
    ucubep1(nip1,:,6) =  vcube(ne:-1:1,1,2); % - 
    
    %  Do the "n+1" V exchanges
    vcubep1(:,nip1,1) =  ucube(1,ne:-1:1,3); % - 
    vcubep1(:,nip1,2) =  vcube(1,:,3);       % - 
    vcubep1(:,nip1,3) =  ucube(1,ne:-1:1,5); % - 
    vcubep1(:,nip1,4) =  vcube(:,1,5);       % - 
    vcubep1(:,nip1,5) =  ucube(1,ne:-1:1,1); % - 
    vcubep1(:,nip1,6) =  vcube(:,1,1);       % - 
    
    %  Get U,V values on the T grid points
    masku = 1.0 - abs(diff(double(ucubep1(:,:,:) == 0.0),1,1));
    maskv = 1.0 - abs(diff(double(vcubep1(:,:,:) == 0.0),1,2));
    diffu = 0.5*diff(ucubep1(:,:,:),1,1);
    diffv = 0.5*diff(vcubep1(:,:,:),1,2);
    uont = ucube(:,:,:) + masku.*diffu;
    vont = vcube(:,:,:) + 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
    T2_acc = T2_acc + tcube.*tcube;
    S2_acc = S2_acc + scube.*scube;
    B2_acc = B2_acc + bcube.*bcube;
    U2_acc = U2_acc + ucube.*ucube;
    V2_acc = V2_acc + vcube.*vcube;
    UV_acc = UV_acc + ucube.*vcube;

    %     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