% Surface Geostrophic Currents % Alex Wineteer % Email: wineteer@jpl.nasa.gov, awinetee@calpoly.edu %{ This code will calculate the surface geostrophic currents based on the input ocean surface heights. INPUTS: surface height folder locations (Eta), size of grid [x,y,z] -- include z, even if it is just 1. grid folder containing lats, dxC, dyC, dxG, dyG OUTPUTS: NONE, though it will write the currents to disk under Ug/Ug directories. Also note that these output binary files will be written in a pcolor-like format -- not the default MITgcm format. THis means that reading them only requires a call to readbin.m Grid Output Locations: (NOTE: actual grid output is flipped up-down, and transposed relative to this grid -- at least, in MATLAB) | | | |____^____|_________|___ | V21 | | | | | | | | |____^____|_________|___ | V11 | | | U11 -> U12 -> | | | |_________|_________|___ * %} clear %folder = ['~dmenemen/llc_2160/regions/global/Eta_daily']; siz = [8640 7560 1]; %fnames = dir([folder '/*.2*']); %numFiles = length(fnames); %currFldrVg = ['~dmenemen/llc_2160/regions/global/Vg_daily']; currFldrVg = ['/nobackupp8/awinetee/global/Vg_noIce_hamming_daily']; eval(['mkdir ' currFldrVg]); %currFldrUg = ['~dmenemen/llc_2160/regions/global/Ug_daily']; currFldrUg = ['/nobackupp8/awinetee/global/Ug_noIce_hamming_daily']; eval(['mkdir ' currFldrUg]); nx = 2160; lats = quikreadpcolor_llc('~dmenemen/llc_2160/grid/YC.data',nx); lons = quikreadpcolor_llc('~dmenemen/llc_2160/grid/XC.data',nx); landV = quikreadpcolor_llc('~dmenemen/llc_2160/grid/hFacS.data',nx); landU = quikreadpcolor_llc('~dmenemen/llc_2160/grid/hFacW.data',nx); [dx,dy] = quikreadpcolor_dxdy_llc('~dmenemen/llc_2160/grid/DXC.data','~dmenemen/llc_2160/grid/DYC.data',nx); [dxG,dyG] = quikreadpcolor_dxGdyG_llc('~dmenemen/llc_2160/grid/DXG.data','~dmenemen/llc_2160/grid/DYG.data',nx); fc = 2*7.2921*10^-5*sind(lats); g = 9.806; rho = 1027.5; SEAICE_rhoIce = 910; SEAICE_rhoSnow = 330; pn='~dmenemen/llc_2160/MITgcm/'; %fnames = dir('~dmenemen/llc_2160/regions/global/Eta_daily/E*'); fnames = dir('/nobackupp8/awinetee/global/Eta_seaIceRemoved_30hLPF_hamming_daily/E*'); %numFiles = length(fnames); strNames = {}; n = 0; for ii = 92160:(80*24):1586400; ii if (ii<1198080) pin=[pn 'run_day49_624/']; else pin=[pn 'run/']; end dyy=ts2dte(ii,45,2011,1,17,30); fnam = ['/nobackupp8/awinetee/global/Eta_seaIceRemoved_30hLPF_hamming_daily/' 'Eta_seaIceRemoved_' int2str((nx)) 'x' int2str((nx*13)) '.' dyy]; %fnamIce = [pin 'SIheff.' myint2str(ii,10) '.data']; %fnamSnow = [pin 'SIhsnow.' myint2str(ii,10) '.data']; if exist(fnam)% && exist(fnamIce) && exist(fnamSnow) n = n + 1; strNames{n} = dy; %SIheff = quikreadpcolor_llc([pin 'SIheff.' myint2str(ii,10) '.data'],nx); % SIhsnow = quikreadpcolor_llc([pin 'SIhsnow.' myint2str(ii,10) '.data'],nx); %% FIX THIS PORTION TO DO MEDIAN FILTERING FOR ICE BUT MAAKE IT DO IT FAST eta = quikreadpcolor_llc(fnam,nx); etaSTD = stdfilt(eta,ones(51)); etaMedFilt = medfilt2(eta,[25,25]); eta(etaSTD > .03&(lats>65|lats<-65)) = etaMedFilt(etaSTD > .03&(lats>65|lats<-65)); % figure % quikpcolor(eta'); % colorbar % caxis([-2 0]); %% % SIheff(isnan(SIheff)) = 0; % SIhsnow(isnan(SIhsnow)) = 0; % iceLoad = SIheff*SEAICE_rhoIce + SIhsnow*SEAICE_rhoSnow; %ssh = eta;% + iceLoad/rho; s = size(eta); ss= s; %eta = ssh; deltaEtaX = [diff(eta,1,1);eta(1,:)-eta(end,:)]; %[eta(2:end,:) - eta(1:end-1,:);eta(1,:)-eta(end,:)]; % diff between cell enters -- therefore centered on boundary deltaEtaY = [eta(:,2:end) - eta(:,1:end-1),nan*ones(s(1),1)]; ffY = [((fc(:,1:end-1) + fc(:,2:end))/2),nan*ones(s(1),1)]; ffX = [(fc(1:end-1,:) + fc(2:end,:)) / 2;(fc(end,:)+fc(1,:))/2]; Us = -(g ./ ffY .* deltaEtaY ./ dy); Vs = g ./ ffX .* deltaEtaX ./ dx; Us = Us.*landV; Vs = Vs.*landU; % switched b/c edges are switched right now % Vs(abs(lats)<15) = nan; % Us(abs(lats)<15) = nan; clear eta ffY ffX deltaEtaX deltaEtaY s= size(Us); Us1 = (Us(:,1:end-1) + Us(:,2:end))/2; Us1(:,end+1) = nan*ones(s(1),1); Umit = Us1(1:end-1,:) + (.5.*dxG(1:end-1,:)./dx(2:end,:)).*(Us1(2:end,:)-Us1(1:end-1,:)); s = size(Umit); Umit(end+1,:) = Us1(end,:) + (.5.*dxG(end,:)./dx(1,:)).*(Us1(1,:)-Us1(end,:)); Vs1 = (Vs(1:end-1,:) + Vs(2:end,:))/2; s = size(Vs1); Vs1(end+1,:) = (Vs(1,:)+Vs(end,:))/2; Vmit = Vs1(:,1:end-1) + (.5.*dyG(:,1:end-1)./dy(:,2:end)).*(Vs1(:,2:end)-Vs1(:,1:end-1)); s = size(Vmit); Vmit(:,end+1) = zeros(s(1),1); Vmit(isnan(Vmit)) = 0; Umit(isnan(Umit)) = 0; Vmit = Vmit.*landV; Umit = Umit.*landU; % land masks s = ss; % % figure % quikpcolor(Umit'); % caxis([-.2 .5]) % figure % quikpcolor(Vmit'); % caxis([-.2 .5]) %pause fileName = [currFldrUg '/Ug_noIce_hamming' int2str(s(1)) 'x' int2str(s(2)) '.' dyy]; writebin(fileName,Umit); fileName = [currFldrVg '/Vg_noIce_hamming' int2str(s(1)) 'x' int2str(s(2)) '.' dyy]; writebin(fileName,Vmit); end end %% % figure(3) % quikpcolor(sqrt(Uout.^2 + Vout.^2)'); % colorbar % caxis([0 1]) % daspect([1 1 1]) % title('Umag 27 Day Avg') % print('-dpng','-r800','MAGmit_27day_global.png') % figure(4) % quikpcolor(Vout'); % colorbar % caxis([-.3 1]) % title('Vmit 27 Day Avg') % daspect([1 1 1]) % print('-dpng','-r800','Vmit_27day_global.png')