% ncdisp('merra2_at_mcmurdo_19800101-20200131.nc') % Global Attributes: % Description = 'MERRA-2 profiles of T, U, V, P, Vortex Marker, Distance from Vortex Edge at McMurdo Station (-77.85, 166.67). Daily at 0Z interpolated to 1 km altitude levels from 20-65 km.' % Author = 'V. Lynn Harvey, file was created using interp_merra22point_theta_mcmurdo.pro' % nz = 'Number of Altitude Levels' % kday = 'Number of Days' % Dimensions: % kday = 14641 % nz = 46 % Variables: % DATES % Size: 14641x1 % Dimensions: kday % Datatype: int32 % ALTITUDE % Size: 46x1 % Dimensions: nz % Datatype: single % U % Size: 14641x46 % Dimensions: kday,nz % Datatype: single % Attributes: % long_name = 'ZONAL WIND' % units = 'M/S' % V % Size: 14641x46 % Dimensions: kday,nz % Datatype: single % Attributes: % long_name = 'MERIDIONAL WIND' % units = 'M/S' % P % Size: 14641x46 % Dimensions: kday,nz % Datatype: single % Attributes: % long_name = 'PRESSURE' % units = 'HPA' % T % Size: 14641x46 % Dimensions: kday,nz % Datatype: single % Attributes: % long_name = 'TEMPERATURE' % units = 'K' % MARK % Size: 14641x46 % Dimensions: kday,nz % Datatype: single % Attributes: % long_name = 'VORTEX MARKER' % units = 'dimensionless' % EDGE_DISTANCE % Size: 14641x46 % Dimensions: kday,nz % Datatype: single % Attributes: % long_name = 'EDGE DISTANCE' % units = 'KM' % VORTEX_EDGE_WIND % Size: 14641x46 % Dimensions: kday,nz % Datatype: single % Attributes: % long_name = 'MEAN VORTEX EDGE WIND SPEED' % units = 'M/S' clc; close all; clear; mainpath = 'C:\Users\l1989\Desktop\PolarVortex\FigureofPolarVortex\';cd(mainpath); %Reading data % read new version data if ~exist('PolarVortexEdge', 'var') ttt = struct('DATE',{},'ALTITUDE',{},'EDGE_DISTANCE',{},'VORTEX_EDGE_WIND',{},'U',{},'V',{},'P',{},'T',{}); PolarVortexEdge = struct2cell(ttt); PolarVortexEdge.DATE = datetime(num2str(ncread('merra2_at_mcmurdo_19800101-20200131.nc','DATES')),'InputFormat','yyyyMMdd'); PolarVortexEdge.ALTITUDE = ncread('merra2_at_mcmurdo_19800101-20200131.nc','ALTITUDE'); PolarVortexEdge.EDGE_DISTANCE = ncread('merra2_at_mcmurdo_19800101-20200131.nc','EDGE_DISTANCE'); PolarVortexEdge.VORTEX_EDGE_WIND = ncread('merra2_at_mcmurdo_19800101-20200131.nc','VORTEX_EDGE_WIND'); PolarVortexEdge.U = ncread('merra2_at_mcmurdo_19800101-20200131.nc','U'); PolarVortexEdge.V = ncread('merra2_at_mcmurdo_19800101-20200131.nc','V'); PolarVortexEdge.P = ncread('merra2_at_mcmurdo_19800101-20200131.nc','P'); PolarVortexEdge.T = ncread('merra2_at_mcmurdo_19800101-20200131.nc','T'); end % read old version data % if ~exist('PolarVortexEdgeOld', 'var') % ttt = struct('DATE',{},'ALTITUDE',{},'EDGE_DISTANCE',{},'VORTEX_EDGE_WIND',{},'U',{},'V',{},'P',{},'T',{}); % PolarVortexEdgeOld = struct2cell(ttt); % PolarVortexEdgeOld.DATE = datetime(num2str(ncread('merra2_at_mcmurdo_20100101-20161231.nc','DATES')),'InputFormat','yyyyMMdd'); % PolarVortexEdgeOld.ALTITUDE = ncread('merra2_at_mcmurdo_20100101-20161231.nc','ALTITUDE'); % PolarVortexEdgeOld.EDGE_DISTANCE = ncread('merra2_at_mcmurdo_20100101-20161231.nc','EDGE_DISTANCE'); % PolarVortexEdgeOld.U = ncread('merra2_at_mcmurdo_20100101-20161231.nc','U'); % PolarVortexEdgeOld.V = ncread('merra2_at_mcmurdo_20100101-20161231.nc','V'); % PolarVortexEdgeOld.P = ncread('merra2_at_mcmurdo_20100101-20161231.nc','P'); % PolarVortexEdgeOld.T = ncread('merra2_at_mcmurdo_20100101-20161231.nc','T'); % end % alititude value, unit: km alt = 45; % time range StartDate = datetime(1999,1,1); EndDate = datetime(2019,12,31); durYear = ceil(years(EndDate-StartDate)); [StartDateInd,~] = find(PolarVortexEdge.DATE == StartDate); [EndDateInd,~] = find(PolarVortexEdge.DATE == EndDate); [~,altInd] = min(abs(PolarVortexEdge.ALTITUDE-alt)); % use for old version % altOld = 40; % [StartDateIndOld,~] = find(PolarVortexEdgeOld.DATE == StartDate); % [EndDateIndOld,~] = find(PolarVortexEdgeOld.DATE == EndDate); % [~,altIndOld] = min(abs(PolarVortexEdgeOld.ALTITUDE-altOld)); plotIndex = 1; % 1 - plot edge distance; 2 - plot edge wind; % x value dn_start = datenum(StartDate); xData = PolarVortexEdge.DATE(StartDateInd:EndDateInd); % y value switch plotIndex case 1 yData = PolarVortexEdge.EDGE_DISTANCE(StartDateInd:EndDateInd,altInd); case 2 yData = PolarVortexEdge.VORTEX_EDGE_WIND(StartDateInd:EndDateInd,altInd); end % use for check old version data % xDataOld = PolarVortexEdgeOld.DATE(StartDateIndOld:EndDateIndOld); % yDataOld = PolarVortexEdgeOld.EDGE_DISTANCE(StartDateIndOld:EndDateIndOld,altIndOld); % load Epm data % load('C:\Users\l1989\Desktop\PolarVortex\Epm\Epm2011-2019final.mat'); % timesort_datetime = datetime(timesort,'convertfrom','datenum'); % timesortTotal = timesort; % n = 1; % l = size(gwpedErr,1); % for nYr = 2011:1:2019 % for nMon = 1:1:12 % index = (timesort_datetime.Year == nYr & timesort_datetime.Month == nMon); % if n == 1 % mean30to50EpmTotal = nanmean(gwped(:,index),1); % ErrorMean30to50EpmTotal = sqrt(sum(gwpedErr(:,index).^2,1))/l; % else % mean30to50EpmTotal = [mean30to50EpmTotal nanmean(gwped(:,index),1)]; % ErrorMean30to50EpmTotal = [ErrorMean30to50EpmTotal sqrt(sum(gwpedErr(:,index).^2,1))/l]; % end % n = n+1; % end % end %Plot 5 years of data scrsz = get(0, 'ScreenSize'); figure('position',[1,300,scrsz(3)*0.95,scrsz(4)*0.6]); subplot('position',[0.05,0.08,0.93,0.8]); % plot(xDataOld,yDataOld/1000,'--r'); switch plotIndex case 1 plot(xData, yData/1000);hold on; plot(xData, zeros(length(xData),1),'--');hold on; ylabel('polar vortex edge distance from McMurdo (10^3 km)','fontsize',16);set(gca,'fontsize',16); case 2 plot(xData, yData);hold on; ylabel('polar vortex edge wind speed (m/s)','fontsize',16);set(gca,'fontsize',16); end [IndJJ,~] = find(PolarVortexEdge.DATE(StartDateInd:EndDateInd).Month == 6); scatter(xData(IndJJ),yData(IndJJ)/1000); ax = gca; ax.XAxis.TickValues = StartDate+years(0:1:durYear)+days(1); ax.XAxis.TickLabelFormat = 'yy'; legend({['at ' num2str(alt) ' km']}); %% monthly mean polar vortex edge MonthlyMeanEdgeDis = []; n = 1; for nYr = PolarVortexEdge.DATE(1).Year:1:PolarVortexEdge.DATE(end).Year for nMon = 1:1:12 Ind = PolarVortexEdge.DATE.Year == nYr & PolarVortexEdge.DATE.Month == nMon; MonthlyMeanEdgeDis(:,n) = nanmean(PolarVortexEdge.EDGE_DISTANCE(Ind,:),1); MonthlyTimeSeries40yrs(n) = datetime(nYr,nMon,15); n = n+1; end end dnMonthlyTimeSeries40yrs = datenum(MonthlyTimeSeries40yrs); %% %calc wind Wind = sqrt(PolarVortexEdge.U.^2+PolarVortexEdge.V.^2); WindAlt = 30; [~,WindAltInd] = min(abs(PolarVortexEdge.ALTITUDE-WindAlt)); figure; plot(datenum(xData)-dn_start, smooth(Wind(StartDateInd:EndDateInd,altInd),14),'LineWidth',2); hold on; xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); ylabel('Wind Speed (m/s)'); set(gca,'xtick',[datenum('20110101','yyyymmdd')-dn_start,datenum('20110201','yyyymmdd')-dn_start,datenum('20110301','yyyymmdd')-dn_start,... datenum('20110401','yyyymmdd')-dn_start,datenum('20110501','yyyymmdd')-dn_start,datenum('20110601','yyyymmdd')-dn_start,... datenum('20110701','yyyymmdd')-dn_start,datenum('20110801','yyyymmdd')-dn_start,datenum('20110901','yyyymmdd')-dn_start,... datenum('20111001','yyyymmdd')-dn_start,datenum('20111101','yyyymmdd')-dn_start,datenum('20111201','yyyymmdd')-dn_start,... datenum('20120101','yyyymmdd')-dn_start,datenum('20120201','yyyymmdd')-dn_start,datenum('20120301','yyyymmdd')-dn_start,... datenum('20120401','yyyymmdd')-dn_start,datenum('20120501','yyyymmdd')-dn_start,datenum('20120601','yyyymmdd')-dn_start,... datenum('20120701','yyyymmdd')-dn_start,datenum('20120801','yyyymmdd')-dn_start,datenum('20120901','yyyymmdd')-dn_start,... datenum('20121001','yyyymmdd')-dn_start,datenum('20121101','yyyymmdd')-dn_start,datenum('20121201','yyyymmdd')-dn_start,... datenum('20130101','yyyymmdd')-dn_start,datenum('20130201','yyyymmdd')-dn_start,datenum('20130301','yyyymmdd')-dn_start,... datenum('20130401','yyyymmdd')-dn_start,datenum('20130501','yyyymmdd')-dn_start,datenum('20130601','yyyymmdd')-dn_start,... datenum('20130701','yyyymmdd')-dn_start,datenum('20130801','yyyymmdd')-dn_start,datenum('20130901','yyyymmdd')-dn_start,... datenum('20131001','yyyymmdd')-dn_start,datenum('20131101','yyyymmdd')-dn_start,datenum('20131201','yyyymmdd')-dn_start,... datenum('20140101','yyyymmdd')-dn_start,datenum('20140201','yyyymmdd')-dn_start,datenum('20140301','yyyymmdd')-dn_start,... datenum('20140401','yyyymmdd')-dn_start,datenum('20140501','yyyymmdd')-dn_start,datenum('20140601','yyyymmdd')-dn_start,... datenum('20140701','yyyymmdd')-dn_start,datenum('20140801','yyyymmdd')-dn_start,datenum('20140901','yyyymmdd')-dn_start,... datenum('20141001','yyyymmdd')-dn_start,datenum('20141101','yyyymmdd')-dn_start,datenum('20141201','yyyymmdd')-dn_start,... datenum('20150101','yyyymmdd')-dn_start,datenum('20150201','yyyymmdd')-dn_start,datenum('20150301','yyyymmdd')-dn_start,... datenum('20150401','yyyymmdd')-dn_start,datenum('20150501','yyyymmdd')-dn_start,datenum('20150601','yyyymmdd')-dn_start,... datenum('20150701','yyyymmdd')-dn_start,datenum('20150801','yyyymmdd')-dn_start,datenum('20150901','yyyymmdd')-dn_start,... datenum('20151001','yyyymmdd')-dn_start,datenum('20151101','yyyymmdd')-dn_start,datenum('20151201','yyyymmdd')-dn_start,... datenum('20160101','yyyymmdd')-dn_start,datenum('20160201','yyyymmdd')-dn_start,datenum('20160301','yyyymmdd')-dn_start,... datenum('20160401','yyyymmdd')-dn_start,datenum('20160501','yyyymmdd')-dn_start,datenum('20160601','yyyymmdd')-dn_start,... datenum('20160701','yyyymmdd')-dn_start,datenum('20160801','yyyymmdd')-dn_start,datenum('20160901','yyyymmdd')-dn_start,... datenum('20161001','yyyymmdd')-dn_start,datenum('20161101','yyyymmdd')-dn_start,datenum('20161201','yyyymmdd')-dn_start,... datenum('20170101','yyyymmdd')-dn_start,datenum('20170201','yyyymmdd')-dn_start,datenum('20170301','yyyymmdd')-dn_start,... datenum('20170401','yyyymmdd')-dn_start,datenum('20170501','yyyymmdd')-dn_start,datenum('20170601','yyyymmdd')-dn_start,... datenum('20170701','yyyymmdd')-dn_start,datenum('20170801','yyyymmdd')-dn_start,datenum('20170901','yyyymmdd')-dn_start,... datenum('20171001','yyyymmdd')-dn_start,datenum('20171101','yyyymmdd')-dn_start,datenum('20171201','yyyymmdd')-dn_start,... datenum('20180101','yyyymmdd')-dn_start,datenum('20180201','yyyymmdd')-dn_start,datenum('20180301','yyyymmdd')-dn_start,... datenum('20180401','yyyymmdd')-dn_start,datenum('20180501','yyyymmdd')-dn_start,datenum('20180601','yyyymmdd')-dn_start,... datenum('20180701','yyyymmdd')-dn_start,datenum('20180801','yyyymmdd')-dn_start,datenum('20180901','yyyymmdd')-dn_start,... datenum('20181001','yyyymmdd')-dn_start,datenum('20181101','yyyymmdd')-dn_start,datenum('20181201','yyyymmdd')-dn_start,... datenum('20190101','yyyymmdd')-dn_start,datenum('20190201','yyyymmdd')-dn_start,datenum('20190301','yyyymmdd')-dn_start,... datenum('20190401','yyyymmdd')-dn_start,datenum('20190501','yyyymmdd')-dn_start,datenum('20190601','yyyymmdd')-dn_start,... datenum('20190701','yyyymmdd')-dn_start,datenum('20190801','yyyymmdd')-dn_start,datenum('20190901','yyyymmdd')-dn_start,... datenum('20191001','yyyymmdd')-dn_start,datenum('20191101','yyyymmdd')-dn_start,datenum('20191201','yyyymmdd')-dn_start,],'xticklabel',... ['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';],'FontSize',15); text(datenum('20110101','yyyymmdd')-dn_start+110,-8,'2011','fontsize',16);text(datenum('20120101','yyyymmdd')-dn_start+110,-8,'2012','fontsize',16); text(datenum('20130101','yyyymmdd')-dn_start+110,-8,'2013','fontsize',16);text(datenum('20140101','yyyymmdd')-dn_start+110,-8,'2014','fontsize',16); text(datenum('20150101','yyyymmdd')-dn_start+110,-8,'2015','fontsize',16);text(datenum('20160101','yyyymmdd')-dn_start+110,-8,'2016','fontsize',16); text(datenum('20170101','yyyymmdd')-dn_start+110,-8,'2017','fontsize',16);text(datenum('20180101','yyyymmdd')-dn_start+110,-8,'2018','fontsize',16); text(datenum('20190101','yyyymmdd')-dn_start+110,-8,'2019','fontsize',16); yyaxis right plot(timesortTotal-dn_start,smooth(mean30to50EpmTotal,7),'LineWidth',2); ylabel('Epm (J/kg)'); legend('30km Wind','Epm'); title('30km wind from Lynns data versus Epm'); %% read QBO index from MERRA2 folderPath = 'F:\New folder\MERRA2\MERRA2_1980-2019_QBO\'; addpath('F:\New folder\MatlabFunc\'); % ncdisp([folderPath 'MERRA2_100.inst6_3d_ana_Np.19800101.SUB.nc']); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Dimensions: % time = 1 (UNLIMITED) % lev = 42 % lat = 5 % lon = 576 % bnds = 2 % Variables: % U % Size: 576x5x42x1 % Dimensions: lon,lat,lev,time % Datatype: single % Attributes: % standard_name = 'eastward_wind' % long_name = 'Eastward wind component' % units = 'm/s' % _FillValue = 999999986991104 % missing_value = 999999986991104 % fmissing_value = 999999986991104 % vmax = 999999986991104 % vmin = -999999986991104 % lat % Size: 5x1 % Dimensions: lat % Datatype: double % Attributes: % standard_name = 'latitude' % long_name = 'latitude' % units = 'degrees_north' % axis = 'Y' % lev % Size: 42x1 % Dimensions: lev % Datatype: double % Attributes: % standard_name = 'air_pressure' % long_name = 'vertical level' % units = 'hPa' % positive = 'down' % axis = 'Z' % lon % Size: 576x1 % Dimensions: lon % Datatype: double % Attributes: % standard_name = 'longitude' % long_name = 'longitude' % units = 'degrees_east' % axis = 'X' % time % Size: 1x1 % Dimensions: time % Datatype: double % Attributes: % standard_name = 'time' % long_name = 'time' % bounds = 'time_bnds' % units = 'minutes since 1980-01-01 00:00:00' % calendar = 'standard' % time_bnds % Size: 2x1 % Dimensions: bnds,time % Datatype: double StartDate = [1980 01]; EndDate = [2019 12]; % MERRA2OutputFileList(StartDate,EndDate,folderPath); lat = ncread([folderPath 'MERRA2_100.inst6_3d_ana_Np.19800101.SUB.nc'],'lat'); lev = ncread([folderPath 'MERRA2_100.inst6_3d_ana_Np.19800101.SUB.nc'],'lev'); lon = ncread([folderPath 'MERRA2_100.inst6_3d_ana_Np.19800101.SUB.nc'],'lon'); ZonalMeanWind = zeros(length(FileList)-2,length(lev)); for nCnt = 1:1:length(FileList)-2 temp = ncread([FileList(nCnt+2).folder '\' FileList(nCnt+2).name],'U'); ZonalMeanWind(nCnt,:) = nanmean(temp(:,3,:),1); end %% plot MERRA2 QBO index - figure 5(b)(c) load('MERRA2QBOIndex.mat'); PressureVal = [20 30]; %unit: hPa [PressureInd,~] = find(lev == PressureVal); yData = mean(ZonalMeanWind(StartDateInd:EndDateInd,PressureInd),2); scrsz = get(0, 'ScreenSize'); figure('position',[1,300,scrsz(3)*0.98,scrsz(4)*0.6]); subplot('position',[0.05,0.59,0.9,0.38]); % subplot(2,1,1); % yyaxis left [IndJJ,~] = find(PolarVortexEdge.DATE(StartDateInd:EndDateInd).Month == 6); % [IndA,~] = find(PolarVortexEdge.DATE(StartDateInd:EndDateInd).Month == 4 &PolarVortexEdge.DATE(StartDateInd:EndDateInd).Day >= 15); % [IndM,~] = find(PolarVortexEdge.DATE(StartDateInd:EndDateInd).Month == 5 &PolarVortexEdge.DATE(StartDateInd:EndDateInd).Day <= 15); % IndJJ = [IndA;IndM]; p1 = plot(xData,yData);hold on; plot(xData,zeros(length(xData),1),'--');hold on; p2 = scatter(xData(IndJJ),yData(IndJJ),'filled'); ylabel('zonal wind at equator (m/s)') ylim([-40 40]); ax = gca; set(ax,'fontsize',16); ax.XAxis.TickValues = StartDate+years(0:1:durYear)+days(2); ax.XAxis.TickLabelFormat = 'yy'; legend({['at ' num2str(PressureVal) ' hPa']}); title('MERRA2 QBO Index'); PressureVal = [40 50]; %unit: hPa [PressureInd,~] = find(lev == PressureVal); yData = mean(ZonalMeanWind(StartDateInd:EndDateInd,PressureInd),2); subplot('position',[0.05,0.14,0.9,0.38]); [IndJJ,~] = find(PolarVortexEdge.DATE(StartDateInd:EndDateInd).Month == 6); % [IndA,~] = find(PolarVortexEdge.DATE(StartDateInd:EndDateInd).Month == 4 &PolarVortexEdge.DATE(StartDateInd:EndDateInd).Day >= 15); % [IndM,~] = find(PolarVortexEdge.DATE(StartDateInd:EndDateInd).Month == 5 &PolarVortexEdge.DATE(StartDateInd:EndDateInd).Day <= 15); % IndJJ = [IndA;IndM]; p1 = plot(xData,yData);hold on; plot(xData,zeros(length(xData),1),'--');hold on; p2 = scatter(xData(IndJJ),yData(IndJJ),'filled'); ylabel('zonal wind at equator (m/s)') ylim([-40 40]); ax = gca; set(ax,'fontsize',16); ax.XAxis.TickValues = StartDate+years(0:1:durYear)+days(2); ax.XAxis.TickLabelFormat = 'yy'; legend({['at ' num2str(PressureVal) ' hPa']}); %% QBO anomaly ?? % [Ind2000,~] = find(PolarVortexEdge.DATE.Year == 2000); % QBOAnomaly = zeros(EndDateInd-StartDateInd+1,1); % AnnualZonalMeanZonalWind = zeros(366,1); % for nDay = 1:1:366 % [Ind,~] = find(PolarVortexEdge.DATE >= StartDate & PolarVortexEdge.DATE <= EndDate & PolarVortexEdge.DATE.Month == PolarVortexEdge.DATE(Ind2000(nDay)).Month & PolarVortexEdge.DATE.Day == PolarVortexEdge.DATE(Ind2000(nDay)).Day); % AnnualZonalMeanZonalWind(nDay) = mean(ZonalMeanWind(Ind,PressureInd),'all'); % QBOAnomaly(Ind) = ZonalMeanWind(Ind,PressureInd)-AnnualZonalMeanZonalWind(nDay); % end %% % figure; % plot(xData,QBOAnomaly(StartDateInd:EndDateInd));hold on; % plot(xData,zeros(length(xData),1));hold on; % [IndJJ,~] = find(PolarVortexEdge.DATE(StartDateInd:EndDateInd).Month == 6); % scatter(xData(IndJJ),QBOAnomaly(StartDateInd+IndJJ)); % %% add polar vortex edge distance - figure 5(a) PressureVal = [20 30 40 50]; %unit: hPa [PressureInd,~] = find(lev == PressureVal); yData = mean(ZonalMeanWind(StartDateInd:EndDateInd,PressureInd),2); scrsz = get(0, 'ScreenSize'); figure('position',[1,300,scrsz(3)*0.98,scrsz(4)*0.6]); subplot('position',[0.05,0.59,0.9,0.38]); % subplot(2,1,1); % yyaxis left [IndJJ,~] = find(PolarVortexEdge.DATE(StartDateInd:EndDateInd).Month == 6); % [IndA,~] = find(PolarVortexEdge.DATE(StartDateInd:EndDateInd).Month == 4 &PolarVortexEdge.DATE(StartDateInd:EndDateInd).Day >= 15); % [IndM,~] = find(PolarVortexEdge.DATE(StartDateInd:EndDateInd).Month == 5 &PolarVortexEdge.DATE(StartDateInd:EndDateInd).Day <= 15); % IndJJ = [IndA;IndM]; p1 = plot(xData,yData);hold on; plot(xData,zeros(length(xData),1),'--');hold on; p2 = scatter(xData(IndJJ),yData(IndJJ),'filled'); ylabel('zonal wind at equator (m/s)') ylim([-40 40]); ax = gca; set(ax,'fontsize',16); ax.XAxis.TickValues = StartDate+years(0:1:durYear)+days(2); ax.XAxis.TickLabelFormat = 'yy'; legend({['at ' num2str(PressureVal) ' hPa']}); title('MERRA2 QBO Index'); yyaxis right yData = PolarVortexEdge.EDGE_DISTANCE(StartDateInd:EndDateInd,altInd); IndJJ = PolarVortexEdge.DATE(StartDateInd:EndDateInd).Month ~= 6; p3 = plot(xData, yData/1000);hold on; yData(IndJJ) = NaN; p4 = plot(xData,yData/1000,'-r','LineWidth',3); rx = gca; rx.YAxis(1).Limits = ([-40, 90]); rx.YAxis(2).Limits = ([-10, 2]); set(rx,'fontsize',16); ylabel('polar vortex edge distance from McMurdo (10^3 km)','fontsize',16);xlabel('Year since 2000'); % legend([p1 p2 p3 p4],{['at ' num2str(PressureVal) ' hPa'] 'June Data' ['at ' num2str(alt) ' km'] 'June 15th'});legend('boxoff'); title(['QBO Index (' num2str(PressureVal) ') hPa versus Polar Vortex Edge Distance at ' num2str(alt) 'km from McMurdo']); %% mid winter QBO index and vortex distance corr - figure 5(d) (e) (f) PressureVal = [20 30 40 50]; %unit: hPa [PressureInd,~] = find(lev == PressureVal); n = 1; for nYr = StartDate.Year:1:EndDate.Year IndexJun = find(PolarVortexEdge.DATE.Year == nYr & PolarVortexEdge.DATE.Month == 6); JunMeanQBO(n) = mean(ZonalMeanWind(IndexJun,PressureInd),'all'); JunMeanDis(n) = nanmean(PolarVortexEdge.EDGE_DISTANCE(IndexJun,altInd)); IndexJul = find(PolarVortexEdge.DATE.Year == nYr & PolarVortexEdge.DATE.Month == 7); JulMeanQBO(n) = mean(ZonalMeanWind(IndexJul,PressureInd),'all'); JulMeanDis(n) = nanmean(PolarVortexEdge.EDGE_DISTANCE(IndexJul,altInd)); IndexWin = [IndexJun;IndexJul]; WinMeanQBO(n) = mean(ZonalMeanWind(IndexWin,PressureInd),'all'); WinMeanDis(n) = nanmean(PolarVortexEdge.EDGE_DISTANCE(IndexWin,altInd)); Index1 = find(PolarVortexEdge.DATE.Year == nYr & PolarVortexEdge.DATE.Month == 6 & PolarVortexEdge.DATE.Day >=5); Index2 = find(PolarVortexEdge.DATE.Year == nYr & PolarVortexEdge.DATE.Month == 7 & PolarVortexEdge.DATE.Day <=5); Index = [Index1;Index2]; MeanQBO(n) = mean(ZonalMeanWind(Index,PressureInd),'all'); MeanDis(n) = nanmean(PolarVortexEdge.EDGE_DISTANCE(Index,altInd)); timesort(n) = datetime(nYr,6,30); n = n+1; end scrsz = get(0, 'ScreenSize'); figure('position',[1,300,scrsz(3)*0.35,scrsz(4)*0.5]); % subplot(2,2,1); [R1,P1] = corrcoef(JunMeanQBO,JunMeanDis/1000); myFit = NonLinearModel.fit(JunMeanQBO,JunMeanDis/1000, 'y ~ b0 + b1*x', [0,1]); p = polyfit(JunMeanQBO,JunMeanDis/1000,1); f = polyval(p,JunMeanQBO); plot(JunMeanQBO,JunMeanDis/1000,'o',JunMeanQBO,f,'-','LineWidth',2); a = get(gca,'XTickLabel'); set(gca,'xticklabel',a,'fontsize',16); ylim([-5 -2]);xlim([-26 14]); xlabel('QBO Index (m/s)','fontsize',16);ylabel('Edge Distance (10^3 km)','fontsize',15); set(gca,'box','on');set(gca,'fontsize',18); title({['QBO @ ' num2str(PressureVal) ' hPa vs. Edge Distance @ ' num2str(alt) ' km '];['Correlation Coeff = ' num2str(R1(1,2)) ' P = ' num2str(1-P1(1,2))]},'fontsize',18); % subplot(2,2,2); % [R2,P2] = corrcoef(JulMeanQBO,JulMeanDis/1000); % p = polyfit(JulMeanQBO,JulMeanDis/1000,1); % f = polyval(p,JulMeanQBO); % plot(JulMeanQBO,JulMeanDis/1000,'o',JulMeanQBO,f,'-'); % xlabel('QBO Index (m/s)');ylabel('Edge distance (10^3 km)'); % set(gca,'box','on'); % title({['2000-2019 Jul QBO@' num2str(PressureVal) 'hPa versus Edge Distance@' num2str(alt) 'km '];['CorrCoef = ' num2str(R2(1,2)) ' ConfLev = ' num2str(1-P2(1,2))]}); % % subplot(2,2,3); % [R3,P3] = corrcoef(MeanQBO,MeanDis/1000); % p = polyfit(MeanQBO,MeanDis/1000,1); % f = polyval(p,MeanQBO); % plot(MeanQBO,MeanDis/1000,'o',MeanQBO,f,'-'); % xlabel('QBO Index (m/s)');ylabel('Edge distance (10^3 km)'); % set(gca,'box','on'); % title({['2000-2019 0605-0705 QBO@' num2str(PressureVal) 'hPa versus Edge Distance@' num2str(alt) 'km '];['CorrCoef = ' num2str(R3(1,2)) ' ConfLev = ' num2str(1-P3(1,2))]}); % % subplot(2,2,4); % [R4,P4] = corrcoef(WinMeanQBO,WinMeanDis/1000); % p = polyfit(WinMeanQBO,WinMeanDis/1000,1); % f = polyval(p,WinMeanQBO); % plot(WinMeanQBO,MeanDis/1000,'o',WinMeanQBO,f,'-'); % xlabel('QBO Index (m/s)');ylabel('Edge distance (10^3 km)'); % set(gca,'box','on'); % title({['2000-2019 Mid Winter QBO@' num2str(PressureVal) 'hPa versus Edge Distance@' num2str(alt) 'km '];['CorrCoef = ' num2str(R4(1,2)) ' ConfLev = ' num2str(1-P4(1,2))]}); %% seprate Eastly and westly phase % typical years between 1980 to 1999: e = [1982 1984 1987 1989 1992 1994 1996 1998]; w = [1980 1985 1990 1995 1997 1999]; % overall years between 1980 to 1999: e = [1981 1982 1984 1986 1987 1988 1989 1991 1992 1994 1996 1998]; %w = [1980 1983 1985 1990 1993 1995 1997 1999]; % typical years between 2000 to 2019: e = [2001 2003 2005 2007 2010 2012 2018]; w = [2002 2004 2006 2008 2013 2019]; % overall years between 2000 to 2019: e = [2000 2001 2003 2005 2007 2009 2010 2012 2014 2015 2017 2018]; %w = [2002 2004 2006 2008 2011 2013 2016 2019]; % typical years between 1980 to 2019: e = [1982 1984 1987 1989 1992 1994 1996 1998 2001 2003 2005 2007 2010 2012 2018]; %w = [1980 1985 1990 1995 1997 1999 2002 2004 2006 2008 2013 2019]; % overall years between 1980 to 2019: e = [1981 1982 1984 1986 1987 1988 1989 1991 1992 1994 1996 1998 2000 2001 2003 2005 2007 2009 2010 2012 2014 2015 2017 2018]; %w = [1980 1983 1985 1990 1993 1995 1997 1999 2002 2004 2006 2008 2011 2013 2016 2019]; % All QBO easterly years = [2000 2001 2003 2005 2007 2009 2010 2012 2014 2015 2017 2018]? % All QBO westerly years = [1999 2002 2004 2006 2008 2011 2013 2016 2019]? %Stable QBO easterly years = [2001 2003 2005 2007 2010 2012 2018]; %Stable QBO westerly years = [1999 2002 2004 2006 2008 2013 2019]; EasterlyYear = [2000 2001 2003 2005 2007 2009 2010 2012 2014 2015 2017 2018]; WesterlyYear = [1999 2002 2004 2006 2008 2011 2013 2016 2019]; EasterlyYearEdgeDis = zeros(365,1); EasterlyYearEdgeWind = zeros(365,1); WesterlyYearEdgeDis = zeros(365,1); WesterlyYearEdgeWind = zeros(365,1); EasterlyEdgeDisStdErr = zeros(365,1); EasterlyEdgeWindStdErr = zeros(365,1); WesterlyEdgeDisStdErr = zeros(365,1); WesterlyEdgeWindStdErr = zeros(365,1); % for nYr = StartDate.Year:1:EndDate.Year % IndTmp = find(EastlyYear == nYr); % if isnan(IndTmp) % WestlyYearEdgeDis(nYr,:) = % else [Ind2011,~] = find(PolarVortexEdge.DATE.Year == 2011); for nDay = 1:1:365 [eastInd,~] = find(PolarVortexEdge.DATE.Year == EasterlyYear & PolarVortexEdge.DATE.Month == PolarVortexEdge.DATE(Ind2011(nDay)).Month & PolarVortexEdge.DATE.Day == PolarVortexEdge.DATE(Ind2011(nDay)).Day); EasterlyYearEdgeDis(nDay) = mean(PolarVortexEdge.EDGE_DISTANCE(eastInd,altInd),'all'); EasterlyYearEdgeWind(nDay) = mean(PolarVortexEdge.VORTEX_EDGE_WIND(eastInd,altInd),'all'); Num = size(PolarVortexEdge.EDGE_DISTANCE(eastInd,altInd),1)*size(PolarVortexEdge.EDGE_DISTANCE(eastInd,altInd),2); EasterlyEdgeDisStdErr(nDay) = std(PolarVortexEdge.EDGE_DISTANCE(eastInd,altInd),0,'all','omitnan')/sqrt(Num); EasterlyEdgeWindStdErr(nDay) = std(PolarVortexEdge.VORTEX_EDGE_WIND(eastInd,altInd),0,'all','omitnan')/sqrt(Num); [westerInd,~] = find(PolarVortexEdge.DATE.Year == WesterlyYear & PolarVortexEdge.DATE.Month == PolarVortexEdge.DATE(Ind2011(nDay)).Month & PolarVortexEdge.DATE.Day == PolarVortexEdge.DATE(Ind2011(nDay)).Day); WesterlyYearEdgeDis(nDay) = mean(PolarVortexEdge.EDGE_DISTANCE(westerInd,altInd),'all'); WesterlyYearEdgeWind(nDay) = mean(PolarVortexEdge.VORTEX_EDGE_WIND(westerInd,altInd),'all'); Num = size(PolarVortexEdge.EDGE_DISTANCE(westerInd,altInd),1)*size(PolarVortexEdge.EDGE_DISTANCE(westerInd,altInd),2); WesterlyEdgeDisStdErr(nDay) = std(PolarVortexEdge.EDGE_DISTANCE(westerInd,altInd),0,'all')/sqrt(Num); WesterlyEdgeWindStdErr(nDay) = std(PolarVortexEdge.VORTEX_EDGE_WIND(westerInd,altInd),0,'all')/sqrt(Num); end %% figure 6 addpath('F:\New folder\MatlabFunc\'); scrsz = get(0, 'ScreenSize'); figure('position',[1,300,scrsz(3)*0.8,scrsz(4)*0.8]); subplot('position',[0.06,0.1,0.4,0.8]); IndEdgeDisE = ~isnan(EasterlyYearEdgeDis); IndEdgeDisW = ~isnan(WesterlyYearEdgeDis); [Ind1,~] = find(IndEdgeDisE == 1); [Ind2,~] = find(IndEdgeDisW == 1); SmthPnts = 7; for nCnt = 1:1:length(EasterlyEdgeDisStdErr) switch nCnt case 1 EasterlyEdgeDisSmthStdErr(nCnt) = sqrt(EasterlyEdgeDisStdErr(1).^2+EasterlyEdgeDisStdErr(2).^2+EasterlyEdgeDisStdErr(3).^2+EasterlyEdgeDisStdErr(4).^2)/4; WesterlyEdgeDisSmthStdErr(nCnt) = sqrt(WesterlyEdgeDisStdErr(1).^2+WesterlyEdgeDisStdErr(2).^2+WesterlyEdgeDisStdErr(3).^2+WesterlyEdgeDisStdErr(4).^2)/4; EasterlyEdgeWindSmthStdErr(nCnt) = sqrt(EasterlyEdgeWindStdErr(1).^2+EasterlyEdgeWindStdErr(2).^2+EasterlyEdgeWindStdErr(3).^2+EasterlyEdgeWindStdErr(4).^2)/4; WesterlyEdgeWindSmthStdErr(nCnt) = sqrt(WesterlyEdgeWindStdErr(1).^2+WesterlyEdgeWindStdErr(2).^2+WesterlyEdgeWindStdErr(3).^2+WesterlyEdgeWindStdErr(4).^2)/4; case 2 EasterlyEdgeDisSmthStdErr(nCnt) = sqrt(EasterlyEdgeDisStdErr(1).^2+EasterlyEdgeDisStdErr(2).^2+EasterlyEdgeDisStdErr(3).^2+EasterlyEdgeDisStdErr(4).^2+EasterlyEdgeDisStdErr(5).^2)/5; WesterlyEdgeDisSmthStdErr(nCnt) = sqrt(WesterlyEdgeDisStdErr(1).^2+WesterlyEdgeDisStdErr(2).^2+WesterlyEdgeDisStdErr(3).^2+WesterlyEdgeDisStdErr(4).^2+WesterlyEdgeDisStdErr(5).^2)/5; EasterlyEdgeWindSmthStdErr(nCnt) = sqrt(EasterlyEdgeWindStdErr(1).^2+EasterlyEdgeWindStdErr(2).^2+EasterlyEdgeWindStdErr(3).^2+EasterlyEdgeWindStdErr(4).^2+EasterlyEdgeWindStdErr(5).^2)/5; WesterlyEdgeWindSmthStdErr(nCnt) = sqrt(WesterlyEdgeWindStdErr(1).^2+WesterlyEdgeWindStdErr(2).^2+WesterlyEdgeWindStdErr(3).^2+WesterlyEdgeWindStdErr(4).^2+WesterlyEdgeWindStdErr(5).^2)/5; case 3 EasterlyEdgeDisSmthStdErr(nCnt) = sqrt(EasterlyEdgeDisStdErr(1).^2+EasterlyEdgeDisStdErr(2).^2+EasterlyEdgeDisStdErr(3).^2+EasterlyEdgeDisStdErr(4).^2+EasterlyEdgeDisStdErr(5).^2+EasterlyEdgeDisStdErr(6).^2)/6; WesterlyEdgeDisSmthStdErr(nCnt) = sqrt(WesterlyEdgeDisStdErr(1).^2+WesterlyEdgeDisStdErr(2).^2+WesterlyEdgeDisStdErr(3).^2+WesterlyEdgeDisStdErr(4).^2+WesterlyEdgeDisStdErr(5).^2+WesterlyEdgeDisStdErr(6).^2)/6; EasterlyEdgeWindSmthStdErr(nCnt) = sqrt(EasterlyEdgeWindStdErr(1).^2+EasterlyEdgeWindStdErr(2).^2+EasterlyEdgeWindStdErr(3).^2+EasterlyEdgeWindStdErr(4).^2+EasterlyEdgeWindStdErr(5).^2+EasterlyEdgeWindStdErr(6).^2)/6; WesterlyEdgeWindSmthStdErr(nCnt) = sqrt(WesterlyEdgeWindStdErr(1).^2+WesterlyEdgeWindStdErr(2).^2+WesterlyEdgeWindStdErr(3).^2+WesterlyEdgeWindStdErr(4).^2+WesterlyEdgeWindStdErr(5).^2+WesterlyEdgeWindStdErr(6).^2)/6; case 365 EasterlyEdgeDisSmthStdErr(nCnt) = sqrt(EasterlyEdgeDisStdErr(362).^2+EasterlyEdgeDisStdErr(363).^2+EasterlyEdgeDisStdErr(364).^2+EasterlyEdgeDisStdErr(365).^2)/4; WesterlyEdgeDisSmthStdErr(nCnt) = sqrt(WesterlyEdgeDisStdErr(362).^2+WesterlyEdgeDisStdErr(363).^2+WesterlyEdgeDisStdErr(364).^2+WesterlyEdgeDisStdErr(365).^2)/4; EasterlyEdgeWindSmthStdErr(nCnt) = sqrt(EasterlyEdgeWindStdErr(362).^2+EasterlyEdgeWindStdErr(363).^2+EasterlyEdgeWindStdErr(364).^2+EasterlyEdgeWindStdErr(365).^2)/4; WesterlyEdgeWindSmthStdErr(nCnt) = sqrt(WesterlyEdgeWindStdErr(362).^2+WesterlyEdgeWindStdErr(363).^2+WesterlyEdgeWindStdErr(364).^2+WesterlyEdgeWindStdErr(365).^2)/4; case 364 EasterlyEdgeDisSmthStdErr(nCnt) = sqrt(EasterlyEdgeDisStdErr(361).^2+EasterlyEdgeDisStdErr(362).^2+EasterlyEdgeDisStdErr(363).^2+EasterlyEdgeDisStdErr(364).^2+EasterlyEdgeDisStdErr(365).^2)/5; WesterlyEdgeDisSmthStdErr(nCnt) = sqrt(WesterlyEdgeDisStdErr(361).^2+WesterlyEdgeDisStdErr(362).^2+WesterlyEdgeDisStdErr(363).^2+WesterlyEdgeDisStdErr(364).^2+WesterlyEdgeDisStdErr(365).^2)/5; EasterlyEdgeWindSmthStdErr(nCnt) = sqrt(EasterlyEdgeWindStdErr(361).^2+EasterlyEdgeWindStdErr(362).^2+EasterlyEdgeWindStdErr(363).^2+EasterlyEdgeWindStdErr(364).^2+EasterlyEdgeWindStdErr(365).^2)/5; WesterlyEdgeWindSmthStdErr(nCnt) = sqrt(WesterlyEdgeWindStdErr(361).^2+WesterlyEdgeWindStdErr(362).^2+WesterlyEdgeWindStdErr(363).^2+WesterlyEdgeWindStdErr(364).^2+WesterlyEdgeWindStdErr(365).^2)/5; case 363 EasterlyEdgeDisSmthStdErr(nCnt) = sqrt(EasterlyEdgeDisStdErr(360).^2+EasterlyEdgeDisStdErr(361).^2+EasterlyEdgeDisStdErr(362).^2+EasterlyEdgeDisStdErr(363).^2+EasterlyEdgeDisStdErr(364).^2+EasterlyEdgeDisStdErr(365).^2)/6; WesterlyEdgeDisSmthStdErr(nCnt) = sqrt(WesterlyEdgeDisStdErr(360).^2+WesterlyEdgeDisStdErr(361).^2+WesterlyEdgeDisStdErr(362).^2+WesterlyEdgeDisStdErr(363).^2+WesterlyEdgeDisStdErr(364).^2+WesterlyEdgeDisStdErr(365).^2)/6; EasterlyEdgeWindSmthStdErr(nCnt) = sqrt(EasterlyEdgeWindStdErr(360).^2+EasterlyEdgeWindStdErr(361).^2+EasterlyEdgeWindStdErr(362).^2+EasterlyEdgeWindStdErr(363).^2+EasterlyEdgeWindStdErr(364).^2+EasterlyEdgeWindStdErr(365).^2)/6; WesterlyEdgeWindSmthStdErr(nCnt) = sqrt(WesterlyEdgeWindStdErr(360).^2+WesterlyEdgeWindStdErr(361).^2+WesterlyEdgeWindStdErr(362).^2+WesterlyEdgeWindStdErr(363).^2+WesterlyEdgeWindStdErr(364).^2+WesterlyEdgeWindStdErr(365).^2)/6; otherwise EasterlyEdgeDisSmthStdErr(nCnt) = sqrt(EasterlyEdgeDisStdErr(nCnt-3).^2+EasterlyEdgeDisStdErr(nCnt-2).^2+EasterlyEdgeDisStdErr(nCnt-1).^2+EasterlyEdgeDisStdErr(nCnt).^2+EasterlyEdgeDisStdErr(nCnt+1).^2+EasterlyEdgeDisStdErr(nCnt+2).^2+EasterlyEdgeDisStdErr(nCnt+3).^2)/7; WesterlyEdgeDisSmthStdErr(nCnt) = sqrt(WesterlyEdgeDisStdErr(nCnt-3).^2+WesterlyEdgeDisStdErr(nCnt-2).^2+WesterlyEdgeDisStdErr(nCnt-1).^2+WesterlyEdgeDisStdErr(nCnt).^2+WesterlyEdgeDisStdErr(nCnt+1).^2+WesterlyEdgeDisStdErr(nCnt+2).^2+WesterlyEdgeDisStdErr(nCnt+3).^2)/7; EasterlyEdgeWindSmthStdErr(nCnt) = sqrt(EasterlyEdgeWindStdErr(nCnt-3).^2+EasterlyEdgeWindStdErr(nCnt-2).^2+EasterlyEdgeWindStdErr(nCnt-1).^2+EasterlyEdgeWindStdErr(nCnt).^2+EasterlyEdgeWindStdErr(nCnt+1).^2+EasterlyEdgeWindStdErr(nCnt+2).^2+EasterlyEdgeWindStdErr(nCnt+3).^2)/7; WesterlyEdgeWindSmthStdErr(nCnt) = sqrt(WesterlyEdgeWindStdErr(nCnt-3).^2+WesterlyEdgeWindStdErr(nCnt-2).^2+WesterlyEdgeWindStdErr(nCnt-1).^2+WesterlyEdgeWindStdErr(nCnt).^2+WesterlyEdgeWindStdErr(nCnt+1).^2+WesterlyEdgeWindStdErr(nCnt+2).^2+WesterlyEdgeWindStdErr(nCnt+3).^2)/7; end end plotMethod = 2; switch plotMethod case 1 p1 = plot(Ind1,smooth(EasterlyYearEdgeDis(IndEdgeDisE)/1000,SmthPnts,'rlowess'),'LineWidth',2);hold on; p3 = errorbar(Ind1,smooth(EasterlyYearEdgeDis(IndEdgeDisE)/1000,SmthPnts,'rlowess'),EasterlyEdgeDisStdErr(IndEdgeDisE)/1000,'+b');hold on; p2 = plot(Ind2,smooth(WesterlyYearEdgeDis(IndEdgeDisW)/1000,SmthPnts,'rlowess'),'LineWidth',2);hold on; p4 = errorbar(Ind2,smooth(WesterlyYearEdgeDis(IndEdgeDisW)/1000,SmthPnts,'rlowess'),WesterlyEdgeDisStdErr(IndEdgeDisW)/1000,'+r'); case 2 p1 = shadedErrorBar(Ind1,smooth(EasterlyYearEdgeDis(IndEdgeDisE)/1000,SmthPnts,'rlowess'),smooth(EasterlyEdgeDisStdErr(IndEdgeDisE)/1000,SmthPnts,'rlowess'),'lineprops','b');hold on; p2 = shadedErrorBar(Ind2,smooth(WesterlyYearEdgeDis(IndEdgeDisW)/1000,SmthPnts,'rlowess'),smooth(WesterlyEdgeDisStdErr(IndEdgeDisW)/1000,SmthPnts,'rlowess'),'lineprops','r'); case 3 p1 = plot(Ind1,smooth(EasterlyYearEdgeDis(IndEdgeDisE)/1000,SmthPnts,'rlowess'),'LineWidth',1,'color','b');hold on; p3 = plot(Ind1,smooth(EasterlyYearEdgeDis(IndEdgeDisE)/1000,SmthPnts,'rlowess')+(EasterlyEdgeDisStdErr(IndEdgeDisE)/1000/SmthPnts),'LineWidth',0.5,'color','b');hold on; p2 = plot(Ind2,smooth(WesterlyYearEdgeDis(IndEdgeDisW)/1000,SmthPnts,'rlowess'),'LineWidth',2);hold on; p4 = errorbar(Ind2,smooth(WesterlyYearEdgeDis(IndEdgeDisW)/1000,SmthPnts,'rlowess'),WesterlyEdgeDisStdErr(IndEdgeDisW)/1000/SmthPnts,'+r'); end xlim([1 365]); xticks([1 32 60 91 121 152 182 213 244 274 305 335]);xticklabels(['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D']); set(gca,'fontsize',18);set(gca,'box','on'); ylabel('Polar Vortex Edge Distance from McMurdo (10^3 km)');xlabel('Month'); switch plotMethod case 1 legend([p1 p3 p2 p4],{'QBO Easterly Mean','Uncertainy in Easterly','QBO Westerly Mean','Uncertainy in Westerly'}); case 2 legend([p1.mainLine p2.mainLine],{'QBO Easterly Mean','QBO Westerly Mean'}); end title(['Vortex Edge from McMurdo @ ' num2str(alt) 'km']); subplot('position',[0.56,0.1,0.4,0.8]); IndEdgeWndE = ~isnan(EasterlyYearEdgeWind); IndEdgeWndW = ~isnan(WesterlyYearEdgeWind); [Ind1,~] = find(IndEdgeWndE == 1); [Ind2,~] = find(IndEdgeWndW == 1); switch plotMethod case 1 p1 = plot(Ind1,smooth(EasterlyYearEdgeWind(IndEdgeWndE),SmthPnts,'rlowess'),'LineWidth',2);hold on; p3 = errorbar(Ind1,smooth(EasterlyYearEdgeWind(IndEdgeWndE),SmthPnts,'rlowess'),EasterlyEdgeWindSmthStdErr(IndEdgeWndE),'+b');hold on; p2 = plot(Ind2,smooth(WesterlyYearEdgeWind(IndEdgeWndW),SmthPnts,'rlowess'),'LineWidth',2);hold on; p4 = errorbar(Ind2,smooth(WesterlyYearEdgeWind(IndEdgeWndW),SmthPnts,'rlowess'),WesterlyEdgeWindSmthStdErr(IndEdgeWndW),'+r'); case 2 p1 = shadedErrorBar(Ind1,smooth(EasterlyYearEdgeWind(IndEdgeWndE),SmthPnts,'rlowess'),smooth(EasterlyEdgeWindStdErr(IndEdgeWndE),SmthPnts,'rlowess'),'lineprops','b');hold on; p2 = shadedErrorBar(Ind2,smooth(WesterlyYearEdgeWind(IndEdgeWndW),SmthPnts,'rlowess'),smooth(WesterlyEdgeWindStdErr(IndEdgeWndW),SmthPnts,'rlowess'),'lineprops','r'); case 3 end xlim([1 365]); xticks([1 32 60 91 121 152 182 213 244 274 305 335]);xticklabels(['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D']); set(gca,'fontsize',18);set(gca,'box','on'); ylabel('Polar Vortex Edge Wind (m/s)');xlabel('Month'); switch plotMethod case 1 legend([p1 p3 p2 p4],{'QBO Easterly Mean','Uncertainy in Easterly','QBO Westerly Mean','Uncertainy in Westerly'}); case 2 case 3 end title(['Polar Vortex Edge Wind @ ' num2str(alt) 'km']); %% wind shear as QBO indices PressureLevel = [50 70]; [PresInd,~] = find(lev == PressureLevel); WindShear = ZonalMeanWind(StartDateInd:EndDateInd,PresInd(1))-ZonalMeanWind(StartDateInd:EndDateInd,PresInd(2)); figure; plot(xData,WindShear);hold on; plot(xData,zeros(length(xData),1));hold on; [IndJJ,~] = find(PolarVortexEdge.DATE(StartDateInd:EndDateInd).Month == 6 & PolarVortexEdge.DATE(StartDateInd:EndDateInd).Day <=5); scatter(xData(IndJJ),WindShear(IndJJ),'filled'); %% sun spots Ra = ncread('american_relative_sunspot_number_monthly.nc','Ra'); RaTime = ncread('american_relative_sunspot_number_monthly.nc','time'); % milliseconds since 1970-01-01 RaTime = RaTime/86400000;% transform to days dn_RaStart = datenum(1970,1,1); dnRaTime = RaTime+dn_RaStart; dtRaTime = datetime(dnRaTime,'convertfrom','datenum'); figure; yyaxis left; plot(dtRaTime,Ra); ylabel('Sunspots number'); yyaxis right; % plot(MonthlyTimeSeries40yrs,MonthlyMeanEdgeDis(alt,:)/1000,'-','LineWidth',2); ylabel('Edge Distance (10^3 km)'); % seperate Easterly/Westerly QBO phase % EasterlyPhase = [1981 1982 1983 1984 1986 1987 1988 1989 1991 1992 1994 1996 1998 2000 2001 2003 2005 2007 2009 2010 2012 2014 2017 2018]; % WesterlyPhase = [1980 1985 1990 1993 1995 1997 1999 2002 2004 2006 2008 2011 2013 2015 2016 2019]; SunSpotsThs = 100; EasterlyYear = [2000 2001 2003 2005 2007 2009 2010 2012 2014 2017 2018]; WesterlyYear = [2002 2004 2006 2008 2011 2013 2015 2016 2019]; [IndEasterly,~] = find(dtRaTime.Month == EasterlyYear); [IndWesterly,~] = find(dtRaTime.Month == 7); figure; plot() %% ESNO 3.4 Anomaly vs. QBO Index % load('C:\Users\l1989\Desktop\PolarVortex\Epm\ESNO34Index.mat'); % dnTimeESNO34Ind = datenum(ESNOTime); load('C:\Users\l1989\Desktop\PolarVortex\Epm\NINO(1+2_3_4_3+4)andAnomaly_1950-2019.mat'); QBOe = [1981 1982 1984 1986 1987 1988 1989 1991 1992 1994 1996 1998 2000 2001 2003 2005 2007 2009 2010 2012 2014 2015 2017 2018]; QBOw = [1980 1983 1985 1990 1993 1995 1997 1999 2002 2004 2006 2008 2011 2013 2016 2019]; [IndNINOStart,~] = find(StartDate.Year == NINO34Anomaly.YR & StartDate.Month == NINO34Anomaly.MON); [IndNINOEnd,~] = find(EndDate.Year == NINO34Anomaly.YR & EndDate.Month == NINO34Anomaly.MON); NINOAnomaly = NINO34Anomaly.ANOM34(IndNINOStart:IndNINOEnd); NINOSmthAnomaly = smooth(NINOAnomaly,5); n = 1; for nYr = StartDate.Year:1:EndDate.Year for nMon = 1:1:12 TimeSortMon(n) = datetime(nYr,nMon,15); n = n+1; end end IndJ = find(TimeSortMon.Month == 6); scrsz = get(0, 'ScreenSize'); figure('position',[1,300,scrsz(3),scrsz(4)*0.5]); subplot('position',[0.05,0.15,0.9,0.73]); yyaxis left plot(TimeSortMon,NINOSmthAnomaly,'linewidth',2);hold on; scatter(TimeSortMon(IndJ),NINOSmthAnomaly(IndJ),50,'linewidth',1.5);hold on; plot(TimeSortMon,zeros(1,length(TimeSortMon)),'--');hold on; plot(TimeSortMon,zeros(1,length(TimeSortMon))+0.4,'-.r');hold on; plot(TimeSortMon,zeros(1,length(TimeSortMon))-0.4,'-.r'); ax = gca; set(ax,'fontsize',16); ax.XAxis.TickValues = StartDate+years(0:1:durYear)+days(2); ax.XAxis.TickLabelFormat = 'yy'; ylabel('NINO 3.4 Anomaly (degC)');ylim([-2 8]);xlabel('Year since 1999'); % legend({['at ' num2str(PressureVal) ' hPa']}); PressureVal = [20 30 40 50]; %unit: hPa [PressureInd,~] = find(lev == PressureVal); yData = mean(ZonalMeanWind(StartDateInd:EndDateInd,PressureInd),2); yyaxis right [IndJJ,~] = find(PolarVortexEdge.DATE(StartDateInd:EndDateInd).Month == 6); n = 1;Inde = [];m = 1; for nYr = StartDate.Year:1:EndDate.Year [IndJun,~] = find(PolarVortexEdge.DATE.Month == 6 & PolarVortexEdge.DATE.Year == nYr); QBOJunMonInd(n) = mean(ZonalMeanWind(IndJun,PressureInd),'all'); if ismember(nYr,QBOe) Inde(m) = n; m = m+1; end n = n+1; end p1 = plot(xData,yData);hold on; plot(xData,zeros(length(xData),1),'--');hold on; p2 = scatter(xData(IndJJ),yData(IndJJ),'filled');ylim([-50 20]); ylabel('QBO Index (m/s)'); title(['NINO 3.4 Anomaly vs. QBO Index ' num2str(PressureVal) ' hPa']); % corr [R1 P1] = corrcoef(QBOJunMonInd,NINOSmthAnomaly(IndJ)); figure; temp = NINOSmthAnomaly(IndJ); p1 = scatter(QBOJunMonInd,NINOSmthAnomaly(IndJ),'linewidth',2);hold on; p2 = scatter(QBOJunMonInd(Inde),temp(Inde),'linewidth',2,'markeredgecolor','r'); set(gca,'box','on');xlabel('QBO Index (m/s)'); ylabel('NINO 3.4 Anomaly (degC)'); title({['1980-2019 QBO 20-50 hPa vs. NINO 3.4 Anomaly correlation']; ['Corr. Coeff = ' num2str(R1(1,2)) ' P = ' num2str(1-P1(1,2))]}); legend([p1 p2],{'QBOw','QBOe'}); %% ESNO 3.4 Anomaly vs. Vortex edge distance load('C:\Users\l1989\Desktop\PolarVortex\Epm\NINO(1+2_3_4_3+4)andAnomaly_1950-2019.mat'); QBOe = [1981 1982 1984 1986 1987 1988 1989 1991 1992 1994 1996 1998 2000 2001 2003 2005 2007 2009 2010 2012 2014 2015 2017 2018]; [IndNINOStart,~] = find(StartDate.Year == NINO34Anomaly.YR & StartDate.Month == NINO34Anomaly.MON); [IndNINOEnd,~] = find(EndDate.Year == NINO34Anomaly.YR & EndDate.Month == NINO34Anomaly.MON); NINOAnomaly = NINO34Anomaly.ANOM34(IndNINOStart:IndNINOEnd); NINOSmthAnomaly = smooth(NINOAnomaly,5); n = 1; for nYr = StartDate.Year:1:EndDate.Year for nMon = 1:1:12 TimeSortMon(n) = datetime(nYr,nMon,15); n = n+1; end end IndJ = find(TimeSortMon.Month == 6); scrsz = get(0, 'ScreenSize'); figure('position',[1,300,scrsz(3),scrsz(4)*0.5]); subplot('position',[0.05,0.15,0.9,0.73]); yyaxis left plot(TimeSortMon,NINOSmthAnomaly,'linewidth',2);hold on; scatter(TimeSortMon(IndJ),NINOSmthAnomaly(IndJ),50,'linewidth',1.5);hold on; plot(TimeSortMon,zeros(1,length(TimeSortMon)),'--');hold on; plot(TimeSortMon,zeros(1,length(TimeSortMon))+0.4,'-.r');hold on; plot(TimeSortMon,zeros(1,length(TimeSortMon))-0.4,'-.r'); ax = gca; set(ax,'fontsize',16); ax.XAxis.TickValues = StartDate+years(0:1:durYear)+days(2); ax.XAxis.TickLabelFormat = 'yy'; ylabel('NINO 3.4 Anomaly (degC)');ylim([-2 8]);xlabel(['Year since ' num2str(StartDate.Year)]); yyaxis right yData = PolarVortexEdge.EDGE_DISTANCE(StartDateInd:EndDateInd,altInd); IndJJ = PolarVortexEdge.DATE(StartDateInd:EndDateInd).Month ~= 6; n = 1;Inde=[];m=1; for nYr = StartDate.Year:1:EndDate.Year [IndJun,~] = find(PolarVortexEdge.DATE.Month == 6 & PolarVortexEdge.DATE.Year == nYr); VortexEdgeDisJunMonInd(n) = mean(PolarVortexEdge.EDGE_DISTANCE(IndJun,altInd),'all')/1000; if ismember(nYr,QBOe) Inde(m) = n; m = m+1; end n = n+1; end p3 = plot(xData, yData/1000);hold on; yData(IndJJ) = NaN; p4 = plot(xData,yData/1000,'-r','LineWidth',3);ylim([-10 2]); ylabel('Vortex Edge Distance (10^3 km)'); title(['NINO 3.4 Anomaly vs. Vortex Edge Distance @ ' num2str(alt) ' km']); % corr [R2 P2] = corrcoef(VortexEdgeDisJunMonInd,NINOSmthAnomaly(IndJ)); figure; p1 = scatter(VortexEdgeDisJunMonInd,NINOSmthAnomaly(IndJ),'linewidth',2);hold on; temp2 = NINOSmthAnomaly(IndJ); p2 = scatter(VortexEdgeDisJunMonInd(Inde),temp2(Inde),'linewidth',2,'markeredgecolor','r'); set(gca,'box','on');xlabel('Vortex Edge Distance (10^3 km)'); ylabel('NINO 3.4 Anomaly (degC)'); title({['1980-2019 Vortex Edge Distance @ 45 km vs. NINO 3.4 Anomaly correlation']; ['Corr. Coeff = ' num2str(R2(1,2)) ' P = ' num2str(1-P2(1,2))]}); legend([p1 p2],{'QBOw','QBOe'}); %% ESNO 3.4 Anomaly vs. Epm load('C:\Users\l1989\Desktop\PolarVortex\Epm\NINO(1+2_3_4_3+4)andAnomaly_1950-2019.mat'); [IndNINOStart,~] = find(StartDate.Year == NINO34Anomaly.YR & StartDate.Month == NINO34Anomaly.MON); [IndNINOEnd,~] = find(EndDate.Year == NINO34Anomaly.YR & EndDate.Month == NINO34Anomaly.MON); NINOAnomaly = NINO34Anomaly.ANOM34(IndNINOStart:IndNINOEnd); NINOSmthAnomaly = smooth(NINOAnomaly,1); n = 1; for nYr = StartDate.Year:1:EndDate.Year for nMon = 1:1:12 TimeSortMon(n) = datetime(nYr,nMon,15); n = n+1; end end IndJ = find(TimeSortMon.Month == 6); dnTimeSortMon = datenum(TimeSortMon); scrsz = get(0, 'ScreenSize'); figure('position',[1,300,scrsz(3),scrsz(4)*0.5]); subplot('position',[0.05,0.15,0.9,0.73]); yyaxis left plot(dnTimeSortMon-dn_start,NINOSmthAnomaly,'linewidth',2);hold on; scatter(dnTimeSortMon(IndJ)-dn_start,NINOSmthAnomaly(IndJ),50,'linewidth',1.5);hold on; plot(dnTimeSortMon-dn_start,zeros(1,length(TimeSortMon)),'--');hold on; plot(dnTimeSortMon-dn_start,zeros(1,length(TimeSortMon))+0.4,'-.r');hold on; plot(dnTimeSortMon-dn_start,zeros(1,length(TimeSortMon))-0.4,'-.r'); ylabel('NINO 3.4 Anomaly (degC)'); ylim([-2 3]); % load 2011-2019 Epm load('C:\Users\l1989\Desktop\PolarVortex\Epm\Epm2011-2019final.mat'); timesort_datetime = datetime(timesort,'convertfrom','datenum'); timesortTotal = timesort; n = 1; l = size(gwpedErr,1); for nYr = 2011:1:2019 for nMon = 1:1:12 index = (timesort_datetime.Year == nYr & timesort_datetime.Month == nMon); if n == 1 mean30to50EpmTotal = nanmean(gwped(:,index),1); ErrorMean30to50EpmTotal = sqrt(sum(gwpedErr(:,index).^2,1))/l; else mean30to50EpmTotal = [mean30to50EpmTotal nanmean(gwped(:,index),1)]; ErrorMean30to50EpmTotal = [ErrorMean30to50EpmTotal sqrt(sum(gwpedErr(:,index).^2,1))/l]; end n = n+1; end end dn_start = datenum(2010,12,31); TimeSeries9yrs = datenum(2011,1,1):1:datenum(2019,12,31); dtSeries9yrs = datetime(TimeSeries9yrs,'convertfrom','datenum'); TimeSeries9yrs = TimeSeries9yrs - dn_start; % 9 years fitting ft = cell(9,1); dataMx = cell(9,1); dataMy = cell(9,1); for nYr = 2011:1:2019 index = find(timesort_datetime.Year == nYr); dataMx{nYr-2010} = timesortTotal(index)-dn_start; % temp = day(datetime(dataMx{nYr-2010},'convertfrom','datenum'),'dayofyear'); dataMy{nYr-2010} = mean30to50EpmTotal(index)'; beta0 = [1,1,180,1,180]; ft{nYr-2010} = fitnlm(dataMx{nYr-2010},mean30to50EpmTotal(index)','y ~ b1 + b2.*cos((2*pi/365).*(x1-b3))+b4.*cos((4*pi/365).*(x1-b5))',beta0); i = i+1; end yyaxis right p2 = plot(timesortTotal(1:425)-dn_start,smooth(mean30to50EpmTotal(1:425),7),'-','LineWidth',2);hold on;% Epm is smoothed by 7 points plot(timesortTotal(426:513)-dn_start,smooth(mean30to50EpmTotal(426:513),7),'-','LineWidth',2);hold on; plot(timesortTotal(516:end)-dn_start,smooth(mean30to50EpmTotal(516:end),7),'-','LineWidth',2); NewTimesortTotal = [timesortTotal(1:513);timesort(516:end)]; NewSmoothEpmTotal = [smooth(mean30to50EpmTotal(1:425),7);smooth(mean30to50EpmTotal(426:513),7);smooth(mean30to50EpmTotal(516:end),7)]; NewEpmTotal = [mean30to50EpmTotal(1:425) mean30to50EpmTotal(426:513) mean30to50EpmTotal(516:end)]; ylim([0,14]);xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); set(gca,'YMinorTick','on'); ylabel('E_{pm} (J/kg)','FontSize',30); set(gca,'xtick',[datenum('20110101','yyyymmdd'),datenum('20110201','yyyymmdd'),datenum('20110301','yyyymmdd'),... datenum('20110401','yyyymmdd'),datenum('20110501','yyyymmdd'),datenum('20110601','yyyymmdd'),... datenum('20110701','yyyymmdd'),datenum('20110801','yyyymmdd'),datenum('20110901','yyyymmdd'),... datenum('20111001','yyyymmdd'),datenum('20111101','yyyymmdd'),datenum('20111201','yyyymmdd'),... datenum('20120101','yyyymmdd'),datenum('20120201','yyyymmdd'),datenum('20120301','yyyymmdd'),... datenum('20120401','yyyymmdd'),datenum('20120501','yyyymmdd'),datenum('20120601','yyyymmdd'),... datenum('20120701','yyyymmdd'),datenum('20120801','yyyymmdd'),datenum('20120901','yyyymmdd'),... datenum('20121001','yyyymmdd'),datenum('20121101','yyyymmdd'),datenum('20121201','yyyymmdd'),... datenum('20130101','yyyymmdd'),datenum('20130201','yyyymmdd'),datenum('20130301','yyyymmdd'),... datenum('20130401','yyyymmdd'),datenum('20130501','yyyymmdd'),datenum('20130601','yyyymmdd'),... datenum('20130701','yyyymmdd'),datenum('20130801','yyyymmdd'),datenum('20130901','yyyymmdd'),... datenum('20131001','yyyymmdd'),datenum('20131101','yyyymmdd'),datenum('20131201','yyyymmdd'),... datenum('20140101','yyyymmdd'),datenum('20140201','yyyymmdd'),datenum('20140301','yyyymmdd'),... datenum('20140401','yyyymmdd'),datenum('20140501','yyyymmdd'),datenum('20140601','yyyymmdd'),... datenum('20140701','yyyymmdd'),datenum('20140801','yyyymmdd'),datenum('20140901','yyyymmdd'),... datenum('20141001','yyyymmdd'),datenum('20141101','yyyymmdd'),datenum('20141201','yyyymmdd'),... datenum('20150101','yyyymmdd'),datenum('20150201','yyyymmdd'),datenum('20150301','yyyymmdd'),... datenum('20150401','yyyymmdd'),datenum('20150501','yyyymmdd'),datenum('20150601','yyyymmdd'),... datenum('20150701','yyyymmdd'),datenum('20150801','yyyymmdd'),datenum('20150901','yyyymmdd'),... datenum('20151001','yyyymmdd'),datenum('20151101','yyyymmdd'),datenum('20151201','yyyymmdd'),... datenum('20160101','yyyymmdd'),datenum('20160201','yyyymmdd'),datenum('20160301','yyyymmdd'),... datenum('20160401','yyyymmdd'),datenum('20160501','yyyymmdd'),datenum('20160601','yyyymmdd'),... datenum('20160701','yyyymmdd'),datenum('20160801','yyyymmdd'),datenum('20160901','yyyymmdd'),... datenum('20161001','yyyymmdd'),datenum('20161101','yyyymmdd'),datenum('20161201','yyyymmdd'),... datenum('20170101','yyyymmdd'),datenum('20170201','yyyymmdd'),datenum('20170301','yyyymmdd'),... datenum('20170401','yyyymmdd'),datenum('20170501','yyyymmdd'),datenum('20170601','yyyymmdd'),... datenum('20170701','yyyymmdd'),datenum('20170801','yyyymmdd'),datenum('20170901','yyyymmdd'),... datenum('20171001','yyyymmdd'),datenum('20171101','yyyymmdd'),datenum('20171201','yyyymmdd'),... datenum('20180101','yyyymmdd'),datenum('20180201','yyyymmdd'),datenum('20180301','yyyymmdd'),... datenum('20180401','yyyymmdd'),datenum('20180501','yyyymmdd'),datenum('20180601','yyyymmdd'),... datenum('20180701','yyyymmdd'),datenum('20180801','yyyymmdd'),datenum('20180901','yyyymmdd'),... datenum('20181001','yyyymmdd'),datenum('20181101','yyyymmdd'),datenum('20181201','yyyymmdd'),... datenum('20190101','yyyymmdd'),datenum('20190201','yyyymmdd'),datenum('20190301','yyyymmdd'),... datenum('20190401','yyyymmdd'),datenum('20190501','yyyymmdd'),datenum('20190601','yyyymmdd'),... datenum('20190701','yyyymmdd'),datenum('20190801','yyyymmdd'),datenum('20190901','yyyymmdd'),... datenum('20191001','yyyymmdd'),datenum('20191101','yyyymmdd'),datenum('20191201','yyyymmdd')]-dn_start,'xticklabel',... ['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';],'FontSize',18); yy = -1.8; text(datenum('20110101','yyyymmdd')-dn_start,yy,'2011','fontsize',20);text(datenum('20120101','yyyymmdd')-dn_start,yy,'2012','fontsize',20); text(datenum('20130101','yyyymmdd')-dn_start,yy,'2013','fontsize',20);text(datenum('20140101','yyyymmdd')-dn_start,yy,'2014','fontsize',20); text(datenum('20150101','yyyymmdd')-dn_start,yy,'2015','fontsize',20);text(datenum('20160101','yyyymmdd')-dn_start,yy,'2016','fontsize',20); text(datenum('20170101','yyyymmdd')-dn_start,yy,'2017','fontsize',20);text(datenum('20180101','yyyymmdd')-dn_start,yy,'2018','fontsize',20); text(datenum('20190101','yyyymmdd')-dn_start,yy,'2019','fontsize',20); title('Smooth Epm vs. NINO 3.4 Anomaly','fontsize',20); n = 1;dtNewTimesortTotal = datetime(NewTimesortTotal,'convertfrom','datenum'); for nYr = 2011:1:2019 Indt = find(dtNewTimesortTotal.Month == 6 & dtNewTimesortTotal.Year == nYr); Indtt = find(dtNewTimesortTotal.Month == 7 & dtNewTimesortTotal.Year == nYr); NewSmthMonEpmTotal6(n) = nanmean(NewSmoothEpmTotal(Indt)); NewSmthMonEpmTotal7(n) = nanmean(NewSmoothEpmTotal(Indtt)); if isempty(Indt) NewMaxSmthMonEpmTotal(n) = NaN; else NewMaxSmthMonEpmTotal(n) = max(NewSmoothEpmTotal(Indt)); end n = n+1; end NewSmthMonEpmTotal = [NewSmthMonEpmTotal6 NewSmthMonEpmTotal7]; QBOe = [2012 2014 2015 2017 2018]; QBOw = [2011 2013 2016 2019]; Indnan = ~isnan(NewSmthMonEpmTotal); [~,IndJun] = find(TimeSortMon.Month == 6 & TimeSortMon.Year >= 2011); [~,IndJul] = find(TimeSortMon.Month == 7 & TimeSortMon.Year >= 2011); IndJJNINO = [IndJun IndJul]; n = 1;IndJune=[];IndJule=[];IndJuneEpm=[];IndJuleEpm=[]; for nYr = 2011:1:2019 if nYr == 2012 || nYr == 2014 || nYr == 2015 || nYr == 2017 || nYr == 2018 [~,IndJune(n)] = find(TimeSortMon.Month == 6 & TimeSortMon.Year == nYr); [~,IndJule(n)] = find(TimeSortMon.Month == 7 & TimeSortMon.Year == nYr); [IndJuneEpm,~] = find(dtNewTimesortTotal.Month == 6 & dtNewTimesortTotal.Year == nYr); [IndJuleEpm,~] = find(dtNewTimesortTotal.Month == 7 & dtNewTimesortTotal.Year == nYr); NewSmthMonEpmTotal66(n) = nanmean(NewSmoothEpmTotal(IndJuneEpm)); NewSmthMonEpmTotal77(n) = nanmean(NewSmoothEpmTotal(IndJuleEpm)); n = n+1; end end IndJJQBOe = [IndJune IndJule]; NewSmthMonEpmTotal67 = [NewSmthMonEpmTotal66 NewSmthMonEpmTotal77]; temp = NINOSmthAnomaly(IndJJNINO); % corr [R3 P3] = corrcoef(NewSmthMonEpmTotal(Indnan),temp(Indnan)); figure; p1 = scatter(NewSmthMonEpmTotal(Indnan),temp(Indnan),'linewidth',2);hold on; p2 = scatter(NewSmthMonEpmTotal67,NINOSmthAnomaly(IndJJQBOe),'linewidth',2,'markeredgecolor','r'); set(gca,'box','on');xlabel('Smoothed Epm (J/kg)'); ylabel('NINO 3.4 Anomaly (degC)'); title({'2011-2019 June-July Smoothed Epm vs. NINO 3.4 Anomaly correlation'; ['Corr. Coeff = ' num2str(R3(1,2)) ' P = ' num2str(1-P3(1,2))]}); legend([p1 p2],{'QBOw','QBOe'}); Indnan(5) = 0;Indnan(14) = 0; [R6 P6] = corrcoef(NewSmthMonEpmTotal(Indnan),temp(Indnan)); figure; p1 = scatter(NewSmthMonEpmTotal(Indnan),temp(Indnan),'linewidth',2);hold on; temp1 = NINOSmthAnomaly(IndJJQBOe);Ind2015 = [1,2,4,5,6,7,9,10]; p2 = scatter(NewSmthMonEpmTotal67(Ind2015),temp1(Ind2015),'linewidth',2,'markeredgecolor','r'); set(gca,'box','on');xlabel('Smoothed Epm (J/kg)'); ylabel('NINO 3.4 Anomaly (degC)'); title({'2011-2019 June-July Smoothed Epm vs. NINO 3.4 Anomaly correlation (without 2015)'; ['Corr. Coeff = ' num2str(R6(1,2)) ' P = ' num2str(1-P6(1,2))]}); legend([p1 p2],{'QBOw','QBOe'}); IndJJQBOw = [378,402,438,474,379,403,439,475];EpmQBOw = [NewSmthMonEpmTotal(1) NewSmthMonEpmTotal(3) NewSmthMonEpmTotal(6) NewSmthMonEpmTotal(9) NewSmthMonEpmTotal(10) NewSmthMonEpmTotal(12) NewSmthMonEpmTotal(15) NewSmthMonEpmTotal(18)]; temp11 = NewSmthMonEpmTotal67; Indxx = ~isnan(temp11); [R7 P7] = corrcoef(NewSmthMonEpmTotal67(Indxx),temp1(Indxx)); %% ESNO 3.4 Anomaly vs. Epm anomaly load('C:\Users\l1989\Desktop\PolarVortex\Epm\EpmAnomaly_Monthly_2011-2019.mat'); [~,IndJun] = find(TimeSortMon.Month == 6 & TimeSortMon.Year >= 2011); [~,IndJul] = find(TimeSortMon.Month == 7 & TimeSortMon.Year >= 2011); IndJJNINO = [IndJun IndJul]; [~,IndJun] = find(xxData.Month == 6); [~,IndJul] = find(xxData.Month == 7); IndJJEpmAnom = [IndJun IndJul]; % corr [R4 P4] = corrcoef(absMonEpmAnomaly(IndJJEpmAnom),NINOSmthAnomaly(IndJJNINO)); figure; scatter(absMonEpmAnomaly(IndJJEpmAnom),NINOSmthAnomaly(IndJJNINO)); set(gca,'box','on');xlabel('Smoothed Epm (J/kg)'); ylabel('NINO 3.4 Anomaly (degC)'); title({['2011-2019 June-July Epm anomaly vs. NINO 3.4 Anomaly correlation']; ['Corr. Coeff = ' num2str(R4(1,2)) ' P = ' num2str(1-P4(1,2))]}); %% ESNO 3.4 Anomaly vs. fitting monthly mean Epm Ind1 = [1,3,5,7,9,11,13,15,17,2,4,6,8,10,12,14,16,18];% [R5 P5] = corrcoef(WinterMonMeanEpm(Ind1),NINOSmthAnomaly(IndJJNINO)); figure; scatter(WinterMonMeanEpm(Ind1),NINOSmthAnomaly(IndJJNINO)); set(gca,'box','on');xlabel('Smoothed Epm (J/kg)'); ylabel('NINO 3.4 Anomaly (degC)'); title({['2011-2019 June-July Fitting Epm vs. NINO 3.4 Anomaly correlation']; ['Corr. Coeff = ' num2str(R5(1,2)) ' P = ' num2str(1-P5(1,2))]}); %% generate June/July annual fitting monthly mean Epm n = 1; for nYr = 2011:1:2019 for nMon = 6:7 Ind = find(dtSeries9yrs.Year == nYr & dtSeries9yrs.Month == nMon); dnAll = datenum(dtSeries9yrs(Ind))-dn_start; temp = ft{nYr-2010}.Coefficients.Estimate(1) + ft{nYr-2010}.Coefficients.Estimate(2).*cos((2*pi/365).*(dnAll-ft{nYr-2010}.Coefficients.Estimate(3)))+ft{nYr-2010}.Coefficients.Estimate(4).*cos((4*pi/365).*(dnAll-ft{nYr-2010}.Coefficients.Estimate(5))); WinterMonMeanEpm(n) = mean(temp); WinterMonMeanEpmStdErr(n) = std(temp)/sqrt(length(temp)); n = n+1; end end %% plot double y axies of position and Epm figure; [ax,h1,h2] = plotyy(PolarVortexEdge.DATE(380:end),... PolarVortexEdge.EDGE_DISTANCE(380:end,16),time1, gwped11,... @plot,@plot); set(h1,'Marker','o','markersize',4); set(h2,'Marker','*','markersize',2); legend([h1 h2],'20km','Epm'); line([PolarVortexEdge.DATE(1) PolarVortexEdge.DATE(end)], [0 0],'color','black'); grid on; grid minor; %% seperate years to subplot if ~exist('Index', 'var') Index(1) = 1; Index(2) = find(PolarVortexEdge.DATE == datetime('2011-01-01')); Index(3) = find(PolarVortexEdge.DATE == datetime('2012-01-01')); Index(4) = find(PolarVortexEdge.DATE == datetime('2013-01-01')); Index(5) = find(PolarVortexEdge.DATE == datetime('2014-01-01')); Index(6) = find(PolarVortexEdge.DATE == datetime('2015-01-01')); Index(7) = find(PolarVortexEdge.DATE == datetime('2016-01-01')); Index(8) = length(PolarVortexEdge.DATE); end figure; for i = 1:7 subplot('position',[0.15 0.88-(i-1)*0.144 0.75 0.125]); plot(PolarVortexEdge.DATE(Index(i):(Index(i+1)-1)), ... PolarVortexEdge.EDGE_DISTANCE(Index(i):(Index(i+1)-1),16),'-o','MarkerSize',2); line([PolarVortexEdge.DATE(Index(i)) PolarVortexEdge.DATE((Index(i+1)-1))], [0 0],... 'color','black','linestyle','--'); xlim([PolarVortexEdge.DATE(Index(i)) PolarVortexEdge.DATE((Index(i+1)-1))]); set(gca,'xtick',(PolarVortexEdge.DATE(Index(i))+calmonths(0:1:11))); dateaxis('x',2); text(PolarVortexEdge.DATE(2),2500,'35km','fontsize',14,'color','r'); ylabel('Distance(km)','fontsize',12); hold on; end %% Statistics of Data figure; load('StatisticData.mat'); test1 = table2array(PolarVortexEdge1); for i = 1:10 subplot('position',[0.15 0.905-(i-1)*0.1 0.75 0.08]); y = reshape(test1(11-i,:),[2,7])'; b = bar(y); width = b.BarWidth; for j=1:7 row = y(j, :); % 0.5 is approximate net width of white spacings per group offset = ((width-0.3) / length(row)) / 2; x = linspace(j-offset, j+offset, length(row)); text(x,row,num2str(row'),'vert','bottom','horiz','center'); end %title([num2str(30+(i-1)*5) 'km']); text(0,max(test1(11-i,:))/2,[num2str(65-(i-1)*5) 'km'],'fontsize',14); text(7.5,max(test1(11-i,:)),'blue is fall','fontsize',12); text(7.5,max(test1(11-i,:))-5,'red is spring','fontsize',12); xticklabels({'2010','2011','2012','2013','2014','2015','2016'}); end %% if ~exist('Index', 'var') Index(1) = 1; Index(2) = find(PolarVortexEdge.DATE == datetime('2011-01-01')); Index(3) = find(PolarVortexEdge.DATE == datetime('2012-01-01')); Index(4) = find(PolarVortexEdge.DATE == datetime('2013-01-01')); Index(5) = find(PolarVortexEdge.DATE == datetime('2014-01-01')); Index(6) = find(PolarVortexEdge.DATE == datetime('2015-01-01')); Index(7) = find(PolarVortexEdge.DATE == datetime('2016-01-01')); Index(8) = length(PolarVortexEdge.DATE); end figure; for i = 1:1:6 subplot(2,3,i); [ax,h1,h2] = plotyy(PolarVortexEdge.DATE(380:end),... PolarVortexEdge.EDGE_DISTANCE(380:end,1),time1, gwped11,... @plot,@scatter); set(h1,'Marker','o','markersize',4); legend([h1 h2],'20km','Epm'); line([PolarVortexEdge.DATE(1) PolarVortexEdge.DATE(end)], [0 0],'color','black'); end %% n=1; for nYr = 1980:1:2019 for nMon = 1:1:12 ESNOTime(n) = datetime(nYr,nMon,15); n = n+1; end end