eddy_flux/c20a/plot_c20a.m

%=======================================================
%
%  $Header:  $
%
%  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_20a" or "c20a" 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  -----------------------


%  ssh eddy
%  cd /r/r0/edhill/eddy_stats/eddy_flux_c20a

%  matlab -nojvm
%  matlab -nojvm -nodisplay

  clear all
  close all


  %=======================================================
  %  Project fields to lower-res Lat-Lon and write as 
  %  NetCDF for viewing with Ingrid
  gvars = { 'XC','YC','DXF','DYF','RA','XG','YG','DXV', ...
            'DYU','RAZ','DXC','DYC','RAW','RAS','DXG','DYG' };

  iface = 1;
  for iface = 1:6
    fname = sprintf('grid/tile%03d.mitgrid', iface);
    gid = fopen(fname, 'r', 'ieee-be');
    tmp = reshape(fread(gid,inf,'real*8',0,'ieee-be'),[511,511,16]);
    fclose(gid);
    %  surf(tmp(:,:,1)), view(2), shading interp
    %  for jj = 1:length(gvars)
    for jj = 1:7
      comm = sprintf('%s(:,:,%d) = tmp(:,:,%d);', ...
                     [gvars{jj}], iface, jj);
      eval(comm);
    end
  end
  %  surf(XC(:,:,1)), view(2), shading interp
  %  subplot(2,1,1),  a = [1:10];     surf(XC(a,a,1)), view(2)
  %  subplot(2,1,2),  a = [500:511];  surf(XC(a,a,1)), view(2)
  %  surf(YC(:,:,1)), view(2), shading interp
  %  surf(XG(:,:,1)), view(2), shading interp
  %  surf(YG(:,:,1)), view(2), shading interp
  i510 = [1:510];
  v510 = { 'XC','YC','DXF','DYF','RA' };
  for i = 1:length(v510)
    eval(sprintf('%s = %s(i510,i510,:);',v510{i},v510{i}));
  end

  %  !echo `ls -1 ave__92_99 | sed -e 's|\.ave||g' \
  %    | sed -e "s|^|'|g" | sed -e "s|$|',|g"`
  fields = { ...
      'DRHODR', 'ETANSQ', 'RHOANOSQ', 'RHOAnoma', 'SALT', ...
      'SALTSQ', 'SFLUX', 'SRELAX', 'TAUX', 'TAUY', 'TFLUX', ...
      'THETA', 'THETASQ', 'TICE', 'TRELAX', 'URHOMASS', ...
      'USLTMASS', 'UTHMASS', 'UVEL', 'UVELMASS', 'UVELSQ', ...
      'UV_VEL_Z', 'VRHOMASS', 'VSLTMASS', 'VTHMASS', ...
      'VVEL', 'VVELMASS', 'VVELSQ', 'WRHOMASS', 'WSLTMASS', ...
      'WTHMASS', 'WU_VEL', 'WVELMASS', 'WVELSQ', 'WV_VEL' };
  

  %=======================================================
  %  Compute UpBp and VpBp

  
