eddy_flux/c22/plot_c22.m

%=======================================================
%
%  $Id: plot_c22.m,v 1.2 2005/03/28 21:42:31 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_22" or "c22" 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


  %==================================================================
  %  Read the tile00?.mitgrid files
  gvars = { 'XC','YC','DXF','DYF','RA','XG','YG','DXV', ...
            'DYU','RAZ','DXC','DYC','RAW','RAS','DXG','DYG' };

  ne = 510;
  nep1 = ne + 1;
  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'),[nep1,nep1,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 = [(nep1-10):nep1];  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
  is = [1:ne];
  vs = { 'XC','YC','DXF','DYF','RA' };
  for i = 1:length(vs)
    eval(sprintf('%s = %s(is,is,:);',vs{i},vs{i}));
  end

  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;
  
  %==================================================================
  %  Project fields to lower-res 1-degree Lat-Lon and write 
  %  as NetCDF for viewing with Ingrid
  %
  fields_3d = { ...
      'DRHODR', 'RHOANOSQ', 'RHOAnoma', 'SALT', 'SALTSQ', ...
      'THETA', 'THETASQ', 'URHOMASS', ...
      'USLTMASS', 'UTHMASS', 'UVEL', 'UVELMASS', 'UVELSQ', ...
      'UV_VEL_Z', 'VRHOMASS', 'VSLTMASS', 'VTHMASS', ...
      'VVEL', 'VVELMASS', 'VVELSQ', 'WRHOMASS', 'WSLTMASS', ...
      'WTHMASS', 'WU_VEL', 'WVELMASS', 'WVELSQ', 'WV_VEL' };
  fields_2d = { ...
      'ETAN', 'ETANSQ', 'EmPmRtave', 'PHIBOT', 'SFLUX', 'SRELAX', ...
      'TAUX', 'TAUY', 'TFLUX', ...
      'TICE', 'TRELAX', 'UICEtave', 'VICEtave' };
  
  ne = 510;
  nf = 6;
  nz = 50;
  nslab = ne*ne*nf;
  adir = 'ave_1992-2004';
  lat = [-90:90];
  lon = [0:360];
  ir = [ 1 2 3 5 10 15 20 25 30 35 40 50 ];
  
  %  ! rm -f cube_22_at1deg.nc
  nc = netcdf(['cube_22_at1deg.nc'], 'clobber');
  nc.reference = [ 'The cube_22 averages from Dimitris Menemenlis' ...
                   ' regridded to 1-deg Lat-Lon' ];
  nc.author = 'Ed Hill <eh3@mit.edu>';
  nc.date = 'Mar 27, 2005';
  nc('X') = length(lon);
  nc('Y') = length(lat);
  nc('Z') = length(ir);
  nc{'X'} = 'X';
  nc{'Y'} = 'Y';
  nc{'Z'} = 'Z';
  nc{'X'}.uniquename = 'X';
  nc{'X'}.long_name = 'longitude';
  nc{'X'}.gridtype = ncint(1);
  nc{'X'}.units = 'degree_east';
  nc{'Y'}.uniquename = 'Y';
  nc{'Y'}.long_name = 'latitude';
  nc{'Y'}.gridtype = ncint(0);
  nc{'Y'}.units = 'degree_north';
  nc{'Z'}.uniquename = 'Z';
  nc{'Z'}.long_name = 'depth';
  nc{'Z'}.gridtype = ncint(0);
  nc{'Z'}.units = 'm';
  nc{'X'}(:) = lon;
  nc{'Y'}(:) = lat;
  nc{'Z'}(:) = R(ir);
  
  ifld = 5;
  fields = union(fields_2d, fields_3d);
  for ifld = 1:length(fields)
    
    id = fields{ifld};
    if ismember(fields{ifld},fields_3d)
      ir = [ 1 2 3 5 10 15 20 25 30 35 40 50 ];
      nc{ id } = { 'Z' 'Y' 'X' };
    else
      ir = [ 1 ];
      nc{ id } = { 'Y' 'X' };
    end
    nc{ id }.missing_value = ncdouble(NaN);
    nc{ id }.FillValue_ = ncdouble(0.0);

    disp([ '  ' fields{ifld} ' :' ]);

    fname = sprintf('%s/%s.ave',adir,fields{ifld});
    fid = fopen(fname,'r','ieee-be');
    
    ii = 1;
    for ii = 1:length(ir)

      iz = ir(ii);
      disp(sprintf('    iz = %g',iz));
      
      fseek(fid,nslab*4*(iz-1),'bof');
      tmp = fread(fid,nslab,'real*4',0,'ieee-be');
      tr = permute(reshape(tmp,[ 510 6 510 ]),[1 3 2]);
      %  surf(tr(:,:,1)), view(2), shading interp
      xc360 = XC + 180;
      trn = tr;
      trn(find(tr == 0.0)) = NaN;
      clear tmp tr
      %  v = sdac_regrid(xc360,YC,trn,lonm,latm);
      v = ll_regrid(xc360,YC,trn,lon,lat);
      %  surf(lon,lat,v'), caxis([25 40]), view(2), shading interp, colorbar

      if length(ir) == 1
        nc{ id }(:,:) = permute(v,[2 1]);
      else
        nc{ id }(ii,:,:) = permute(v,[2 1]);
      end
      
    end
    
    fclose(fid);
  
  end

  nc = close(nc);


  %  ! ncdump cube_22_at1deg.nc | more
  %  ! scp cube_22_at1deg.nc channel.mit.edu:/home/edhill/INGRID_PEOPLE/EH3/eddy_flux/cube_22/


  %=======================================================
  %  Compute [uvw]'[tsb]'
  
  clear all
  close all
  
  gvars = { 'XC','YC','DXF','DYF','RA','XG','YG','DXV', ...
            'DYU','RAZ','DXC','DYC','RAW','RAS','DXG','DYG' };
  ne = 510;
  nep1 = ne + 1;
  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'),[nep1,nep1,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 = [(nep1-10):nep1];  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
  is = [1:ne];
  vs = { 'XC','YC','DXF','DYF','RA' };
  for i = 1:length(vs)
    eval(sprintf('%s = %s(is,is,:);',vs{i},vs{i}));
  end

  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;

  n1 = ne - 1;
  dux = zeros(size(XC));
  duy = zeros(size(XC));
  dvx = zeros(size(XC));
  dvy = zeros(size(XC));
  dux(:,:,:) = diff(XG(:,1:ne,:),1,1);
  dvx(:,:,:) = diff(XG(1:ne,:,:),1,2);
  duy(:,:,:) = diff(YG(:,1:ne,:),1,1);
  dvy(:,:,:) = diff(YG(1:ne,:,:),1,2);
  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);

  lpath = 'ave_1992-2004/';
  u__id = fopen( [lpath 'UVEL.ave'], 'r', 'ieee-be');      %   1
  v__id = fopen( [lpath 'VVEL.ave'], 'r', 'ieee-be');      %   2
  %w__id = fopen( [lpath 'WVEL.ave'], 'r', 'ieee-be');
  u2_id = fopen( [lpath 'UVELSQ.ave'], 'r', 'ieee-be');    %   3
  v2_id = fopen( [lpath 'VVELSQ.ave'], 'r', 'ieee-be');    %   4
  w2_id = fopen( [lpath 'WVELSQ.ave'], 'r', 'ieee-be');    %   5
  um_id = fopen( [lpath 'UVELMASS.ave'], 'r', 'ieee-be');  %   6
  vm_id = fopen( [lpath 'VVELMASS.ave'], 'r', 'ieee-be');  %   7
  wm_id = fopen( [lpath 'WVELMASS.ave'], 'r', 'ieee-be');  %   8
  t__id = fopen( [lpath 'THETA.ave'], 'r', 'ieee-be');     %   9
  t2_id = fopen( [lpath 'THETASQ.ave'], 'r', 'ieee-be');   %  10
  s__id = fopen( [lpath 'SALT.ave'], 'r', 'ieee-be');      %  11
  s2_id = fopen( [lpath 'SALTSQ.ave'], 'r', 'ieee-be');    %  12
  b__id = fopen( [lpath 'RHOAnoma.ave'], 'r', 'ieee-be');  %  13
  b2_id = fopen( [lpath 'RHOANOSQ.ave'], 'r', 'ieee-be');  %  14
  ut_id = fopen( [lpath 'UTHMASS.ave'], 'r', 'ieee-be');   %  15
  vt_id = fopen( [lpath 'VTHMASS.ave'], 'r', 'ieee-be');   %  16
  wt_id = fopen( [lpath 'WTHMASS.ave'], 'r', 'ieee-be');   %  17
  us_id = fopen( [lpath 'USLTMASS.ave'], 'r', 'ieee-be');  %  18
  vs_id = fopen( [lpath 'VSLTMASS.ave'], 'r', 'ieee-be');  %  19
  ws_id = fopen( [lpath 'WSLTMASS.ave'], 'r', 'ieee-be');  %  20
  ub_id = fopen( [lpath 'URHOMASS.ave'], 'r', 'ieee-be');  %  21
  vb_id = fopen( [lpath 'VRHOMASS.ave'], 'r', 'ieee-be');  %  22
  wb_id = fopen( [lpath 'WRHOMASS.ave'], 'r', 'ieee-be');  %  23
  dr_id = fopen( [lpath 'DRHODR.ave'], 'r', 'ieee-be');    %  24
  iids = [ u__id v__id u2_id v2_id w2_id um_id vm_id wm_id ... 
           t__id t2_id s__id s2_id b__id b2_id ...
           ut_id vt_id wt_id us_id vs_id ws_id ...
           ub_id vb_id wb_id dr_id ];

  %  ! rm -rf primes_92_04 ; mkdir primes_92_04
  opath = 'primes_92_04/';
  up2___id = fopen([opath 'up2'], 'wb', 'ieee-be');
  vp2___id = fopen([opath 'vp2'], 'wb', 'ieee-be');
  wp2___id = fopen([opath 'wp2'], 'wb', 'ieee-be');
  tp2___id = fopen([opath 'tp2'], 'wb', 'ieee-be');
  sp2___id = fopen([opath 'sp2'], 'wb', 'ieee-be');
  bp2___id = fopen([opath 'bp2'], 'wb', 'ieee-be');
  uptp__id = fopen([opath 'uptp'], 'wb', 'ieee-be');
  vptp__id = fopen([opath 'vptp'], 'wb', 'ieee-be');
  wptp__id = fopen([opath 'wptp'], 'wb', 'ieee-be');
  upsp__id = fopen([opath 'upsp'], 'wb', 'ieee-be');
  vpsp__id = fopen([opath 'vpsp'], 'wb', 'ieee-be');
  wpsp__id = fopen([opath 'wpsp'], 'wb', 'ieee-be');
  upbp__id = fopen([opath 'upbp'], 'wb', 'ieee-be');
  vpbp__id = fopen([opath 'vpbp'], 'wb', 'ieee-be');
  wpbp__id = fopen([opath 'wpbp'], 'wb', 'ieee-be');
  vpbpdzid = fopen([opath 'vpbp_dbdz'], 'wb', 'ieee-be');
  str___id = fopen([opath 'stress'], 'wb', 'ieee-be');

  comm = [ 'permute(reshape(fread( id ,nslab,''real*4'',0,' ...
           '''ieee-be''),[ne 6 ne]),[1 3 2]);' ];
  readslab = inline(comm,'id','nslab','ne');
  
  ne = 510;
  nslab = 6 * ne * ne;
  nztot = 50;
  iz = 1;
  for iz = 1:nztot

    disp(sprintf('  iz = %d',iz)); 
    offset = (iz - 1)*nslab*4;
    for iid = 1:length(iids)
      fseek(iids(iid), offset, 'bof');
    end
    u  = readslab(u__id,nslab,ne);   u2 = readslab(u2_id,nslab,ne);
    v  = readslab(v__id,nslab,ne);   v2 = readslab(v2_id,nslab,ne);
    um = readslab(um_id,nslab,ne);
    vm = readslab(vm_id,nslab,ne);
    wm = readslab(wm_id,nslab,ne);
    %  surf(v2(:,:,1)), view(2), shading interp
    if (iz < nztot)
      wmp1 = readslab(wm_id,nslab,ne);
    else
      wmp1 = zeros(size(wm));
    end
    wmp05 = (wm + wmp1)/2.0;
    w2 = readslab(w2_id,nslab,ne);
    t  = readslab(t__id,nslab,ne);   t2 = readslab(t2_id,nslab,ne);
    s  = readslab(s__id,nslab,ne);   s2 = readslab(s2_id,nslab,ne);
    b  = readslab(b__id,nslab,ne);   b2 = readslab(b2_id,nslab,ne);

    %  "simple squared" quantities
    up2 = u2 - u.^2;
    vp2 = v2 - v.^2;
    wp2 = w2 - wm.^2;
    fwrite(up2___id, up2, 'real*4');
    fwrite(vp2___id, vp2, 'real*4');
    fwrite(wp2___id, wp2, 'real*4');
    clear up2 vp2 wp2 u2 v2 w2
    tp2 = t2 - t.^2;
    sp2 = s2 - s.^2;
    bp2 = b2 - b.^2;
    fwrite(tp2___id, tp2, 'real*4');
    fwrite(sp2___id, sp2, 'real*4');
    fwrite(bp2___id, bp2, 'real*4');
    clear tp2 sp2 bp2 t2 s2 b2
    ut = readslab(ut_id,nslab,ne);   vt = readslab(vt_id,nslab,ne);
    us = readslab(us_id,nslab,ne);   vs = readslab(vs_id,nslab,ne);
    ub = readslab(ub_id,nslab,ne);   vb = readslab(vb_id,nslab,ne);
    drhodr = readslab(dr_id,nslab,ne);
    [ tonu, tonv ] = mass_on_u_v(t);
    [ sonu, sonv ] = mass_on_u_v(s);
    [ bonu, bonv ] = mass_on_u_v(b);
    uptp = ut - u .* tonu;    vptp = vt - v .* tonv;
    upsp = us - u .* sonu;    vpsp = vs - v .* sonv;
    upbp = ub - u .* bonu;    vpbp = vb - v .* bonv;

    %  llupbp = upbp .* llux + vpbp .* llvx;
    llvpbp = upbp .* lluy + vpbp .* llvy;
    if iz > 1
      %  ave_llupbp = (llupbp + old_llupbp)/2.0;
      ave_llvpbp = (llvpbp + old_llvpbp)/2.0;
      tmp = drhodr;
      tmp(find(tmp == 0.0)) = 1.0;
      vpbpdbdz = ave_llvpbp ./ tmp;
      fwrite(vpbpdzid, vpbpdbdz, 'real*4');
      fac = 1000 * 2*pi/(24*3600);
      stress = fac * sin(pi*YC/180) .* vpbpdbdz;
      fwrite(str___id, stress, 'real*4');
    end
    %  old_llupbp = llupbp;
    old_llvpbp = llvpbp;

    fwrite(uptp__id, uptp, 'real*4');
    fwrite(vptp__id, vptp, 'real*4');
    fwrite(upsp__id, upsp, 'real*4');
    fwrite(vpsp__id, vpsp, 'real*4');
    fwrite(upbp__id, upbp, 'real*4');
    fwrite(vpbp__id, vpbp, 'real*4');
    clear uptp vptp upsp vpsp upbp vpbp vpbpdbdz
    wt = readslab(wt_id,nslab,ne);    wptp = wt - wmp05 .* t;
    ws = readslab(ws_id,nslab,ne);    wpsp = ws - wmp05 .* s;
    wb = readslab(wb_id,nslab,ne);    wpbp = wb - wmp05 .* b;
    fwrite(wptp__id, wptp, 'real*4');
    fwrite(wpsp__id, wpsp, 'real*4');
    fwrite(wpbp__id, wpbp, 'real*4');
  end

  ne = 510;
  nz = 1;
  nr = 50;
  nrm1 = nr - 1;
  nlat = 181;  nlatm1 = nlat - 1;
  %  save indicies for zonal averages
  hvals = linspace(-90,90,nlat);
  i = 2;
  for i = 2:nlat
    inds = find(hvals(i-1)<YC & YC<hvals(i));
    comm = sprintf('inds%04d = uint32(inds);',i-1);
    eval(comm);
  end

  comm = [ 'reshape(fread( id ,nslab,''real*4'',0,' ...
           '''ieee-be''),[ne ne 6]);' ];
  readcubelev = inline(comm,'id','nslab','ne');

  %  Zonally average the ll_vpbp
  acc = zeros(nlatm1, nrm1);
  num = zeros(size(acc));
  zidu = fopen('primes_92_04/upbp', 'r', 'ieee-be');
  zidv = fopen('primes_92_04/vpbp', 'r', 'ieee-be');
  for iz = 1:50,
    disp(sprintf('iz = %d',iz));
    fseek(zidu, (iz - 1)*(ne*ne*6)*4, 'bof');
    fseek(zidv, (iz - 1)*(ne*ne*6)*4, 'bof');
    upbp = readcubelev(zidu,nslab,ne);
    vpbp = readcubelev(zidv,nslab,ne);
    llupbp = upbp .* llux + vpbp .* llvx;
    llvpbp = upbp .* lluy + vpbp .* llvy;
    %  surf(llvpbp(:,:,1)), view(2), shading interp, caxis([-.1 .1])
    for jj = 1:nlatm1
      eval( sprintf('clear inds; inds = inds%04d;',jj) );
      tmp = llvpbp(inds);
      nzinds = find(tmp ~= 0.0);
      num(jj,iz) = length(nzinds);
      acc(jj,iz) = sum(tmp(nzinds));
    end
  end
  fclose(zidu);
  fclose(zidv);
  llzvpbp = acc ./ num;
  %  surf(llupbp(:,:,1)'), view(2), shading interp, caxis([-.1 .1])
  %  surf(llvpbp(:,:,1)'), view(2), shading interp, caxis([-.1 .1])
  %  surf(llzvpbp'), view(2), shading interp, caxis([-.1 .1])
  %  save primes_92_04/za_llvpbp.mat llzvpbp llupbp llvpbp

  %  zonally average ll_vpbp_dbdz
  acc = zeros(nlatm1, nrm1);
  num = zeros(size(acc));
  zid = fopen('primes_92_04/vpbp_dbdz', 'r', 'ieee-be');
  iz = 1;
  for iz = 1:49,
    disp(sprintf('iz = %d',iz));
    fseek(zid, (iz - 1)*(ne*ne*6)*4, 'bof');
    vpbpdbdz = readcubelev(zid,nslab,ne);
    %  surf(vpbpdbdz(:,:,1)), view(2), shading interp, caxis([-50 50])
    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);
  za_ll_vpbp_dbdz = acc ./ num;
  %  save primes_92_04/za_ll_vpbp_dbdz.mat za_ll_vpbp_dbdz 
  
  %  zonally average stress
  acc = zeros(nlatm1, nrm1);
  num = zeros(size(acc));
  zid = fopen('primes_92_04/stress', 'r', 'ieee-be');
  iz = 1;
  for iz = 1:49,
    disp(sprintf('iz = %d',iz));
    fseek(zid, (iz - 1)*(ne*ne*6)*4, 'bof');
    stress = readcubelev(zid,nslab,ne);
    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);
  za_stress = acc ./ num;
  %  save primes_92_04/stress.mat  za_stress
  
  %  Average the stress over SSH-contours (streamlines) 
  eta_id = fopen('ave_1992-2004/ETAN.ave', 'r', 'ieee-be');
  comm = [ 'permute(reshape(fread( id ,nslab,''real*4'',0,' ...
           '''ieee-be''),[ne 6 ne]),[1 3 2]);' ];
  readslab = inline(comm,'id','nslab','ne');
  etan = readslab(eta_id,nslab,ne);
  fclose(eta_id);
  %  surf(YC(:,:,1)'), view(2),shading interp,colorbar
  %  surf(etan(:,:,1)'), view(2),shading interp,colorbar
  nssh = 50;
  sshvals = linspace(-2.3,1.2,nssh);
  i = 2;
  for i = 2:nssh
    inds = find(sshvals(i-1)<etan & etan<sshvals(i) & YC<-20);
    comm = sprintf('sshinds%04d = uint32(inds);',i-1);
    eval(comm);
  end
  zid = fopen('primes_92_04/stress', 'r', 'ieee-be');
  iz = 1;
  for iz = 1:49,
    disp(sprintf('iz = %d',iz));
    fseek(zid, (iz - 1)*(ne*ne*6)*4, 'bof');
    stress = readcubelev(zid,nslab,ne);
    for jj = 1:(nssh-1)
      eval( sprintf('clear inds; inds = sshinds%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);
  ssha_stress = acc ./ num;
  %  surf(log(ssha_stress')), view(2),shading interp,colorbar
  %  save primes_92_04/ssha_stress.mat  ssha_stress


  clear all
  close all
  
  %==================================================================
  %  Read the tile00?.mitgrid files
  gvars = { 'XC','YC','DXF','DYF','RA','XG','YG','DXV', ...
            'DYU','RAZ','DXC','DYC','RAW','RAS','DXG','DYG' };

  ne = 510;
  nep1 = ne + 1;
  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'),[nep1,nep1,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 = [(nep1-10):nep1];  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
  is = [1:ne];
  vs = { 'XC','YC','DXF','DYF','RA' };
  for i = 1:length(vs)
    eval(sprintf('%s = %s(is,is,:);',vs{i},vs{i}));
  end

  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;

  n1 = ne - 1;
  dux = zeros(size(XC));
  duy = zeros(size(XC));
  dvx = zeros(size(XC));
  dvy = zeros(size(XC));
  dux(:,:,:) = diff(XG(:,1:ne,:),1,1);
  dvx(:,:,:) = diff(XG(1:ne,:,:),1,2);
  duy(:,:,:) = diff(YG(:,1:ne,:),1,1);
  dvy(:,:,:) = diff(YG(1:ne,:,:),1,2);
  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);

  %==================================================================
  %  Project fields to lower-res 1-degree Lat-Lon and write 
  %  as NetCDF for viewing with Ingrid
  %
  ne = 510;
  nf = 6;
  nz = 50;
  nslab = ne*ne*nf;
  adir = 'primes_92_04';
  lat = [-90:90];
  lon = [0:360];
  ir = [ 1 2 3 5 10 15 20 25 30 35 40 ];
  xc360 = XC + 180;
  
  %  ! rm -f cube_22_primes_at1deg.nc
  nc = netcdf(['cube_22_primes_at1deg.nc'], 'clobber');
  nc.reference = [ 'The cube_22 primes from Dimitris Menemenlis' ...
                   ' regridded to 1-deg Lat-Lon' ];
  nc.author = 'Ed Hill <eh3@mit.edu>';
  nc.date = 'March 27, 2005';
  nc('X') = length(lon);
  nc('Y') = length(lat);
  nc('Z') = length(ir);
  nc{'X'} = 'X';
  nc{'Y'} = 'Y';
  nc{'Z'} = 'Z';
  nc{'X'}.uniquename = 'X';
  nc{'X'}.long_name = 'longitude';
  nc{'X'}.gridtype = ncint(1);
  nc{'X'}.units = 'degree_east';
  nc{'Y'}.uniquename = 'Y';
  nc{'Y'}.long_name = 'latitude';
  nc{'Y'}.gridtype = ncint(0);
  nc{'Y'}.units = 'degree_north';
  nc{'Z'}.uniquename = 'Z';
  nc{'Z'}.long_name = 'depth';
  nc{'Z'}.gridtype = ncint(0);
  nc{'Z'}.units = 'm';
  nc{'X'}(:) = lon;
  nc{'Y'}(:) = lat;
  nc{'Z'}(:) = R(ir);
  
  f_s_3d = { {'tp2'}, {'sp2'}, {'bp2'}, {'vpbp_dbdz'} };
  
  ifg = 1;
  for ifg = 1:length(f_s_3d)

    acell = f_s_3d{ifg};
    tname = acell{1};
    disp([ '  ' tname ' :' ]);

    fname = sprintf('%s/%s',adir,tname);
    fid = fopen(fname,'r','ieee-be');
    
    id = tname;
    nc{ id } = { 'Z' 'Y' 'X' };
    nc{ id }.missing_value = ncdouble(NaN);
    nc{ id }.FillValue_ = ncdouble(0.0);

    ii = 1;
    for ii = 1:length(ir)

      iz = ir(ii);
      disp(sprintf('    iz = %3d   R(iz) = %g',iz,R(iz)));
      
      fseek(fid,nslab*4*(iz-1),'bof');
      tmp = fread(fid,nslab,'real*4',0,'ieee-be');
      tr = reshape(tmp,[ 510 510 6 ]);
      %  surf(tr(:,:,1)), view(2), shading interp
      trn = tr;
      trn(find(tr == 0.0)) = NaN;
      clear tmp tr
      %  v = sdac_regrid(xc360,YC,trn,lonm,latm);
      v = ll_regrid(xc360,YC,trn,lon,lat);
      %  surf(lon,lat,v'), caxis([25 40]), view(2), shading interp, colorbar

      nc{ id }(ii,:,:) = permute(v,[2 1]);
      
    end
    
    fclose(fid);
  
  end

  id = 'sum_up2_vp2';
  nc{ id } = { 'Z' 'Y' 'X' };
  nc{ id }.missing_value = ncdouble(NaN);
  nc{ id }.FillValue_ = ncdouble(0.0);
  fidu = fopen(sprintf('%s/%s',adir,'up2'),'r','ieee-be');
  fidv = fopen(sprintf('%s/%s',adir,'vp2'),'r','ieee-be');
  for ii = 1:length(ir)
    iz = ir(ii);
    disp(sprintf('    iz = %3d   R(iz) = %g',iz,R(iz)));
    fseek(fidu,nslab*4*(iz-1),'bof');
    fseek(fidv,nslab*4*(iz-1),'bof');
    tru = reshape(fread(fidu,nslab,'real*4',0,'ieee-be'),[510 510 6]);
    trv = reshape(fread(fidv,nslab,'real*4',0,'ieee-be'),[510 510 6]);
    trnu = tru;  trnu(find(tru == 0.0)) = NaN;
    trnv = trv;  trnv(find(trv == 0.0)) = NaN;
    clear tmp tru trv
    lluv = ll_regrid(xc360,YC,trnu+trnv,lon,lat);
    nc{ id }(ii,:,:) = permute(lluv,[2 1]);
  end  
  fclose(fidu);
  fclose(fidv);

  f_v_3d = { {'up2','vp2'}, ...
             {'uptp','vptp'}, {'upsp','vpsp'}, {'upbp','vpbp'} };
  for ip = 1:length(f_v_3d)
    cell = f_v_3d{ip};
    idu = cell{1};
    idv = cell{2};
    disp(['  ' idu '  ' idv]);
    nc{ idu } = { 'Z' 'Y' 'X' };
    nc{ idu }.missing_value = ncdouble(NaN);
    nc{ idu }.FillValue_ = ncdouble(0.0);
    nc{ idv } = { 'Z' 'Y' 'X' };
    nc{ idv }.missing_value = ncdouble(NaN);
    nc{ idv }.FillValue_ = ncdouble(0.0);
    fidu = fopen(sprintf('%s/%s',adir,idu),'r','ieee-be');
    fidv = fopen(sprintf('%s/%s',adir,idv),'r','ieee-be');
    for ii = 1:length(ir)
      iz = ir(ii);
      disp(sprintf('    iz = %3d   R(iz) = %g',iz,R(iz)));
      fseek(fidu,nslab*4*(iz-1),'bof');
      fseek(fidv,nslab*4*(iz-1),'bof');
      tru = reshape(fread(fidu,nslab,'real*4',0,'ieee-be'),[510 510 6]);
      trv = reshape(fread(fidv,nslab,'real*4',0,'ieee-be'),[510 510 6]);
      trnu = tru;  trnu(find(tru == 0.0)) = NaN;
      trnv = trv;  trnv(find(trv == 0.0)) = NaN;
      clear tmp tru trv
      llru = trnu .* llux + trnv .* llvx;
      llrv = trnu .* lluy + trnv .* llvy;
      llu = ll_regrid(xc360,YC,llru,lon,lat);
      llv = ll_regrid(xc360,YC,llrv,lon,lat);
      nc{ idu }(ii,:,:) = permute(llu,[2 1]);
      nc{ idv }(ii,:,:) = permute(llv,[2 1]);
    end
  end
  fclose(fidu);
  fclose(fidv);
  
  nc = close(nc);

  %  ! scp cube_22_primes_at1deg.nc  channel.mit.edu:/home/edhill/INGRID_PEOPLE/EH3/eddy_flux/cube_22/
  %  ! mv cube_22_primes_at1deg.nc primes_92_04


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

  
  %=======================================================
  %  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  -----------------------

  