%=======================================================
%
%  $Id: plot_all_10.m,v 1.3 2004/08/20 20:36:49 edhill Exp $
%
%  Ed Hill
%

%  The following are the MatLAB commands used to create the various
%  plots related to eddy fluxes using average velocities and densities
%  (called bouyancy or "b" in many of the variables) from Dimitris'
%  "cube_5" integration.

%  Groups of commands contained within the following "UNUSED" comments
%  we not used in this analyssis.  They are "left over" from the
%  previous temperature-based calculations and have been kept for
%  reference purposes only.

%-------  UNUSED  -----------------------
%  ...Commands...
%-------  UNUSED  -----------------------


matlab -nojvm
matlab -nojvm -nodisplay

  %=======================================================
  %  Compute 10-yr averages from the ten 1-yr averages

  clear all ; close all
  ne = 510;
  ntot = ne*ne*6 * 50;
  by = 1994;
  ey = 2003;
  vnam = [ 'U_'; 'V_'; 'T_'; 'S_'; 'B_'; 'B1'; ...
           'UT'; 'VT'; 'US'; 'VS'; 'UB'; 'VB' ];
  nm = 1;
  for nm = 1:size(vnam,1)
    sum = zeros(ntot,1);
    eval(['pat = ''' sprintf('%s',vnam(nm,:)) '_tave_%4.4d'';']);
    disp(sprintf('pat = "%s"',pat));
    iy = by;
    for iy = by:ey
      disp(['  opening:  ' sprintf(pat,iy)]);
      fid = fopen(sprintf(pat,iy), 'r', 'ieee-be');
      a = fread(fid, ntot, 'real*4');
      fclose(fid);
      sum = sum + a;
    end
    a = sum ./ (ey-by+1);
    ido = fopen(sprintf([pat '-%4.4d'],by,ey), 'w', 'ieee-be');
    fwrite(ido, a, 'real*4');
    fclose(ido);
  end
  
  %=======================================================
  %  Compute the "primes"
  
  clear all;  close all
  
  vnam = [ 'U_'; 'V_'; 'T_'; 'S_'; 'B_'; 'B1'; ...
           'UT'; 'VT'; 'US'; 'VS'; 'UB'; 'VB' ];
  nztot = 50;
  ne = 510;
  nps = 6 * ne * ne;
  iz = 2;
  iz = 1;
  for iz = 1:nztot
    
    disp(sprintf('iz = %d', iz));
    
    offset = (iz - 1)*nps*4;
    nm = 1;
    for nm = 1:size(vnam,1)
      eval(['pat = ''' sprintf('%s',vnam(nm,:)) ...
            '_tave_1994-2003'';']);
      fid = fopen(pat, 'r', 'ieee-be');
      fseek(fid, offset, 'bof');
      a = reshape(fread(fid, nps, 'real*4'),ne,ne,6);
      eval([lower(vnam(nm,:)) ' = a;' ]);
      fclose(fid);
    end
    clear a
    %  surf(u_(:,:,1)), view(2), shading interp, axis equal
    %  surf(v_(:,:,1)), view(2), shading interp, axis equal
    %  surf(t_(:,:,1)), view(2), shading interp, axis equal
    %  s_(find(s_==0.0)) = NaN;
    %  surf(s_(:,:,1)), view(2), shading interp, axis equal
    %  b_(find(b_==0.0)) = NaN;
    %  surf(b_(:,:,1)), view(2), shading interp, axis equal
    %  b1(find(b1==0.0)) = NaN;
    %  surf(b1(:,:,1)), view(2), shading interp, axis equal
    %  surf(ut(:,:,1)), view(2), shading interp, axis equal
    %  surf(vt(:,:,1)), view(2), shading interp, axis equal
    %  surf(us(:,:,1)), view(2), shading interp, axis equal
    %  surf(vs(:,:,1)), view(2), shading interp, axis equal
    %  surf(ub(:,:,1)), view(2), shading interp, axis equal
    %  surf(vb(:,:,1)), view(2), shading interp, axis equal

    %  tmask   = double(t_ ~= 0.0);
    
    %  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) = t_(:,:,icf);
      scubep1(2:nip1,2:njp1,icf) = s_(:,:,icf);
      bcubep1(2:nip1,2:njp1,icf) = b_(:,:,icf);
    end

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

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

    %  Do the upwind-edge B exchanges
    bcubep1(1,2:njp1,1) = b_(ni:-1:1,nj,5);   % - 
    bcubep1(2:nip1,1,1) = b_(1:ni,nj,6);      % -
    bcubep1(1,2:njp1,2) = b_(ni,1:nj,1);      % -
    bcubep1(2:nip1,1,2) = b_(ni,nj:-1:1,6);   % -
    bcubep1(1,2:njp1,3) = b_(ni:-1:1,nj,1);   % -
    bcubep1(2:nip1,1,3) = b_(1:ni,nj,2);      % -
    bcubep1(1,2:njp1,4) = b_(ni,1:nj,3);      % -
    bcubep1(2:nip1,1,4) = b_(ni,nj:-1:1,2);   % -
    bcubep1(1,2:njp1,5) = b_(ni:-1:1,nj,3);   % -
    bcubep1(2:nip1,1,5) = b_(1:ni,nj,4);      % -
    bcubep1(1,2:njp1,6) = b_(ni,1:nj,5);      % -
    bcubep1(2:nip1,1,6) = b_(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 = t_(:,:,:) + masku.*diffu;
    tonv = t_(:,:,:) + 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 = s_(:,:,:) + masku.*diffu;
    sonv = s_(:,:,:) + 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 = b_(:,:,:) + masku.*diffu;
    bonv = b_(:,:,:) + maskv.*diffv;
    
    %  Compute the primes
    uptp = ut - u_ .* tonu;
    vptp = vt - v_ .* tonv;
    upsp = us - u_ .* sonu;
    vpsp = vs - v_ .* sonv;
    upbp = ub - u_ .* bonu;
    vpbp = vb - v_ .* bonv;
    
    %  Write the results
    for f = [ 't' 's' 'b' ]
      for uv = [ 'u' 'v' ]
        %  disp(['   writing:  ' uv f]);
        ido = fopen(['out/' uv 'p' f 'p_1994-2003'], 'a', 'ieee-be');
        eval(['fwrite(ido, ' uv 'p' f 'p, ''real*4'');']);
        fclose(ido);
      end
    end
    %     f = [ 't_'; 's_'; 'b_'; 'b1'; 'u_'; 'v_'; ];
    %     for fi = 1:size(f,1)
    %       disp(['   writing:  ' f(fi,:) ]);
    %       ido = fopen(['out/' f(fi,:) '_1994-2003'], 'a', 'ieee-be');
    %       eval(['fwrite(ido, ' f(fi,:) ', ''real*4'');']);
    %       fclose(ido);
    %     end
    
  end
  
  %=======================================================
  %  Plot u'B' and v'B' on a Lat-Lon grid

  clear all, close all
  tx = 85;
  ty = 85;
  nt = 216;
  cx = 510;
  cy = 510;
  nz = 1;
  ne = 510;
  nps = ne * ne * 6;

  %  XCS YCS XGS YGS
  fnam = [ 'XC' ; 'YC'; 'XG'; 'YG' ];
  for in = 1:4
    uid = fopen([fnam(in,:) '.data'], 'r', 'ieee-be');
    phi = unmangleJPL1( reshape(fread(uid, nps, 'real*4'), ...
                                tx*nt,ty), ne, tx );
    fclose(uid);
    eval([lower(fnam(in,:)) ' = phi;']);
    a = zeros(6*510,510);
    for i = 1:6
      xb = (i-1)*510 + 1;  xe = xb + 510 - 1;
      yb = 1;              ye = 510;
      a(xb:xe,yb:ye) = phi(:,:,i);
    end
    eval([lower(fnam(in,:)) 's = a;'])
  end
  xcs = xcs + -360.0*(xcs > 180.0);
  xgs = xgs + -360.0*(xgs > 180.0);
  clear phi a

  %=======================================================
  %  Calculation rotations for Lat-Lon vector alignment
  
  ne = 510;
  n1 = ne - 1;
  dux = zeros(size(xg));
  duy = zeros(size(xg));
  dvx = zeros(size(xg));
  dvy = zeros(size(xg));
  dux(1:n1,:,:) = diff(xg,1,1);  dux(ne,:,:) = dux(n1,:,:);
  dvx(:,1:n1,:) = diff(xg,1,2);  dvx(:,ne,:) = dvx(:,n1,:);
  duy(1:n1,:,:) = diff(yg,1,1);  duy(ne,:,:) = duy(n1,:,:);
  dvy(:,1:n1,:) = diff(yg,1,2);  dvy(:,ne,:) = dvy(:,n1,:);
  dux = dux + 360*double(dux < 180); 
  dux = dux - 360*double(dux > 180);  %  [ min(min(dux)) max(max(dux)) ]
  duy = duy + 360*double(duy < 180);
  duy = duy - 360*double(duy > 180);  %  [ min(min(duy)) max(max(duy)) ]
  dvx = dvx + 360*double(dvx < 180);
  dvx = dvx - 360*double(dvx > 180);  %  [ min(min(dvx)) max(max(dvx)) ]
  dvy = dvy + 360*double(dvy < 180);
  dvy = dvy - 360*double(dvy > 180);  %  [ min(min(dvy)) max(max(dvy)) ]
  llux = dux ./ sqrt(dux.^2 + duy.^2);
  lluy = duy ./ sqrt(dux.^2 + duy.^2);
  llvx = dvx ./ sqrt(dvx.^2 + dvy.^2);
  llvy = dvy ./ sqrt(dvx.^2 + dvy.^2);

  %  surf(xg(:,:,2)), view(2), axis equal, shading interp, colorbar
  %  surf(dux(:,:,2)), view(2), axis equal, shading interp, colorbar
  %  surf(duy(:,:,2)), view(2), axis equal, shading interp, colorbar

  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;
  Rmid = R(1:(length(R)-1)) + 0.5*diff(R);

  ne = 510;
  nps = ne * ne * 6;
        
  clow  = [ -1.0 ];  % -0.2;
  chigh = [  1.0 ];  %  0.2;
  
  iz = 1;
  iz = 15;
  iz = 25;
  for iz = [ 1:50 ]

    disp(sprintf('iz = %d',iz));
    offset = (iz - 1)*nps*4;
    %  Load the "primes"
    for f = [ 't' 's' 'b' ]
      for uv = [ 'u' 'v' ]
        % disp(['   reading:  ' uv f]);
        fid = fopen(['out/' uv 'p' f 'p_1994-2003'],'r','ieee-be');
        fseek(fid, offset, 'bof');
        a = reshape(fread(fid, nps, 'real*4'),ne,ne,6);
        fclose(fid);
        eval([uv 'p' f 'p = a;']);
      end
    end

    %  Get LatLon-oriented vectors at the CubeSphere points
    lluptp = uptp .* llux + vptp .* llvx;
    llvptp = uptp .* lluy + vptp .* llvy;
    llupsp = upsp .* llux + vpsp .* llvx;
    llvpsp = upsp .* lluy + vpsp .* llvy;
    llupbp = upbp .* llux + vpbp .* llvx;
    llvpbp = upbp .* lluy + vpbp .* llvy;
    % surf(xg(:,:,2),yg(:,:,2),lluptp(:,:,2)), view(2),shading interp,colorbar
    % surf(xg(:,:,2),yg(:,:,2),llupbp(:,:,2)), view(2),shading interp,colorbar
    
    %  save the LatLon-rotated "primes"
    for f = [ 't' 's' 'b' ]
      for uv = [ 'u' 'v' ]
        % disp(['   writing:  ' uv f]);
        fid = fopen(['out/ll' uv 'p' f 'p_1994-2003'], ...
                    'a','ieee-be');
        eval(['fwrite(fid,ll' uv 'p' f 'p,''real*4'');']);
        fclose(fid);
      end
    end

    %  Look at Lat-Lon projections of the "primes"
    clim = [ 0.6 0.15 0.0015 ];
    fi = 0;
    for f = [ 't' 's' 'b' ]
      
      fi = fi + 1;
      eval(['llu = up' f 'p;']);
      eval(['llv = vp' f 'p;']);
      
      pllu = zeros(6*510,510);
      pllv = zeros(6*510,510);
      for i = 1:6
        xb = (i-1)*510 + 1;  xe = xb + 510 - 1;
        yb = 1;              ye = 510;
        pllu(xb:xe,yb:ye) = llu(:,:,i);
        pllv(xb:xe,yb:ye) = llv(:,:,i);
      end
      pllu(find(pllu == 0.0)) = NaN;
      pllv(find(pllv == 0.0)) = NaN;
      shift=-1;
      for uv = [ 'u' 'v' ]
        eval(['a = pll' uv ';']);
        clo = -clim(fi);
        chi =  clim(fi);
        grph_CS(a,xcs,ycs,xgs,ygs,[clo],[chi],shift);
        title([ uv '''' upper(f) ''' at ' ...
                sprintf('%g',R(iz)) ...
                'm depth on the 510x510x6 cubesphere ' ...
                'for 1994--2003  ["cube5"]']);
        print('-painters', '-dpng', '-r150', ...
              [ uv 'p' upper(f) 'p_' ...
                sprintf('%02d_%07.1f',iz,R(iz)) 'm_150.png'])
      end
      
    end
  end
  
  %  Save and print the Lat-Lon U,V values
  % !rm -f out/ll[uv]_1994-2003
  iz = 1;
  for iz = [ 1:50 ]

    disp(sprintf('iz = %d',iz));
    offset = (iz - 1)*nps*4;
    %  Load the U,V values
    for uv = [ 'u' 'v' ]
      fid = fopen([upper(uv) '__tave_1994-2003'],'r','ieee-be');
      fseek(fid, offset, 'bof');
      a = reshape(fread(fid, nps, 'real*4'),ne,ne,6);
      fclose(fid);
      eval([uv ' = a;']);
    end

    %  Get LatLon-oriented vectors at the CubeSphere points
    llu = u .* llux + v .* llvx;
    llv = u .* lluy + v .* llvy;
    
    %  save the LatLon-rotated U,V
    for uv = [ 'u' 'v' ]
      % disp(['   writing:  ' uv f]);
      fid = fopen(['out/ll' uv '_1994-2003'],'a','ieee-be');
      eval(['fwrite(fid,ll' uv ',''real*4'');']);
      fclose(fid);
    end

    clim = 1.0;
    pllu = zeros(6*510,510);
    pllv = zeros(6*510,510);
    for i = 1:6
      xb = (i-1)*510 + 1;  xe = xb + 510 - 1;
      yb = 1;              ye = 510;
      pllu(xb:xe,yb:ye) = llu(:,:,i);
      pllv(xb:xe,yb:ye) = llv(:,:,i);
    end
    pllu(find(pllu == 0.0)) = NaN;
    pllv(find(pllv == 0.0)) = NaN;
    shift=-1;
    for uv = [ 'u' 'v' ]
      eval(['a = pll' uv ';']);
      clo = -clim;
      chi =  clim;
      grph_CS(a,xcs,ycs,xgs,ygs,[clo],[chi],shift);
      title([ uv ' at ' ...
              sprintf('%g',R(iz)) ...
              'm depth on the 510x510x6 cubesphere ' ...
              'for 1994--2003  ["cube5"]']);
      print('-painters', '-dpng', '-r150', ...
            [ uv '_' sprintf('%02d_%07.1f',iz,R(iz)) 'm_150.png'])
    end
    
  end

  %=======================================================
  %  Calculate :  (v'B')/(dB/dz)
  
  bid    = fopen('B__tave_1994-2003', 'r', 'ieee-be');
  bm1id  = fopen('B1_tave_1994-2003', 'r', 'ieee-be');
  upbpid = fopen('out/upbp_1994-2003', 'r', 'ieee-be');
  vpbpid = fopen('out/vpbp_1994-2003', 'r', 'ieee-be');
  
  iz = 1;
  offset = (iz - 1)*nps*4;
  fseek(upbpid, offset, 'bof');
  fseek(vpbpid, offset, 'bof');
  upbp = reshape(fread(upbpid, nps, 'real*4'),ne,ne,6);
  vpbp = reshape(fread(vpbpid, nps, 'real*4'),ne,ne,6);
  llvpbp = upbp .* lluy + vpbp .* llvy;
  llvpbpm1 = llvpbp;
  ne = 510;  tx = 85;  ty = 85;  nt = 216;

  iz = 2;
  iz = 9;
  iz = 10;
  % !rm out/stress_1994-2003 out/vpbpdbdz_1994-2003
  for iz = 2:50,
    izm1 = iz - 1;
    disp(sprintf('iz = %d',iz));

    offset = (iz - 1)*nps*4;
    fseek(upbpid, offset, 'bof');
    fseek(vpbpid, offset, 'bof');
    upbp = reshape(fread(upbpid, nps, 'real*4'),ne,ne,6);
    vpbp = reshape(fread(vpbpid, nps, 'real*4'),ne,ne,6);

    %  v'B' on ll coords
    llvpbp = upbp .* lluy + vpbp .* llvy;
    
    %  B and Bm1
    fseek(bid, offset, 'bof');
    b = reshape(fread(bid, nps, 'real*4'),ne,ne,6);
    fseek(bm1id, offset, 'bof');
    bm1 =reshape(fread(bm1id, nps, 'real*4'),ne,ne,6);

    %  ( v'B' )/( dB/dz )
    dbdz = (bm1 - b)/(R(iz) - R(izm1));
    ind0 = find(dbdz==0.0);
    dbdz(ind0) = 1.0;
    rdbdz = 1.0 / (dbdz);
    rdbdz(ind0) = 0.0;
    vpbpdbdz = 0.5*(llvpbp + llvpbpm1) .* rdbdz;
    
    %  Write the results
    mid = fopen('out/vpbpdbdz_1994-2003', 'a', 'ieee-be');
    fwrite(mid, vpbpdbdz, 'real*4');
    fclose(mid);
    
    %  Plot results
    clo = [ -10.0 ];  % -20;
    chi = [  10.0 ];  %  20;
    ll = zeros(6*510,510);
    for i = 1:6
      xb = (i-1)*510 + 1;  xe = xb + 510 - 1;
      yb = 1;              ye = 510;
      ll(xb:xe,yb:ye) = vpbpdbdz(:,:,i);
    end
    shift=-1;
    grph_CS(sq(ll),xcs,ycs,xgs,ygs,[clo],[chi],shift);
    title([ '(v''B'')/(dB/dz) at ' ...
            sprintf('%g',Rmid(iz)) ...
            'm depth on the 510x510x6 cubesphere for 1994--2003  ["cube5"]']);
    % print('-painters', '-dpng', '-r650', ...
    %       ['vpTpdTdz_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_650.png'])
    print('-painters', '-dpng', '-r150', ...
          ['vpBpdBdz_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_150.png'])

    %  Stress: calc, write, and plot
    stress = 1000 * (2*2*pi/(24*3600)*sin(pi*yc/180));
    stress = stress .* vpbpdbdz;
    sid = fopen('out/stress_1994-2003', 'a', 'ieee-be');
    fwrite(sid, stress, 'real*4');
    fclose(sid);
    clim = 1.0;  % clim = 0.2
    clo = [ -clim ];  % [ -1 ];
    chi = [  clim ];  % [  1 ];
    ll = zeros(6*510,510);
    for i = 1:6
      xb = (i-1)*510 + 1;  xe = xb + 510 - 1;
      yb = 1;              ye = 510;
      ll(xb:xe,yb:ye) = stress(:,:,i);
    end
    shift=-1;
    grph_CS(sq(ll),xcs,ycs,xgs,ygs,[clo],[chi],shift);
    title([ 'Stress at ' ...
            sprintf('%g',Rmid(iz)) ...
            'm depth on the 510x510x6 cubesphere for 1994--2003  ["cube5"]']);
    % print('-painters', '-dpng', '-r650', ...
    %       ['stress_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_650.png'])
    print('-painters', '-dpng', '-r150', ...
          ['stress_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_150_clim_0.2.png'])


    %  Next level
    llvpbpm1 = llvpbp;

  end
  
  
  %=======================================================
  %  Zonally average vpbpdbdz and stress

  clear all ; close all
  tx = 85;
  ty = 85;
  nt = 216;
  cx = 510;
  cy = 510;
  nz = 1;
  ne = 510;
  nps = ne * ne * 6;

  %  XCS YCS XGS YGS
  fnam = [ 'XC' ; 'YC'; 'XG'; 'YG' ];
  for in = 1:4
    uid = fopen([fnam(in,:) '.data'], 'r', 'ieee-be');
    phi = unmangleJPL1( reshape(fread(uid, nps, 'real*4'), ...
                                tx*nt,ty), ne, tx );
    fclose(uid);
    eval([lower(fnam(in,:)) ' = phi;']);
    a = zeros(6*510,510);
    for i = 1:6
      xb = (i-1)*510 + 1;  xe = xb + 510 - 1;
      yb = 1;              ye = 510;
      a(xb:xe,yb:ye) = phi(:,:,i);
    end
    eval([lower(fnam(in,:)) 's = a;'])
  end
  xcs = xcs + -360.0*(xcs > 180.0);
  xgs = xgs + -360.0*(xgs > 180.0);
  clear phi a

  nz = 1;
  nr = 50;
  nrm1 = nr - 1;
  nlat = 181;  nlatm1 = nlat - 1;

  hvals = linspace(-90,90,nlat);
  %  save indicies for zonal averages
  i = 2;
  for i = 2:nlat
    inds = find(hvals(i-1)<yg & yg<hvals(i));
    comm = sprintf('inds%04d = uint32(inds);',i-1);
    eval(comm);
  end

  %  Zonally average Lat-Lon [uv]p[tsb]p
  for f = [ 't' 's' 'b' ]
    for uv = [ 'u' 'v' ]
      acc = zeros(nlatm1, nrm1);
      num = zeros(size(acc));
      ne = 510;
      nps = ne * ne * 6;
      fname = [ 'out/ll' uv 'p' f 'p_1994-2003' ];
      disp(['fname = "' fname '"']);
      zid = fopen(fname, 'r', 'ieee-be');
      iz = 1;
      for iz = 1:50,
        disp(sprintf('   iz = %d',iz));
        offset = (iz - 1)*nps*4;
        fseek(zid, offset, 'bof');
        var = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
        
        for jj = 1:nlatm1
          eval( sprintf('clear inds; inds = inds%04d;',jj) );
          tmp = var(inds);
          nzinds = find(tmp ~= 0.0);
          num(jj,iz) = length(nzinds);
          acc(jj,iz) = sum(tmp(nzinds));
        end
      end
      fclose(zid);
      eval(['z' uv 'p' f 'p = acc ./ num;']);
      eval(['save z' uv 'p' f 'p z' uv 'p' f 'p']);
    end
  end
  
  %  zonally average T,S,B
  for f = [ 't' 's' 'b' ]
    acc = zeros(nlatm1, 50);
    num = zeros(size(acc));
    ne = 510;
    nps = ne * ne * 6;
    fname = [upper(f) '__tave_1994-2003'];
    disp(fname);
    zid = fopen(fname, 'r', 'ieee-be' );
    iz = 1;
    for iz = 1:50,
      disp(sprintf('   iz = %d',iz));
      offset = (iz - 1)*nps*4;
      fseek(zid, offset, 'bof');
      var = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
      jj = 22;
      for jj = 1:nlatm1
        eval( sprintf('clear inds; inds = inds%04d;',jj) );
        tmp = var(inds);
        nzinds = find(tmp ~= 0.0);
        num(jj,iz) = length(nzinds);
        acc(jj,iz) = sum(tmp(nzinds));
      end
    end
    fclose(zid);
    eval(['zonal' upper(f) ' = acc ./ num;']);
    eval(['save zonal' upper(f) ' zonal' upper(f)]);
  end
  
  %  zonally average Lat-Lon u,v
  for f = [ 'u' 'v' ]
    acc = zeros(nlatm1, 50);
    num = zeros(size(acc));
    ne = 510;
    nps = ne * ne * 6;
    fname = ['out/ll' f '_1994-2003'];
    disp(fname);
    zid = fopen(fname, 'r', 'ieee-be' );
    iz = 1;
    for iz = 1:50,
      disp(sprintf('   iz = %d',iz));
      offset = (iz - 1)*nps*4;
      fseek(zid, offset, 'bof');
      var = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
      jj = 22;
      for jj = 1:nlatm1
        eval( sprintf('clear inds; inds = inds%04d;',jj) );
        tmp = var(inds);
        nzinds = find(tmp ~= 0.0);
        num(jj,iz) = length(nzinds);
        acc(jj,iz) = sum(tmp(nzinds));
      end
    end
    fclose(zid);
    eval(['zonal' upper(f) ' = acc ./ num;']);
    eval(['save zonal' upper(f) ' zonal' upper(f)]);
  end
  
  %  zonally average vpbpdbdz
  acc = zeros(nlatm1, nrm1);
  num = zeros(size(acc));
  ne = 510;
  nps = ne * ne * 6;
  zid = fopen('out/vpbpdbdz_1994-2003', 'r', 'ieee-be');
  iz = 1;
  for iz = 1:49,
    disp(sprintf('iz = %d',iz));
    offset = (iz - 1)*nps*4;
    fseek(zid, offset, 'bof');
    vpbpdbdz = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
    
    for jj = 1:nlatm1
      eval( sprintf('clear inds; inds = inds%04d;',jj) );
      tmp = vpbpdbdz(inds);
      nzinds = find(tmp ~= 0.0);
      num(jj,iz) = length(nzinds);
      acc(jj,iz) = sum(tmp(nzinds));
    end
  end
  fclose(zid);
  zvpbpdbdz = acc ./ num;
  save zvpbpdbdz zvpbpdbdz
  
  %  zonally average stress
  acc = zeros(nlatm1, nrm1);
  num = zeros(size(acc));
  ne = 510;
  nps = ne * ne * 6;
  zid = fopen('out/stress_1994-2003', 'r', 'ieee-be');
  iz = 1;
  for iz = 1:49,
    disp(sprintf('iz = %d',iz));
    offset = (iz - 1)*nps*4;
    fseek(zid, offset, 'bof');
    stress = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
    jj = 22;
    for jj = 1:nlatm1
      eval( sprintf('clear inds; inds = inds%04d;',jj) );
      tmp = stress(inds);
      nzinds = find(tmp ~= 0.0);
      num(jj,iz) = length(nzinds);
      acc(jj,iz) = sum(tmp(nzinds));
    end
  end
  fclose(zid);
  stress = acc ./ num;
  save stress stress
  
  
  %  zonally average stress IN A SPECIFIC LONGITUDE RANGE
  hvals = linspace(-90,90,nlat);
  %  save indicies for zonal averages
  i = 2;
  for i = 2:nlat
    inds = find(hvals(i-1)<yg & yg<hvals(i) & 210.0<xg & xg<230.0);
    comm = sprintf('inds%04d = uint32(inds);',i-1);
    eval(comm);
  end
  acc = zeros(nlatm1, nrm1);
  num = zeros(size(acc));
  ne = 510;
  nps = ne * ne * 6;
  zid = fopen('out/stress_1994-2003', 'r', 'ieee-be');
  iz = 1;
  for iz = 1:49,
    disp(sprintf('iz = %d',iz));
    offset = (iz - 1)*nps*4;
    fseek(zid, offset, 'bof');
    stress = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
    jj = 22;
    for jj = 1:nlatm1
      eval( sprintf('clear inds; inds = inds%04d;',jj) );
      tmp = stress(inds);
      nzinds = find(tmp ~= 0.0);
      num(jj,iz) = length(nzinds);
      acc(jj,iz) = sum(tmp(nzinds));
    end
  end
  fclose(zid);
  s_210_230 = acc ./ num;
  save s_210_230 s_210_230
  
  %-------  UNUSED  -----------------------
  %  Vertical gradient of zonally averaged B (dB/dz)
  load zonalT
  load allR
  nz = size(zonalT,1);
  dzaTdz = zeros(nz,size(zonalT,2)-1);
  for iz = 2:50
    dzaTdz(:,iz-1) = (zonalT(:,iz) - zonalT(:,iz-1)) ./ (R(iz) - R(iz-1));
  end
  %  save dzaTdz dzaTdz
  %-------  UNUSED  -----------------------

  
  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;
  Rmid = R(1:(length(R)-1)) + 0.5*diff(R);
  hmid = [ -89.5:1:89.5 ];
  save allR delR R Rmid hmid
  
  %-------  UNUSED  -----------------------
  surf(hmid,-Rmid,num'), view(2), shading interp, colorbar
  caxis([ -1 1 ])
  caxis('manual') 
  surf(hmid,-Rmid,mu'), view(2), shading interp, colorbar
  title('Zonally Averaged Stress for 1994--2003 [cube5]')
  xlabel('Latitude [deg]')
  zlabel('Depth [m]')
  axis([ -90 90 -5500 200 ])
  
  %  print -dps -painters t_001.ps
  print('-painters', '-dpng', '-r150', 'za_str_94--03_r150.png')
  print('-painters', '-dpng', '-r650', 'za_str_94--03_r650.png')
  %-------  UNUSED  -----------------------

  
  %=======================================================
  %  Put it all in NetCDF

  clear all ; close all

  load allR
  
  load zonalU  ; i1 = 'U';   u1 = 'm/s';
  load zonalV  ; i2 = 'V';   u2 = 'm/s';
  load zonalT  ; i3 = 'T';   u3 = 'deg C';
  load zonalS  ; i4 = 'S';   u4 = 'PSU';
  load zonalB  ; i5 = 'B';   u5 = 'm/(s^2)';
  load zuptp   ; i10 = 'uT'; u10 = 'm*(deg C)/s';
  load zvptp   ; i11 = 'vT'; u11 = 'm*(deg C)/s';
  load zupsp   ; i12 = 'uS'; u12 = 'm*(PSU)/s';
  load zvpsp   ; i13 = 'vS'; u13 = 'm*(PSU)/s';
  load zupbp   ; i14 = 'uB'; u14 = 'm^2/(s^3)';
  load zvpbp   ; i15 = 'vB'; u15 = 'm^2/(s^3)';
  load zvpbpdbdz ; i20 = 'zvpbpdbdz'; u20 = 'm^2/s';
  load stress    ; i21 = 'stress';    u21 = 'N/(m^2)';
  load s_210_230 ; i30 = 'stress_210_230';  u30 = 'N/(m^2)';
  
  !rm -f Zave_94-03.nc
  nc = netcdf('Zave_94-03.nc', 'clobber');
  nc.reference = 'Zonal averages from Dimitris "cube 5" integrations';
  nc.author = 'Ed Hill <eh3@mit.edu>';
  nc.date = 'Aug 19, 2004';
  nc('X') = length(hmid);
  nc('Y') = length(R);
  nc('Ym1') = length(R) - 1;
  nc{'X'}  = 'X';
  nc{'Y'}  = 'Y';
  nc{'Ym1'}  = 'Ym1';
  for ii = [ 1:5 10:15 ]
    eval(sprintf('nc{ i%d } = { ''Y'', ''X'' };',ii));
  end
  nc{ i20 } = { 'Ym1', 'X' };
  nc{ i21 } = { 'Ym1', 'X' };
  nc{ i30 } = { 'Ym1', 'X' };
  nc{'X'}.uniquename = 'X';
  nc{'X'}.long_name = 'latitude';
  nc{'X'}.gridtype = ncint(0);
  nc{'X'}.units = 'degree_north';
  nc{'Y'}.uniquename = 'Y';
  nc{'Y'}.long_name = 'depth';
  nc{'Y'}.gridtype = ncint(1);
  nc{'Y'}.units = 'm';
  nc{'Ym1'}.uniquename = 'Ym1';
  nc{'Ym1'}.long_name = 'depth';
  nc{'Ym1'}.gridtype = ncint(1);
  nc{'Ym1'}.units = 'm';
  for ii = [ 1:5 10:15 20:21 30 ]
    eval(sprintf('nc{ i%d }.units = u%d;',ii,ii));
    eval(sprintf('nc{ i%d }.long_name = i%d;',ii,ii));
    eval(sprintf('nc{ i%d }.missing_value = ncdouble(NaN);',ii));
    eval(sprintf('nc{ i%d }.FillValue_ = ncdouble(0.);',ii));
  end
  nc{'X'}(:) = hmid;
  nc{'Y'}(:) = -R;
  nc{'Ym1'}(:) = -Rmid;
  nc{ i1 }(:) = permute(zonalU,[ 2 1 ]);
  nc{ i2 }(:) = permute(zonalV,[ 2 1 ]);
  nc{ i3 }(:) = permute(zonalT,[ 2 1 ]);
  nc{ i4 }(:) = permute(zonalS,[ 2 1 ]);
  nc{ i5 }(:) = permute(zonalB,[ 2 1 ]);
  nc{ i10 }(:) = permute(zuptp,[ 2 1 ]);
  nc{ i11 }(:) = permute(zvptp,[ 2 1 ]);
  nc{ i12 }(:) = permute(zupsp,[ 2 1 ]);
  nc{ i13 }(:) = permute(zvpsp,[ 2 1 ]);
  nc{ i14 }(:) = permute(zupbp,[ 2 1 ]);
  nc{ i15 }(:) = permute(zvpbp,[ 2 1 ]);
  nc{ i20 }(:) = permute(zvpbpdbdz,[ 2 1 ]);
  nc{ i21 }(:) = permute(stress,[ 2 1 ]);
  nc{ i30 }(:) = permute(s_210_230,[ 2 1 ]);
  nc = close(nc);

  %  !scp Zave_94-03.nc ingrid.mit.edu:INGRID_PEOPLE/EH3/eddy_flux/cube_5/

  
  %=======================================================
  %  Create a new NetCDF file containing the 
  %  temperature-averaged fields:  stress, U, v'b'

  
  clear all ; close all
  tx = 85;
  ty = 85;
  nt = 216;
  cx = 510;
  cy = 510;
  nz = 1;
  ne = 510;
  nps = ne * ne * 6;

  %  XCS YCS XGS YGS
  fnam = [ 'XC' ; 'YC'; 'XG'; 'YG' ];
  for in = 1:4
    uid = fopen([fnam(in,:) '.data'], 'r', 'ieee-be');
    phi = unmangleJPL1( reshape(fread(uid, nps, 'real*4'), ...
                                tx*nt,ty), ne, tx );
    fclose(uid);
    eval([lower(fnam(in,:)) ' = phi;']);
    a = zeros(6*510,510);
    for i = 1:6
      xb = (i-1)*510 + 1;  xe = xb + 510 - 1;
      yb = 1;              ye = 510;
      a(xb:xe,yb:ye) = phi(:,:,i);
    end
    eval([lower(fnam(in,:)) 's = a;'])
  end
  xcs = xcs + -360.0*(xcs > 180.0);
  xgs = xgs + -360.0*(xgs > 180.0);
  clear phi a

  nr = 50;
  nrm1 = nr - 1;

  %  Get the SST
  f = [ 't' ];
  ne = 510;
  nps = ne * ne * 6;
  fname = [upper(f) '__tave_1994-2003'];
  zid = fopen(fname, 'r', 'ieee-be' );
  iz = 1;
  offset = (iz - 1)*nps*4;
  fseek(zid, offset, 'bof');
  sst = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
  fclose(zid);
  sst(find(sst==0.0)) = NaN;
  %  surf(sst(:,:,1)), shading interp, view(2)
  %  [ squeeze(min(min(sst)))' ; squeeze(max(max(sst)))' ]'
  %  [ squeeze(min(min(yg)))' ; squeeze(max(max(yg)))' ]'

  %  Get the SST-based indicies
  ntemp = 51;
  tval = linspace(-3,33,ntemp);
  for i = 2:length(tval)
    %  Northern Hemisphere (yg > -10)
    inds = find(yg>-10.0 & tval(i-1)<sst & sst<tval(i));
    eval(sprintf('ni%04d = uint32(inds);',i-1));
    %  Southern Hemisphere (yg <  10)
    inds = find(yg<10.0  & tval(i-1)<sst & sst<tval(i));
    eval(sprintf('si%04d = uint32(inds);',i-1));
  end
  tmid = tval(1:(length(tval)-1)) + diff(tval);
  save tval tval tmid
  
  %  "temperature-average" the stress field
  nlatm1 = length(tval) - 1;
  accn = zeros(nlatm1, nrm1);
  accs = zeros(nlatm1, nrm1);
  nnum = zeros(size(nacc));
  snum = zeros(size(nacc));
  ne = 510;
  nps = ne * ne * 6;
  zid = fopen('out/stress_1994-2003', 'r', 'ieee-be');
  iz = 1;
  for iz = 1:48,
    disp(sprintf('iz = %d',iz));
    offset = (iz - 1)*nps*4;
    fseek(zid, offset, 'bof');
    stress = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
    jj = 22;
    for jj = 1:nlatm1
      eval( sprintf('clear ni; ni = ni%04d;',jj) );
      eval( sprintf('clear si; si = si%04d;',jj) );
      ntmp = stress(ni);
      stmp = stress(si);
      nz_ni = find(ntmp ~= 0.0);
      nz_si = find(stmp ~= 0.0);
      nnum(jj,iz) = length(nz_ni);
      snum(jj,iz) = length(nz_si);
      nacc(jj,iz) = sum(ntmp(nz_ni));
      sacc(jj,iz) = sum(stmp(nz_si));
    end
  end
  fclose(zid);
  ntstress = nacc ./ nnum;
  ststress = sacc ./ snum;
  save t_stress ntstress ststress

  
  %  "temperature-average" the U field
  nlatm1 = length(tval) - 1;
  accn = zeros(nlatm1, nrm1);
  accs = zeros(nlatm1, nrm1);
  nnum = zeros(size(nacc));
  snum = zeros(size(nacc));
  ne = 510;
  nps = ne * ne * 6;
  zid = fopen('out/llu_1994-2003', 'r', 'ieee-be');
  iz = 1;
  for iz = 1:49,
    disp(sprintf('iz = %d',iz));
    offset = (iz - 1)*nps*4;
    fseek(zid, offset, 'bof');
    stress = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
    jj = 22;
    for jj = 1:nlatm1
      eval( sprintf('clear ni; ni = ni%04d;',jj) );
      eval( sprintf('clear si; si = si%04d;',jj) );
      ntmp = stress(ni);
      stmp = stress(si);
      nz_ni = find(ntmp ~= 0.0);
      nz_si = find(stmp ~= 0.0);
      nnum(jj,iz) = length(nz_ni);
      snum(jj,iz) = length(nz_si);
      nacc(jj,iz) = sum(ntmp(nz_ni));
      sacc(jj,iz) = sum(stmp(nz_si));
    end
  end
  fclose(zid);
  ntu = nacc ./ nnum;
  stu = sacc ./ snum;
  save t_u ntu stu

  
  %  "temperature-average" the llvpbp field
  nlatm1 = length(tval) - 1;
  accn = zeros(nlatm1, nrm1);
  accs = zeros(nlatm1, nrm1);
  nnum = zeros(size(nacc));
  snum = zeros(size(nacc));
  ne = 510;
  nps = ne * ne * 6;
  zid = fopen('out/llvpbp_1994-2003', 'r', 'ieee-be');
  iz = 1;
  for iz = 1:49,
    disp(sprintf('iz = %d',iz));
    offset = (iz - 1)*nps*4;
    fseek(zid, offset, 'bof');
    stress = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
    jj = 22;
    for jj = 1:nlatm1
      eval( sprintf('clear ni; ni = ni%04d;',jj) );
      eval( sprintf('clear si; si = si%04d;',jj) );
      ntmp = stress(ni);
      stmp = stress(si);
      nz_ni = find(ntmp ~= 0.0);
      nz_si = find(stmp ~= 0.0);
      nnum(jj,iz) = length(nz_ni);
      snum(jj,iz) = length(nz_si);
      nacc(jj,iz) = sum(ntmp(nz_ni));
      sacc(jj,iz) = sum(stmp(nz_si));
    end
  end
  fclose(zid);
  ntllvpbp = nacc ./ nnum;
  stllvpbp = sacc ./ snum;
  save t_llvpbp  ntllvpbp stllvpbp

  
  clear all ; close all
  
  load allR
  load tval
  load t_stress;
    i1 = 'nt_stress';    u1 = 'N/(m^2)';
    i2 = 'st_stress';    u2 = 'N/(m^2)';
  load t_llvpbp
    i3 = 'nt_vpbp';      u3 = 'm^2/(s^3)';
    i4 = 'st_vpbp';      u4 = 'm^2/(s^3)';
  load t_u
    i5 = 'nt_u';         u5 = 'deg C';
    i6 = 'st_u';         u6 = 'deg C';

  !rm -f SSTave_94-03.nc
  nc = netcdf('SSTave_94-03.nc', 'clobber');
  nc.reference = 'SST-based averages from Dimitris "cube 5" integrations';
  nc.author = 'Ed Hill <eh3@mit.edu>';
  nc.date = 'Aug 20, 2004';
  nc('X') = length(tmid);
  nc('Y') = length(Rmid);
  nc('Ym1') = length(Rmid) - 1;
  Rmidm1 = Rmid([1:(length(Rmid) - 1)]);
  nc{'X'}  = 'X';
  nc{'Y'}  = 'Y';
  nc{'Ym1'}  = 'Ym1';
  for ii = [ 1:2 ]
    eval(sprintf('nc{ i%d } = { ''Ym1'', ''X'' };',ii));
  end
  for ii = [ 3:6 ]
    eval(sprintf('nc{ i%d } = { ''Y'', ''X'' };',ii));
  end
  nc{'X'}.uniquename = 'X';
  nc{'X'}.long_name = 'SST';
  nc{'X'}.gridtype = ncint(0);
  nc{'X'}.units = 'deg C';
  nc{'Y'}.uniquename = 'Y';
  nc{'Y'}.long_name = 'depth';
  nc{'Y'}.gridtype = ncint(1);
  nc{'Y'}.units = 'm';
  nc{'Ym1'}.uniquename = 'Ym1';
  nc{'Ym1'}.long_name = 'depth';
  nc{'Ym1'}.gridtype = ncint(1);
  nc{'Ym1'}.units = 'm';
  for ii = [ 1:6 ]
    eval(sprintf('nc{ i%d }.units = u%d;',ii,ii));
    eval(sprintf('nc{ i%d }.long_name = i%d;',ii,ii));
    eval(sprintf('nc{ i%d }.missing_value = ncdouble(NaN);',ii));
    eval(sprintf('nc{ i%d }.FillValue_ = ncdouble(0.);',ii));
  end
  nc{'X'}(:) = tmid;
  nc{'Y'}(:) = -Rmid;
  nc{'Ym1'}(:) = -Rmidm1;
  nc{ i1 }(:) = permute(ntstress,[ 2 1 ]);
  nc{ i2 }(:) = permute(ststress,[ 2 1 ]);
  nc{ i3 }(:) = permute(ntllvpbp,[ 2 1 ]);
  nc{ i4 }(:) = permute(stllvpbp,[ 2 1 ]);
  nc{ i5 }(:) = permute(ntu,[ 2 1 ]);
  nc{ i6 }(:) = permute(stu,[ 2 1 ]);
  nc = close(nc);

  %  !scp SSTave_94-03.nc ingrid.mit.edu:INGRID_PEOPLE/EH3/eddy_flux/cube_5/