%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Purpose: % % Calculates the wind power input, including mean and eddy terms % This is one of several routines developed for the work that culminated in % Scott and Xu, DSR 2008 % % Uses: Peruse this routine to figure out "exactly" what we did! % Learn some effecient methods to analyze hundreds of data files. % % Many different data products for wind stress and ocean currents (mean and % anomaly) were compared. Data choices are made through flags. % Using this single routine ensured uniformity of method. Automatically % generating output file names based upon the flags reduces the risk of % user error. % This won't run on your machine because it won't find the appropriate data % files. % Details: % NO adjustment of NCEP stress for ocean current effect % see get_u_tau1_*.m for this effect % When using QuikSCAT, mask out < 52 weeks of obs % because near ice edge QuikSCAT has problems %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Coded by Rob Scott, 2007 through 2008 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Choose parameters and data flag_noise = 0; %add random noise %flag_smooth = 0; % don't smooth %flag_smooth = 1; % smooth flag_smooth = 0; % flag_gen_ave = 0; % don't generate flag_gen_ave = 1; % generate ave %------------------ % currents: %------------------ % Allow for switching between: % flag = 0 % Maximenko MDT and AVISO SSHA % flag = 1; % AVISO absolute currents upd % flag = 2; % AVISO absolute currents ref % flag = 3 % GRACE-Tellus MDT and AVISO u_g,v_g anomalies % flag = 4; % AVISO SSHA with Maximenko % flag = 5; % AVISO SSHA TP/J1 only reference with Rio05 mean flag = 5; %------------------ % winds: %------------------ % flag_stress = 1; % use NCEP1 winds % flag_stress = 2; % use NCEP2 winds % flag_stress = 3; % use QuikSCAT Liu new release % flag_stress = 4; % use GSSFT2 % flag_stress = 5; % use ERA-40 % flag_stress = 6; % use QiukSCAT Large and Pond flag_stress = 3; %------------------ % dates: %------------------ % flag_date = 0; % 2000-2005 inclusive % flag_date = 2; % 2000-2005 inclusive, scrambled stress dates % flag_date = 1; % 1994-1999 inclusive flag_date = 0; %------------------------------------------------------------------------- % current data switch flag_date case 0 filelist = 'ncfile_names_2000_2005.txt'; date_range = '20000105_20051228'; case 1 filelist = 'ncfile_names_1994_1999.txt'; date_range = '19940105_19991229'; case 2 filelist = 'ncfile_names_2000_2005_scrambled.txt'; date_range = '20000105_20051228'; end switch flag case 0 disp('Using Nikolais MDT ') load('/disk/cg2/rscott/Mean_Dynamic_Topography/Nikolai/mdt_nikolai_aviso_grid.mat'); [ubar,vbar] = ssh2vel(mdt,NbLongitudes,NbLatitudes); disp('Using upd AVISO anomalies u,v') filenamelen = 37; dir = '/disk/staff/rscott/Aviso/dt_upd_msla_UV'; outfile_fin = strcat('Pin_u_tau_',date_range,'_stress',int2str(flag_stress),... '_Nikolai_AVISO_upd_anom_smooth',int2str(flag_smooth),... '_noise',int2str(flag_noise),'.mat') outfile_uvbar = strcat('uv_bar_',date_range,... '_Maximenko_mean_AVISO_upd_anom_smooth',int2str(flag_smooth),... '_noise',int2str(flag_noise),'.mat') case 1 disp('Using updated AVISO absolute u,v') filenamelen = 37; % dir = '/disk/cg2/rscott/dt_upd_madt_UV'; dir = '/disk/staff/rscott/Aviso/dt_upd_madt_UV'; outfile_fin = strcat('Pin_u_tau_',date_range,'_stress',int2str(flag_stress),... '_aviso_upd_smooth',int2str(flag_smooth),... '_noise',int2str(flag_noise),'.mat') outfile_uvbar = strcat('uv_bar_',date_range,'_upd_smooth',int2str(flag_smooth),... '_noise',int2str(flag_noise),'.mat') disp('Using ref 2-satellite AVISO absolute u,v') filenamelen = 37; dir = '/disk/staff/rscott/Aviso/dt_ref_madt_UV'; outfile_fin = strcat('Pin_u_tau_',date_range,'_stress',int2str(flag_stress),... '_aviso_ref_smooth',int2str(flag_smooth),'.mat') outfile_uvbar = strcat('uv_bar_',date_range,'_ref_smooth',int2str(flag_smooth),'.mat') case 3 disp('Using upd AVISO anomalies u,v') filenamelen = 37; dir = '/disk/staff/rscott/Aviso/dt_upd_msla_UV'; outfile_fin = strcat('Pin_u_tau_',date_range,'_stress',int2str(flag_stress),... '_GT_mean_AVISO_upd_anom_sm',int2str(flag_smooth),... '_noise',int2str(flag_noise),'.mat') outfile_uvbar = strcat('uv_bar_',date_range,... '_GT_mean_AVISO_upd_anom_smooth',int2str(flag_smooth),... '_noise',int2str(flag_noise),'.mat') load ../Mean_Dynamic_Topography/Grace_1yr/GGM02_uv_aviso_grid_v6.mat case 4 disp('Using Nikolais MDT ') load('/disk/cg2/rscott/Mean_Dynamic_Topography/Nikolai/mdt_nikolai_aviso_grid.mat'); % [ubar,vbar] = ssh2vel(mdt,NbLongitudes,NbLatitudes); disp('Using upd AVISO ssh anomalies') filenamelen = 36; dir = '/disk/cg2/rscott/dt_upd_msla_ssha'; outfile_fin = strcat('Pin_u_tau_',date_range,'_stress',int2str(flag_stress),... '_Nikolai_AVISO_upd_anom_smooth',int2str(flag_smooth),... '_noise',int2str(flag_noise),'.mat') outfile_uvbar = strcat('uv_bar_',date_range,... '_Maximenko_mean_AVISO_upd_anom_smooth',int2str(flag_smooth),... '_noise',int2str(flag_noise),'.mat') case 5 disp('Using AVISO Rio05 MDT ') % load('/disk/cg2/rscott/Mean_Dynamic_Topography/Nikolai/mdt_nikolai_aviso_grid.mat'); load('/disk/cg2/rscott/Mean_Dynamic_Topography/Aviso/mdt_rio05_aviso_grid.mat') mdt = dtarray_aviso; disp('Using ref AVISO ssh anomalies from TP & J1') filenamelen = 36; dir = '/disk/staff/rscott/Aviso/dt_ref_msla_ssha_tpj1'; outfile_fin = strcat('Pin_u_tau_',date_range,'_stress',int2str(flag_stress),... '_AVISO_TPJ1ref_anom_smooth',int2str(flag_smooth),... '_noise',int2str(flag_noise),'.mat') outfile_uvbar = strcat('uv_bar_',date_range,... '_Rio05_mean_AVISO_TPJ1ref_smooth',int2str(flag_smooth),... '_noise',int2str(flag_noise),'.mat') otherwise error('set flag') end % stress data % stress data switch flag_stress case 1 disp('Bad idea: using NCEP1 /NCAR reanalysis ') load('NCEP/stress_wkly_20000105_20051228_ncep.mat') load('NCEP/u10m_wkly_20000105_20051228_ncep.mat') case 2 disp('Good idea: using NCEP 2 DOE reanalysis') load('/disk/staff/rscott/NCEP2/stress_wkly_20000105_20051228_ncep2.mat') load('/disk/staff/rscott/NCEP2/u10m_wkly_20000105_20051228_ncep2.mat') case 3 disp('Very good idea: using QuikSCAT new release Liu algorithm') load('/disk/staff/rscott/QuikSCAT/Stress/stress_Liu_6yrave.mat') lon = (360/1440)*([0:1439]'+0.5); lat = (180/720) *([0:719]'+0.5)-90; case 4 disp('Interesting choice: using GSSFT2 -- SSM/I data') %stress_GSSFT2_1994_1999_ave.mat' load('/disk/staff/rscott/GSSFT2/stress_GSSFT2_wkly.mat') case 5 disp('Interesting choice: using ERA-40 data') load('/disk/staff/rscott/ERA40/stress_weekly_era40_1994_1999.mat') case 6 load('/disk/staff/rscott/QuikSCAT/Stress/stress_Liu_6yrave.mat') count_tau_tot = count_tot; clear taux_bar tauy_bar count_tot disp('Very good idea: using QuikSCAT new release Large Pond algorithm') dir_stress = '/disk/staff/rscott/QuikSCAT/Winds/LargePond_stress'; lon = (360/1440)*([0:1439]'+0.5); lat = (180/720) *([0:719]'+0.5)-90; otherwise error('choose flag_stress flag') end outfile_tau = strcat('tau_bar_stress',int2str(flag_stress),'.mat') lon_tau = lon; lat_tau = lat; clear lon lat pause(15) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % don't change below this line: % initial arrays utaux_bar = zeros(915,1080); vtauy_bar = zeros(915,1080); utaux_sqbar = zeros(915,1080); vtauy_sqbar = zeros(915,1080); ngood = zeros(915,1080); unity = ones(915,1080); if flag_gen_ave == 1 u_bar = zeros(915,1080); v_bar = zeros(915,1080); uv_count = zeros(915,1080); taux_bar = zeros(915,1080); tauy_bar = zeros(915,1080); n_tau = zeros(915,1080); end nx = 1082; ny = 915; % geophysical constants: omega = 2*pi/(3600*23+56*60); g=9.8; a = 6371e3*2*pi/360; % all file names: files = textread(strcat(dir,'/',filelist),'%s'); nfiles = size(files,1); % get first file % dt_upd_global_merged_madt_uv_20051214_20051214_20060519.nc infile = char(files(1)) outfile = strcat(infile(1:filenamelen),'.mat'); %getnetcdffile(strcat(dir,'/',infile),strcat(dir,'/',outfile)) if flag < 5 load(strcat(dir,'/',outfile)) else NbLongitudes = getnc(strcat(dir,'/',infile),'NbLongitudes'); NbLatitudes = getnc(strcat(dir,'/',infile),'NbLatitudes'); end % make target grid % will use for regridding stress below [xi,yi] = meshgrid(NbLongitudes,NbLatitudes); [Lon,Lat] = meshgrid(NbLongitudes,NbLatitudes); [ny,nx]=size(Lon); lon_periodic = [Lon(:,nx)-360*ones(ny,1) Lon 360*ones(ny,1)+Lon(:,1)]; geophysical_constants dx = mperdeg*cos(pi*Lat/180) .* (lon_periodic(:,3:nx+2) - lon_periodic(:,1:nx))/2; dA = dx.^2; % loop over files for n = 1:nfiles infile = char(files(n)) outfile = strcat(infile(1:filenamelen),'.mat'); % getnetcdffile(strcat(dir,'/',infile),strcat(dir,'/',outfile)) if flag < 5 load(strcat(dir,'/',outfile)) else Grid_0001 = getnc(strcat(dir,'/',infile),'Grid_0001'); end switch flag case {1,2} % convert cm/s --> m/s: orient so first index corresponds to latitude: u = Grid_0001'/100; u(u>1e17) = NaN; v = Grid_0002'/100; v(v>1e17) = NaN; case {0,3} % convert cm/s --> m/s: orient so first index corresponds to latitude: u = Grid_0001'/100; u(u>1e17) = NaN; v = Grid_0002'/100; v(v>1e17) = NaN; u = u + ubar; v = v + vbar; case {4,5} % convert cm --> m: orient so first index corresponds to latitude: ssha = Grid_0001'/100; ssha(ssha>1e17) = NaN; noise = flag_noise * randn(size(ssha)); [u,v] = ssh2vel(mdt+ssha+noise,NbLongitudes,NbLatitudes); otherwise error('step flag dumbie') end if flag_smooth == 1 u = smooth_2d(u); v = smooth_2d(v); end if flag_gen_ave == 1 guv = find(~isnan(u+v)); u_bar(guv) = u_bar(guv)+ u(guv); v_bar(guv) = v_bar(guv)+ v(guv); uv_count(guv) = uv_count(guv) + unity(guv); end %---------------------------------------------------------------- % get the stress switch flag_stress case {1,2,5} taux = squeeze(taux_wkly(n,:,:)); tauy = squeeze(tauy_wkly(n,:,:)); tau = sqrt(taux.^2 + tauy.^2); nx_ncep = length(lon_tau); lon_cyc = [ lon_tau(nx_ncep)-360; lon_tau; 360+lon_tau(1)]; taux = [taux(:,nx_ncep) taux taux(:,1)]; tauy = [tauy(:,nx_ncep) tauy tauy(:,1)]; taux_aviso = interp2(lon_cyc,lat_tau,taux,NbLongitudes,NbLatitudes); tauy_aviso = interp2(lon_cyc,lat_tau,tauy,NbLongitudes,NbLatitudes); case 3 % QuikSCAT Liu % load weekly stress file load(strcat('/disk/staff/rscott/QuikSCAT/Stress/tau_Liu_wk_',int2str(n),'.mat')) tau_wk = NaN*ones(720,1440); tau_wk(good) = taux; tau_wk(count_tot<52) = NaN; % near ice edge, has problems taux_aviso = regrid_map(lon_tau,lat_tau,tau_wk,NbLongitudes,NbLatitudes); tau_wk(good) = tauy; tau_wk(count_tot<52) = NaN; % near ice edge, has problems tauy_aviso = regrid_map(lon_tau,lat_tau,tau_wk,NbLongitudes,NbLatitudes); case 4 % GSSFT2 % regrid weekly stress field taux = squeeze(taux_fixed(n,:,:)); taux_aviso = regrid_map(lon_tau,lat_tau,taux,NbLongitudes,NbLatitudes); tauy = squeeze(tauy_fixed(n,:,:)); tauy_aviso = regrid_map(lon_tau,lat_tau,tauy,NbLongitudes,NbLatitudes); case 6 % QuikSCAT Large Pond % load weekly stress file load(strcat(dir_stress,'/tau_LP_wk_',int2str(n),'.mat')) tau_wk = NaN*ones(720,1440); tau_wk(good) = taux; tau_wk(count_tau_tot<52) = NaN; % near ice edge, has problems taux_aviso = regrid_map(lon_tau,lat_tau,tau_wk,NbLongitudes,NbLatitudes); tau_wk(good) = tauy; tau_wk(count_tau_tot<52) = NaN; % near ice edge, has problems tauy_aviso = regrid_map(lon_tau,lat_tau,tau_wk,NbLongitudes,NbLatitudes); otherwise error('set stress flag') end if flag_gen_ave == 1 good = find(~isnan(taux_aviso+tauy_aviso)); taux_bar(good) = taux_bar(good) + taux_aviso(good); tauy_bar(good) = tauy_bar(good) + tauy_aviso(good); n_tau(good) = n_tau(good) + unity(good); end % get Pin good = find(~isnan(u + v + taux_aviso + tauy_aviso)); ngood(good) = ngood(good) + unity(good); utaux_bar(good) = utaux_bar(good) + u(good).*taux_aviso(good); vtauy_bar(good) = vtauy_bar(good) + v(good).*tauy_aviso(good); utaux_sqbar(good) = utaux_sqbar(good) + (u(good).*taux_aviso(good)).^2; vtauy_sqbar(good) = vtauy_sqbar(good) + (v(good).*tauy_aviso(good)).^2; Pin_flux_den_ts(n) = sum( (u(good).*taux_aviso(good) + ... v(good).*tauy_aviso(good)).*dA(good))/sum(dA(good)); end utaux_bar(ngood>0) = utaux_bar(ngood>0)./ngood(ngood>0); vtauy_bar(ngood>0) = vtauy_bar(ngood>0)./ngood(ngood>0); % mean currents: if flag_gen_ave u_bar(uv_count>0) = u_bar(uv_count>0) ./uv_count(uv_count>0); v_bar(uv_count>0) = v_bar(uv_count>0) ./uv_count(uv_count>0); save(outfile_uvbar,'NbLongitudes','NbLatitudes','u_bar','v_bar','uv_count') taux_bar(n_tau>0) = taux_bar(n_tau>0)./n_tau(n_tau>0); tauy_bar(n_tau>0) = tauy_bar(n_tau>0)./n_tau(n_tau>0); save(outfile_tau,'NbLongitudes','NbLatitudes','taux_bar','tauy_bar','n_tau') end [Lon,Lat] = meshgrid(NbLongitudes,NbLatitudes); [ny,nx]=size(Lon); lon_periodic = [Lon(:,nx)-360*ones(ny,1) Lon 360*ones(ny,1)+Lon(:,1)]; geophysical_constants dx = mperdeg*cos(pi*Lat/180) .* (lon_periodic(:,3:nx+2) - lon_periodic(:,1:nx))/2; dA = dx.^2; save(outfile_fin,'NbLongitudes','NbLatitudes','utaux_bar','vtauy_bar','ngood',... 'dA','xi','yi','Pin_flux_den_ts') %good = find(abs(yi)>3 & ngood>0); %pinx = sum(utaux_bar(good).*dA(good).*uv_count(good)/313) %piny = sum(vtauy_bar(good).*dA(good).*uv_count(good)/313) good = find(ngood>0); pinx = sum(utaux_bar(good).*dA(good)) piny = sum(vtauy_bar(good).*dA(good)) Ain = sum(dA(good)) Pin = pinx + piny Pin/Ain load('/disk/cg2/rscott/Bottom_Topography/topography.mat') good = find(ngood>0 & TopoInterp < -1000); % Zonal integral: pin = utaux_bar + vtauy_bar; pin(ngood<52) = NaN; pin_dx = pin .* dx; pin_z = NaN*ones(length(NbLatitudes),1); for n = 1:length(NbLatitudes) gz = find(~isnan(pin_dx(n,:))); pin_z(n) = sum(pin_dx(n,gz))*mperdeg; end