%======================================================= % % $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/high_res_cube/eddy_flux/plot_UVB.m,v 1.1 2004/08/11 15:49:15 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 12-yr averages from the 1-yr averages clear all ; close all ne = 510; ntot = ne*ne*6 * 50; by = 1992; ey = 2003; p1 = 'U_ave_%4.4d'; p2 = 'V_ave_%4.4d'; p3 = 'b-ave-%4.4d'; p4 = 'bm1-ave-%4.4d'; p5 = 'ub-ave-%4.4d'; p6 = 'vb-ave-%4.4d'; ft = 1; for ft = 1:6 sum = zeros(ntot,1); eval(sprintf('pat = p%d;',ft)); disp(sprintf('pat = "%s"',pat)); iy = by; for iy = by:ey disp(sprintf(' iy = %d',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 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)