%======================================================= % % $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/high_res_cube/eddy_flux/c20a/plot_c20a.m,v 1.8 2005/03/23 21:46:02 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_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 %================================================================== % 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 % % !echo `ls -1 ave__92_99 | sed -e 's|\.ave||g' \ % | sed -e "s|^|'|g" | sed -e "s|$|',|g"` 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 = { ... 'ETANSQ', 'SFLUX', 'SRELAX', 'TAUX', 'TAUY', 'TFLUX', ... 'TICE', 'TRELAX' }; % fid = fopen('ave__92_99/DRHODR.ave','r','ieee-be'); % tmp = fread(gid,inf,'real*4',0,'ieee-be'); % fclose(fid); ne = 510; nf = 6; nz = 50; nslab = ne*ne*nf; adir = 'ave__92_99'; lat = [-90:90]; lon = [0:360]; % ! rm -f cube_20a_*.nc ifld = 5; fields = union(fields_2d, fields_3d); for ifld = 1:length(fields) ir = [ 1 2 3 5 10 15 20 25 30 35 40 50 ]; if ismember(fields{ifld},fields_2d) ir = [ 1 ]; end disp([ ' ' fields{ifld} ' :' ]); nc = netcdf(['cube_20a_' fields{ifld} '.nc'], 'clobber'); nc.reference = [ 'Results from Dimitris Menemenlis''' ... '"cube 20a" integrations' ]; nc.author = 'Ed Hill '; nc.date = 'Feb 24, 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); fname = sprintf('%s/%s.ave',adir,fields{ifld}); fid = fopen(fname,'r','ieee-be'); id = fields{ifld}; 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 = %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 nc{ id }(ii,:,:) = permute(v,[2 1]); end fclose(fid); nc = close(nc); end % ! ncdump cube_20a.nc | more % ! rm -rf output_netCDF ; mkdir output_netCDF % ! scp output_netCDF/*.nc channel.mit.edu:/home/edhill/INGRID_PEOPLE/EH3/eddy_flux/cube_20a/ %======================================================= % 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__92_04/'; 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) 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 50 ]; xc360 = XC + 180; % ! rm -f cube_20a_primes.nc nc = netcdf(['cube_20a_primes.nc'], 'clobber'); nc.reference = [ 'Results from Dimitris Menemenlis''' ... '"cube 20a" integrations' ]; nc.author = 'Ed Hill '; nc.date = 'Feb 28, 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_20a_primes.nc channel.mit.edu:/home/edhill/INGRID_PEOPLE/EH3/eddy_flux/cube_20a/primes/ % ! mv cube_20a_primes.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)