%-------  UNUSED  -----------------------

  
  %=======================================================
  %  Compute u'b' and v'b'
  
  clear all;  close all
  
  bfid  = fopen( 'b-ave-1992-2003', 'r', 'ieee-be');
  ufid  = fopen( 'U_ave_1992-2003', 'r', 'ieee-be');
  vfid  = fopen( 'V_ave_1992-2003', 'r', 'ieee-be');
  ubfid = fopen('ub-ave-1992-2003', 'r', 'ieee-be');
  vbfid = fopen('vb-ave-1992-2003', 'r', 'ieee-be');
  
  afid = [ ufid  vfid ];
  av   = [ 'ua'; 'va' ];
  
  nztot = 50;
  ne = 510;
  nps = 6 * ne * ne;
  iz = 1;
  for iz = 1:nztot
    
    disp(sprintf('iz = %d', iz));

    offset = (iz - 1)*nps*4;
    id = 1;
    for id = 1:length(afid)
      fseek(afid(id), offset, 'bof');
      %  Convert the horrible JPL mdsio layout to six saner faces
      ne = 510;  tx = 85;  ty = 85;  nt = 216;
      phi = unmangleJPL1( reshape(fread(afid(id), nps, 'real*4'), ...
                                  tx*nt,ty), ne, tx );
      comm = [ av(id,:) ' = phi;' ];
      %  disp(comm);
      eval(comm);
    end
    clear phi
    %  surf(ua(:,:,1)), view(2), shading interp, axis equal
    %  surf(va(:,:,1)), view(2), shading interp, axis equal

    fseek(bfid, offset, 'bof');
    b = reshape(fread(bfid, nps, 'real*4'),ne,ne,6);
    %  surf(b(:,:,1)), view(2), shading interp, axis equal
    fseek(ubfid, offset, 'bof');
    ub = reshape(fread(ubfid, nps, 'real*4'),ne,ne,6);
    %  surf(ub(:,:,1)), view(2), shading interp, axis equal
    fseek(vbfid, offset, 'bof');
    vb = reshape(fread(vbfid, nps, 'real*4'),ne,ne,6);
    %  surf(vb(:,:,1)), view(2), shading interp, axis equal
    
    %  Fill the (ne+1)*(ne+1) grid with B values
    ni = ne;    nip1 = ni + 1;
    nj = ne;    njp1 = nj + 1;
    bcubep1 = zeros( [ nip1 njp1 6 ] );
    for icf = 1:6
      bcubep1(2:nip1,2:njp1,icf) = b(:,:,icf);
    end

    %  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 B values on the U grid points
    masku = 1.0 - abs(diff(double(bcubep1(2:nip1,:,:) == 0.0),1,2));
    diffu = 0.5*diff(bcubep1(2:nip1,:,:),1,2);
    bonu = b(:,:,:) + masku.*diffu;
    
    %  Get B values on the V grid points
    maskv = 1.0 - abs(diff(double(bcubep1(:,2:nip1,:) == 0.0),1,1));
    diffv = 0.5*diff(bcubep1(:,2:nip1,:),1,1);
    bonv = b(:,:,:) + maskv.*diffv;
    
    %  Compute u'B' and v'B'
    upbp = ub - ua .* bonu;
    vpbp = vb - va .* bonv;
    
    %  Write the results
    ubid = fopen('upbp-1992-2003', 'a', 'ieee-be');
    vbid = fopen('vpbp-1992-2003', 'a', 'ieee-be');
    fwrite(ubid, upbp, 'real*4');
    fwrite(vbid, vpbp, 'real*4');
    fclose(ubid);
    fclose(vbid);
    
  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;
    
  %  u'T' and v'T'
  upbpid = fopen('upbp-1992-2003', 'r', 'ieee-be');
  vpbpid = fopen('vpbp-1992-2003', 'r', 'ieee-be');

  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));
    %  iz = 21;
    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);

  
    %  Look at u'b' and v'b' on ll coords
    % surf(xg(:,:,2), yg(:,:,2), upbp(:,:,2)), view(2), shading interp, colorbar
    % surf(xg(:,:,3), yg(:,:,3), upbp(:,:,3)), view(2), shading interp, colorbar
    % surf(xg(:,:,4), yg(:,:,4), upbp(:,:,4)), view(2), shading interp, colorbar
    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), llvptp(:,:,2)), view(2), shading
    % interp, colorbar
    
    %  save for later use
    upbpllid = fopen('upbpll-1992-2003', 'a', 'ieee-be');
    vpbpllid = fopen('vpbpll-1992-2003', 'a', 'ieee-be');
    fwrite(upbpllid, llupbp, 'real*4');
    fwrite(vpbpllid, llvpbp, 'real*4');
    fclose(upbpllid);
    fclose(vpbpllid);

    %  Look at lat-lon projections of u'T' and v'T'
    pllupbp = zeros(6*510,510);
    pllvpbp = zeros(6*510,510);
    for i = 1:6
      xb = (i-1)*510 + 1;  xe = xb + 510 - 1;
      yb = 1;              ye = 510;
      pllupbp(xb:xe,yb:ye) = llupbp(:,:,i);
      pllvpbp(xb:xe,yb:ye) = llvpbp(:,:,i);
    end
    shift=-1;
    grph_CS(sq(pllupbp),xcs,ycs,xgs,ygs,[clow],[chigh],shift);
    title([ 'u''B'' at ' ...
            sprintf('%g',R(iz)) ...
            'm depth on the 510x510x6 cubesphere for 1992--2003  ["cube5"]']);
%    print('-painters', '-dpng', '-r650', ...
%          ['upTp_' sprintf('%02d_%07.1f',iz,R(iz)) 'm_650.png'])
    print('-painters', '-dpng', '-r150', ...
          ['upBp_' sprintf('%02d_%07.1f',iz,R(iz)) 'm_150.png'])

    grph_CS(sq(pllvpbp),xcs,ycs,xgs,ygs,[clow],[chigh],shift);
    title([ 'v''B'' at ' ...
            sprintf('%g',R(iz)) ...
            'm depth on the 510x510x6 cubesphere for 1992--2003  ["cube5"]']);
%    print('-painters', '-dpng', '-r650', ...
%          ['vpTp_' sprintf('%02d_%07.1f',iz,R(iz)) 'm_650.png'])
    print('-painters', '-dpng', '-r150', ...
          ['vpBp_' sprintf('%02d_%07.1f',iz,R(iz)) 'm_150.png'])

  end
  
  
  %=======================================================
  %  Calculate :  (v'B')/(dB/dz)
    
  bid   = fopen(  'b-ave-1992-2003', 'r', 'ieee-be');
  bm1id = fopen('bm1-ave-1992-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 = 15;
  iz = 16;
  for iz = 2:50,
    izm1 = iz - 1;
    disp(sprintf('iz = %d',iz));
    %  iz = 21;
    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);
    offm1 = (iz - 2)*nps*4;
    fseek(bm1id, offm1, 'bof');
    bm1 =reshape(fread(bm1id, nps, 'real*4'),ne,ne,6);

    %  ( v'B' )/( dB/dz )
    dbdz = (b - bm1)/(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('vpbpdbdz-1992-2003', 'a', 'ieee-be');
    fwrite(mid, vpbpdbdz, 'real*4');
    fclose(mid);
    
    %  Plot results
    clow  = [ -10 ];  % -20;
    chigh = [  10 ];  %  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,[clow],[chigh],shift);
    title([ '(v''B'')/(dB/dz) at ' ...
            sprintf('%g',Rmid(iz)) ...
            'm depth on the 510x510x6 cubesphere for 1992--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*pi/(24*3600)*sin(pi*yc/180));
    stress = stress .* vpbpdbdz;
    sid = fopen('stress-b-1992-2003', 'a', 'ieee-be');
    fwrite(sid, stress, 'real*4');
    fclose(sid);
    clow  = [];  % [ -1 ];
    chigh = [];  % [  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,[clow],[chigh],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.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 vpbp
  acc = zeros(nlatm1, nrm1);
  num = zeros(size(acc));
  ne = 510;
  nps = ne * ne * 6;
  zid = fopen('vpbpll-1992-2003', 'r', 'ieee-be');
  iz = 1;
  for iz = 1:50,
    disp(sprintf('iz = %d',iz));
    offset = (iz - 1)*nps*4;
    fseek(zid, offset, 'bof');
    vpbp = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
    
    for jj = 1:nlatm1
      eval( sprintf('clear inds; inds = inds%04d;',jj) );
      tmp = vpbp(inds);
      nzinds = find(tmp ~= 0.0);
      num(jj,iz) = length(nzinds);
      acc(jj,iz) = sum(tmp(nzinds));
    end
  end
  fclose(zid);
  zvpbp = acc ./ num;
  %  save zvpbp zvpbp

  %  zonally average vpbpdbdz
  acc = zeros(nlatm1, nrm1);
  num = zeros(size(acc));
  ne = 510;
  nps = ne * ne * 6;
  zid = fopen('vpbpdbdz-1992-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('stress-b-1992-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
  
  %-------  UNUSED  -----------------------
  %  zonally average u
  acc = zeros(nlatm1, 50);
  num = zeros(size(acc));
  ne = 510;
  nps = ne * ne * 6;
  zid = fopen('Ull_ave_1994--2003', 'r', 'ieee-be');
  iz = 1;
  for iz = 1:50,
    disp(sprintf('iz = %d',iz));
    offset = (iz - 1)*nps*4;
    fseek(zid, offset, 'bof');
    U = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
    jj = 22;
    for jj = 1:nlatm1
      eval( sprintf('clear inds; inds = inds%04d;',jj) );
      tmp = U(inds);
      nzinds = find(tmp ~= 0.0);
      num(jj,iz) = length(nzinds);
      acc(jj,iz) = sum(tmp(nzinds));
    end
  end
  fclose(zid);
  zonalu = acc ./ num;
  %  save zonalu zonalu
  %-------  UNUSED  -----------------------

  %  zonally average B
  acc = zeros(nlatm1, 50);
  num = zeros(size(acc));
  ne = 510;
  nps = ne * ne * 6;
  zid = fopen('b-ave-1992-2003', 'r', 'ieee-be');
  iz = 1;
  for iz = 1:50,
    disp(sprintf('iz = %d',iz));
    offset = (iz - 1)*nps*4;
    fseek(zid, offset, 'bof');
    B = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
    jj = 22;
    for jj = 1:nlatm1
      eval( sprintf('clear inds; inds = inds%04d;',jj) );
      tmp = B(inds);
      nzinds = find(tmp ~= 0.0);
      num(jj,iz) = length(nzinds);
      acc(jj,iz) = sum(tmp(nzinds));
    end
  end
  fclose(zid);
  zonalB = acc ./ num;
  %  save zonalB zonalB

  %-------  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')
  
  %  plot and print bar(U)
  load zonalu
  surf(hmid,-R,zonalu'), view(2), shading interp, colorbar
  contour(hmid,-R,mu'), colorbar
  ac = 25;
  ac = [-0.3:.005:.3];
  [c,h] = contourf(hmid,-R,mu',ac); colorbar
  hold on, contour(hmid,-R,mu',ac), hold off
  title('Zonally Averaged U for 1994--2003 [cube5]')
  xlabel('Latitude [deg]')
  zlabel('Depth [m]')
  axis([ -90 90 -6000 200 ])
  print('-painters', '-dpng', '-r150', 'Ull_94--03_r150.png')
  print('-painters', '-dpng', '-r650', 'Ull_94--03_r650.png')
  %-------  UNUSED  -----------------------

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

  clear all ; close all
  load allR
  % load zonalu
  load zvpbpdbdz
  load stress
  load zonalB
  % load dzaTdz
  load zvpbp

  % id0 = 'U';
  id1 = 'zvpbpdbdz';
  id2 = 'stress';
  id3 = 'B';
  % id4 = 'dzaTdz';
  id5 = 'zvpbp';
  % units0 = 'm/s';
  units1 = 'm^2/s';
  units2 = 'N/m^2';
  units3 = 'kg/m^3';
  % units4 = 'deg C/m';
  units5 = 'm*kg/m^3';

  !rm -f Zonal_ave_92-03.nc
  nc = netcdf('Zonal_ave_92-03.nc', 'clobber');
  nc.reference = 'Zonal averages from Dimitris "cube 5" integrations';
  nc.author = 'Ed Hill <eh3@mit.edu>';
  nc.date = 'Aug 11, 2004';
  
  nc('X') = length(hmid);
  nc('Y') = length(R);
  nc('Ym1') = length(R) - 1;
  nc{'X'}  = 'X';
  nc{'Y'}  = 'Y';
  nc{'Ym1'}  = 'Ym1';
  % nc{ id0 } = { 'Y', 'X' };
  nc{ id1 } = { 'Ym1', 'X' };
  nc{ id2 } = { 'Ym1', 'X' };
  nc{ id3 } = { 'Y', 'X' };
  % nc{ id4 } = { 'Ym1', 'X' };
  nc{ id5 } = { 'Y', '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';
  % nc{ id0 }.units = units0;
  % nc{ id0 }.long_name = id0;
  % nc{ id0 }.missing_value = ncdouble(NaN);
  % nc{ id0 }.FillValue_ = ncdouble(0.);
  nc{ id1 }.units = units1;
  nc{ id1 }.long_name = id1;
  nc{ id1 }.missing_value = ncdouble(NaN);
  nc{ id1 }.FillValue_ = ncdouble(0.);
  nc{ id2 }.units = units2;
  nc{ id2 }.long_name = id2;
  nc{ id2 }.missing_value = ncdouble(NaN);
  nc{ id2 }.FillValue_ = ncdouble(0.);
  nc{ id3 }.units = units3;
  nc{ id3 }.long_name = id3;
  nc{ id3 }.missing_value = ncdouble(NaN);
  nc{ id3 }.FillValue_ = ncdouble(0.);
  % nc{ id4 }.units = units4;
  % nc{ id4 }.long_name = id4;
  % nc{ id4 }.missing_value = ncdouble(NaN);
  % nc{ id4 }.FillValue_ = ncdouble(0.);
  nc{ id5 }.units = units5;
  nc{ id5 }.long_name = id5;
  nc{ id5 }.missing_value = ncdouble(NaN);
  nc{ id5 }.FillValue_ = ncdouble(0.);
  nc{'X'}(:) = hmid;
  nc{'Y'}(:) = -R;
  nc{'Ym1'}(:) = -Rmid;
  % nc{ id0 }(:) = permute(zonalu,[ 2 1 ]);
  nc{ id1 }(:) = permute(zvpbpdbdz,[ 2 1 ]);
  nc{ id2 }(:) = permute(stress,[ 2 1 ]);
  nc{ id3 }(:) = permute(zonalB,[ 2 1 ]);
  % nc{ id4 }(:) = permute(dzaTdz,[ 2 1 ]);
  nc{ id5 }(:) = permute(zvpbp,[ 2 1 ]);
  nc = close(nc);

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