% wind rotation angle from MERRA2 %% read and save MERRA2 data //no need to run after the first time clc;close all;clear; mainpath = 'F:\New folder\MERRA2\'; cd(mainpath); % path = 'MERRA2_20110101-20191231_Diurnal_UV_Temp_ModelLevel_-90to0Latitude\'; % ncdisp([path 'MERRA2_400.inst6_3d_ana_Nv.20110101.SUB.nc']); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Source: % E:\New folder\MERRA2_20110101-20190630_Diurnal_UV_Temp_ModelLevel_aroundMcMurdo\MERRA2_400.inst6_3d_ana_Nv.20110101.SUB.nc % Format: % netcdf4 % Global Attributes: % CDI = 'Climate Data Interface version 1.6.9 (http://mpimet.mpg.de/cdi)' % history = 'Wed Jun 12 05:02:44 2019: cdo -s -L -f nc4 -timmean -sellonlatbox,160.0,180.0,-90.0,-60.0 -selname,T,U,V /tmpdata/regridder/s_1560315737_33632/MERRA2_400.inst6_3d_ana_Nv.20110101.nc4 /tmpdata/regridder/s_1560315737_33632/MERRA2_400.inst6_3d_ana_Nv.20110101.SUB.nc4' % Conventions = 'CF-1' % ks = 40 % nstep = 15760 % ptop = 1 % pint = 15039.3027 % ak = [1 2.00000024 3.27000046 4.75850105 6.60000134 8.93450165 11.9703016 15.9495029 21.134903 27.8526058 36.5041084 47.5806084 61.6779099 79.5134125 101.944023 130.051025 165.079025 208.49704 262.021057 327.643066 407.657104 504.680115 621.680115 761.984192 929.294189 1127.69019 1364.34021 1645.71033 1979.1604 2373.04053 2836.78052 3381.00073 4017.54102 4764.39111 5638.79102 6660.34131 7851.23145 9236.57227 10866.3018 12783.7031 15039.3027 17693.0039 20119.2012 21686.502 22436.3008 22389.8008 21877.5977 21214.998 20325.8984 19309.6953 18161.8965 16960.8965 15625.9961 14290.9951 12869.5938 11895.8623 10918.1709 9936.52148 8909.99219 7883.42188 7062.19824 6436.26367 5805.32129 5169.61084 4533.90088 3898.20093 3257.08081 2609.20068 1961.31055 1313.48035 659.375244 4.80482578 0] % bk = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 8.1754e-09 0.00696 0.02801 0.06372 0.1136 0.15622 0.20035 0.24674 0.2944 0.34338 0.39289 0.44374 0.49459 0.5463 0.58104 0.61582 0.65063 0.6859 0.72117 0.74938 0.77064 0.79195 0.8133 0.83466 0.85602 0.87743 0.89891 0.92039 0.94187 0.96341 0.98495 1] % History = 'Original file generated: Mon Feb 16 04:38:36 2015 GMT' % Comment = 'GMAO filename: d5124_m2_jan10.inst6_3d_ana_Nv.20110101_.nc4' % Filename = 'MERRA2_400.inst6_3d_ana_Nv.20110101.nc4' % Institution = 'NASA Global Modeling and Assimilation Office' % References = 'http://gmao.gsfc.nasa.gov' % Format = 'NetCDF-4/HDF-5' % SpatialCoverage = 'global' % VersionID = '5.12.4' % TemporalRange = '1980-01-01 -> 2016-12-31' % identifier_product_doi_authority = 'http://dx.doi.org/' % ShortName = 'M2I6NVANA' % GranuleID = 'MERRA2_400.inst6_3d_ana_Nv.20110101.nc4' % ProductionDateTime = 'Original file generated: Mon Feb 16 04:38:36 2015 GMT' % LongName = 'MERRA2 inst6_3d_ana_Nv: 3d,6-Hourly,Instantaneous,Model-Level,Analysis,Analyzed Meteorological Fields' % Title = 'MERRA2 inst6_3d_ana_Nv: 3d,6-Hourly,Instantaneous,Model-Level,Analysis,Analyzed Meteorological Fields' % SouthernmostLatitude = '-90.0' % NorthernmostLatitude = '90.0' % WesternmostLongitude = '-180.0' % EasternmostLongitude = '179.375' % LatitudeResolution = '0.5' % LongitudeResolution = '0.625' % DataResolution = '0.5 x 0.625 (72 native layers)' % Contact = 'http://gmao.gsfc.nasa.gov' % identifier_product_doi = '10.5067/IUUF4WB9FT4W' % RangeBeginningDate = '2011-01-01' % RangeBeginningTime = '00:00:00.000000' % RangeEndingDate = '2011-01-01' % RangeEndingTime = '18:00:00.000000' % CDO = 'Climate Data Operators version 1.6.9 (http://mpimet.mpg.de/cdo)' % Dimensions: % time = 1 (UNLIMITED) % lev = 72 % lat = 61 % lon = 33 % bnds = 2 % Variables: % T % Size: 33x61x72x1 % Dimensions: lon,lat,lev,time % Datatype: single % Attributes: % standard_name = 'air_temperature' % long_name = 'Air temperature' % units = 'K' % _FillValue = 999999986991104 % missing_value = 999999986991104 % fmissing_value = 999999986991104 % vmax = 999999986991104 % vmin = -999999986991104 % U % Size: 33x61x72x1 % 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 % V % Size: 33x61x72x1 % Dimensions: lon,lat,lev,time % Datatype: single % Attributes: % standard_name = 'northward_wind' % long_name = 'Northward wind component' % units = 'm/s' % _FillValue = 999999986991104 % missing_value = 999999986991104 % fmissing_value = 999999986991104 % vmax = 999999986991104 % vmin = -999999986991104 % lat % Size: 61x1 % Dimensions: lat % Datatype: double % Attributes: % standard_name = 'latitude' % long_name = 'latitude' % units = 'degrees_north' % axis = 'Y' % lev % Size: 72x1 % Dimensions: lev % Datatype: double % Attributes: % long_name = 'vertical level' % units = 'layer' % positive = 'down' % axis = 'Z' % lon % Size: 33x1 % 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 2011-01-01 00:00:00' % calendar = 'standard' % time_bnds % Size: 2x1 % Dimensions: bnds,time % Datatype: double %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % load altitude(lev),lon,lat a=load('MERRA2_20110101-20191231_Diurnal_UV_Temp_ModelLevel_-90to0Latitude\MERRA2_ModelLevel 72.txt'); alt = struct('ModelLev',a(:,1),'Pressure',a(:,2)); alt.Altitude = atmospalt(alt.Pressure*100); lon = ncread('MERRA2_20110101-20191231_Diurnal_UV_Temp_ModelLevel_-90to0Latitude\MERRA2_400.inst6_3d_ana_Nv.20110101.SUB.nc','lon'); lat = ncread('MERRA2_20110101-20191231_Diurnal_UV_Temp_ModelLevel_-90to0Latitude\MERRA2_400.inst6_3d_ana_Nv.20110101.SUB.nc','lat'); % read nc files - lat: 1-41 is -80 to - 60 lon: 545-576 is 160 to 179.375 % lev: 9-56 is 60 to 3 km lon = lon(545:576); lat = lat(1:41); alt = alt.Altitude(9:56); alt = alt./1000; if ~exist('F:\New folder\MERRA2\UV(-80to60s)(160to179.375)(3to60km)_20110101to20191231.mat','file') DateStart = datetime([2011 01 01]); DateEnd = datetime([2019 12 31]); mon=months(datestr(DateStart),datestr(DateEnd))+1; filename0 = 'MERRA2_20110101-20191231_Diurnal_UV_Temp_ModelLevel_-90to0Latitude\MERRA2_400.inst6_3d_ana_Nv.';filename1 = '.SUB.nc'; nCnt = 1; for i=1:1:mon Days = eomday(year(DateStart),month(DateStart)); for j = 0:1:(Days-1) DateTemp = DateStart+days(j); str = datestr(DateTemp,'yyyymmdd'); Utemp = ncread([filename0 str filename1],'U'); Vtemp = ncread([filename0 str filename1],'V'); if size(Utemp,2) == 161 U{nCnt} = Utemp(545:576,1:41,9:56); V{nCnt} = Vtemp(545:576,1:41,9:56); else U{nCnt} = Utemp(545:576,21:61,9:56); V{nCnt} = Vtemp(545:576,21:61,9:56); end nCnt = nCnt+1; Time(nCnt) = datenum(str,'yyyymmdd'); end DateStart = DateStart + days(Days); disp(i); end timesortU = datetime(Time,'convertfrom','datenum'); save('F:\New folder\MERRA2\UV(-80to60s)(160to179.375)(3to60km)_20110101to20191231.mat','U','V','Time','timesort'); end %% calculate Wind rotation angle = max angle(11 to 30km) - min angle(11 to 30km) %{ Wind rotaion angle is defined as the biggest wind angle difference between 11 to 30km. So first calculate the wind angle for all altitude, then using the maximum angle minus the minimum angle. Smoothing is hiding in the plotting region. %} if ~exist('U','var') && ~exist('V','var') load('F:\New folder\MERRA2\UV(-80to60s)(160to179.375)(3to60km)_20110101to20191231.mat'); % load('F:\New folder\MERRA2\UV(-80to-60s)(-80to-40)(3to30km)_20110101to20191231.mat'); % load('F:\New folder\MERRA2\UV(-80to-60s)(-180to-160)(3to30km)_20110101to20191231.mat'); end % deal with the outlier for nCnt = 1:1:length(U) nanindex=(U{nCnt}>9999 | U{nCnt}<-9999); U{nCnt}(nanindex) = NaN; end %lat: 1-41 is -80 to - 60 lon: 545-576 is 160 to 179.375 lev: 9-56 is 60 to 3 km latRg = 5:7; lonRg = 10:12; altRg = 18:35;%[18 48]; zAlt = floor(alt(altRg(end))):0.5:ceil(alt(altRg(1))); method = 3; % 1 - calc wind rot angel between raw wind at different location % 2 - take mean wind of different location then calc wind rot % angel % 3 - calc wind rot angel using boundaries % 4 - calc WRA contour (lat, time) switch method case 1 %calc wind angle ang = cell(length(U),1); for nCnt = 1:1:length(U) ang{nCnt} = atan2(V{nCnt}, U{nCnt}); end %lev: 26 to 43 are 11 to 30km WinRotAng = zeros(1,length(U)); for nCnt = 1:1:length(U) WinRotAng(nCnt) = (max(ang{nCnt}(lonRg,latRg,altRg),[],'all')-min(ang{nCnt}(lonRg,latRg,altRg),[],'all'))/pi*180; end index=WinRotAng>360; WinRotAng(index)=360; ZonalWind = zeros(1,length(U)); for nCnt = 1:1:length(U) ZonalWind(nCnt) = sqrt(mean(U{nCnt}(lonRg,latRg,altRg),3).^2 + mean(V{nCnt}(lonRg,latRg,altRg),3).^2); end case 2 MeanU = zeros(length(zAlt),length(U));MeanV = zeros(length(zAlt),length(V)); for nCnt = 1:1:length(U) tempU = reshape(nanmean(U{nCnt}(lonRg,latRg,altRg),[1 2]),[length(altRg),1]); MeanU(:,nCnt) = interp1(alt(altRg),tempU,zAlt,'spline')'; tempV = reshape(nanmean(V{nCnt}(lonRg,latRg,altRg),[1 2]),[length(altRg),1]); MeanV(:,nCnt) = interp1(alt(altRg),tempV,zAlt,'spline')'; end ang = atan2(MeanV,MeanU); WinRotAng = (max(ang,[],1)-min(ang,[],1))/pi*180; index=WinRotAng>360; WinRotAng(index)=360; case 3 MeanU = zeros(length(altRg),length(U));MeanV = zeros(length(altRg),length(V)); for nCnt = 1:1:length(U) MeanU(:,nCnt) = nanmean(U{nCnt}(lonRg,latRg,altRg),[1 2]); MeanV(:,nCnt) = nanmean(V{nCnt}(lonRg,latRg,altRg),[1 2]); end ang = atan2(MeanV,MeanU)/pi*180; BoundaryL = zeros(1,length(U));BoundaryR = zeros(1,length(U)); for nCnt = 1:1:length(U) BoundaryL(nCnt) = ang(1,nCnt); BoundaryR(nCnt) = ang(1,nCnt); for n = 2:1:size(ang,1) temp = ang(n,nCnt)-ang(n-1,nCnt); if temp < 0 && temp > -180 %&& ang(n,nCnt) < BoundaryL(nCnt) if ang(n,nCnt) - BoundaryL(nCnt) < 0 && ang(n,nCnt) - BoundaryL(nCnt) > -180 BoundaryL(nCnt) = ang(n,nCnt); elseif ang(n,nCnt) - BoundaryL(nCnt) > 180 BoundaryL(nCnt) = ang(n,nCnt); end elseif temp < -180%&& (BoundaryL(nCnt)<=BoundaryR(nCnt)) if ang(n,nCnt) - BoundaryR(nCnt) > 0 && ang(n,nCnt) - BoundaryR(nCnt) < 180 BoundaryR(nCnt) = ang(n,nCnt); elseif ang(n,nCnt) - BoundaryL(nCnt) < -180 BoundaryR(nCnt) = ang(n,nCnt); end elseif temp > 0 && temp < 180%&& (BoundaryL(nCnt)<=BoundaryR(nCnt)) if ang(n,nCnt) - BoundaryR(nCnt) > 0 && ang(n,nCnt) - BoundaryR(nCnt) < 180 BoundaryR(nCnt) = ang(n,nCnt); elseif ang(n,nCnt) - BoundaryL(nCnt) < -180 BoundaryR(nCnt) = ang(n,nCnt); end elseif temp > 180 %&& (BoundaryL(nCnt)<=BoundaryR(nCnt)) if ang(n,nCnt) - BoundaryL(nCnt) < 0 && ang(n,nCnt) - BoundaryL(nCnt) > -180 BoundaryL(nCnt) = ang(n,nCnt); elseif ang(n,nCnt) - BoundaryL(nCnt) > 180 BoundaryL(nCnt) = ang(n,nCnt); end end end end WinRotAng = BoundaryR-BoundaryL; index=WinRotAng<0; WinRotAng(index)=WinRotAng(index)+360; case 4 % contour for nLev = 2:1:length(lat)-1 latRg = (nLev-1):1:(nLev+1); MeanU = zeros(length(altRg),length(U));MeanV = zeros(length(altRg),length(V)); for nCnt = 1:1:length(U) MeanU(:,nCnt) = nanmean(U{nCnt}(lonRg,latRg,altRg),[1 2]); MeanV(:,nCnt) = nanmean(V{nCnt}(lonRg,latRg,altRg),[1 2]); end ang = atan2(MeanV,MeanU)/pi*180; BoundaryL = zeros(1,length(U));BoundaryR = zeros(1,length(U)); for nCnt = 1:1:length(U) BoundaryL(nCnt) = ang(1,nCnt); BoundaryR(nCnt) = ang(1,nCnt); for n = 2:1:size(ang,1) temp = ang(n,nCnt)-ang(n-1,nCnt); if temp < 0 && temp > -180 %&& ang(n,nCnt) < BoundaryL(nCnt) if ang(n,nCnt) - BoundaryL(nCnt) < 0 && ang(n,nCnt) - BoundaryL(nCnt) > -180 BoundaryL(nCnt) = ang(n,nCnt); elseif ang(n,nCnt) - BoundaryL(nCnt) > 180 BoundaryL(nCnt) = ang(n,nCnt); end elseif temp < -180%&& (BoundaryL(nCnt)<=BoundaryR(nCnt)) if ang(n,nCnt) - BoundaryR(nCnt) > 0 && ang(n,nCnt) - BoundaryR(nCnt) < 180 BoundaryR(nCnt) = ang(n,nCnt); elseif ang(n,nCnt) - BoundaryL(nCnt) < -180 BoundaryR(nCnt) = ang(n,nCnt); end elseif temp > 0 && temp < 180%&& (BoundaryL(nCnt)<=BoundaryR(nCnt)) if ang(n,nCnt) - BoundaryR(nCnt) > 0 && ang(n,nCnt) - BoundaryR(nCnt) < 180 BoundaryR(nCnt) = ang(n,nCnt); elseif ang(n,nCnt) - BoundaryL(nCnt) < -180 BoundaryR(nCnt) = ang(n,nCnt); end elseif temp > 180 %&& (BoundaryL(nCnt)<=BoundaryR(nCnt)) if ang(n,nCnt) - BoundaryL(nCnt) < 0 && ang(n,nCnt) - BoundaryL(nCnt) > -180 BoundaryL(nCnt) = ang(n,nCnt); elseif ang(n,nCnt) - BoundaryL(nCnt) > 180 BoundaryL(nCnt) = ang(n,nCnt); end end end end tempWRA = BoundaryR-BoundaryL; index=tempWRA<0; tempWRA(index)=tempWRA(index)+360; WinRotAng(nLev-1,:) = tempWRA; end end TimeSeries9yrs = datenum(datetime(2011,1,1):days(1):datetime(2019,12,31)); %% 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 %% plot Epm(7 points smooth) vs wind rotaion angle ---figure 9(a)(b)(c)(d) 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 SmthYdata = smooth(WinRotAng,14); p1 = plot(Time-dn_start,SmthYdata,'LineWidth',1.5);hold on;% Wind rotation angle is smoothed by 14 points(days) dtTime = datetime(Time,'convertfrom','datenum'); [~,IndJJ] = find(dtTime.Month == 6 | dtTime.Month == 7); % scatter(Time(IndJJ),SmthYdata(IndJJ),'filled'); ylim([0,360]);yticks([0 60 120 180 240 300 360]);ylabel('Wind Rotation Angle (deg.)','FontSize',30); 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,20]);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',12); text(datenum('20110101','yyyymmdd')-dn_start,-3,'2011','fontsize',20);text(datenum('20120101','yyyymmdd')-dn_start,-3,'2012','fontsize',20); text(datenum('20130101','yyyymmdd')-dn_start,-3,'2013','fontsize',20);text(datenum('20140101','yyyymmdd')-dn_start,-3,'2014','fontsize',20); text(datenum('20150101','yyyymmdd')-dn_start,-3,'2015','fontsize',20);text(datenum('20160101','yyyymmdd')-dn_start,-3,'2016','fontsize',20); text(datenum('20170101','yyyymmdd')-dn_start,-3,'2017','fontsize',20);text(datenum('20180101','yyyymmdd')-dn_start,-3,'2018','fontsize',20); text(datenum('20190101','yyyymmdd')-dn_start,-3,'2019','fontsize',20); % legend('Wind Rotation Angle','E_{pm} (J/kg)','fontsize',20); % monthly mean wind rot angle IndJun = [];IndJul = [];TimeSortJun = [];TimeSortJul = []; MonthlyMeanJunWindRotAng = zeros(9,1);MonthlyMeanJulWindRotAng = zeros(9,1); stdMonthlyMeanJunWindRotAng = zeros(9,1);stdMonthlyMeanJulWindRotAng = zeros(9,1);stdJJMeanWindRotAng = zeros(9,1); stdMonthlyMeanWindRotAng = zeros(108,1); for nYr = 2011:1:2019 [~,IndJun] = find(dtTime.Year == nYr & dtTime.Month == 6); [~,IndJul] = find(dtTime.Year == nYr & dtTime.Month == 7); TimeSortJun(nYr-2010) = datenum(datetime(nYr,6,15)); TimeSortJul(nYr-2010) = datenum(datetime(nYr,7,15)); MonthlyMeanJunWindRotAng(nYr-2010) = mean(SmthYdata(IndJun)); MonthlyMeanJulWindRotAng(nYr-2010) = mean(SmthYdata(IndJul)); stdMonthlyMeanJunWindRotAng(nYr-2010) = std(SmthYdata(IndJun),0,'all'); stdMonthlyMeanJulWindRotAng(nYr-2010) = std(SmthYdata(IndJul),0,'all'); stdJJMeanWindRotAng(nYr-2010) = std([SmthYdata(IndJun);SmthYdata(IndJul)],0,'all'); end n = 1; for nYr = 2011:1:2019 for nMon = 1:1:12 [~,IndMon] = find(dtTime.Year == nYr & dtTime.Month == nMon); stdMonthlyMeanWindRotAng(n) = std(SmthYdata(IndMon),0,'all'); TimesortMon(n) = datenum(datetime(nYr,nMon,15)); n = n+1; end end yyaxis left p3 = scatter(TimeSortJun-dn_start,stdMonthlyMeanJunWindRotAng.*8,80,'m','filled');hold on; % p4 = scatter(TimeSortJul-dn_start,stdMonthlyMeanJulWindRotAng.*8,'k','filled');hold on; % p4 = scatter(TimeSortJun-dn_start,stdJJMeanWindRotAng,'g','filled'); % p5 = plot(TimesortMon-dn_start,stdMonthlyMeanWindRotAng.*8,'-*b'); set(gca,'TickDir','out');set(gca,'TickLength',[0.006 0.0025]);set(gca,'fontsize',16); legend([p1 p2 p3],{'14 points smth Rot Ang'; '7 points smooth Epm'; 'Std of June'});%'Std of Monthly Wind Rot Ang' strlon1 = num2str(lon(lonRg(1))); strlon2 = num2str(lon(lonRg(end))); strlat1 = num2str(lat(latRg(1))); strlat2 = num2str(lat(latRg(end))); stralt1 = num2str(alt(altRg(1))); stralt2 = num2str(alt(altRg(end))); switch method case 1 title(['Wind Rot Angle: lon: ' strlon1 ' - ' strlon2 ' lat: ' strlat1 ' - ' strlat2 ' alt: ' stralt2 ' - ' stralt1 ' km']); case 2 title(['Wind Rot Angle(Mean for diff location): lon: ' strlon1 ' - ' strlon2 ' lat: ' strlat1 ' - ' strlat2 ' alt: ' stralt2 ' - ' stralt1 ' km']); case 3 title(['Wind Rot Angle: lon: ' strlon1 ' - ' strlon2 ' lat: ' strlat1 ' - ' strlat2 ' alt: ' stralt2 ' - ' stralt1 ' km'],'position',[(datenum('20150701','yyyymmdd')-dn_start),375],'FontSize',18); end % Epm anomaly vs. std wind rot ang load('C:\Users\l1989\Desktop\PolarVortex\Epm\EpmAnomaly_Monthly_2011-2019.mat'); Indxx = find(xxData.Month == 6); Indyy = find(xxData.Month == 7); dtNewTimesortTotal = datetime(NewTimesortTotal,'convertfrom','datenum'); n = 1;NewSmthMonEpmTotal = [];NewMaxSmthMonEpmTotal = [];NewMonEpmTotal = []; for nYr = 2011:1:2019 for nMon = 1:1:12 Indt = find(dtNewTimesortTotal.Month == nMon & dtNewTimesortTotal.Year == nYr); NewSmthMonEpmTotal(n) = nanmean(NewSmoothEpmTotal(Indt)); NewMonEpmTotal(n) = nanmean(NewEpmTotal(Indt)); if isempty(Indt) NewMaxSmthMonEpmTotal(n) = NaN; else NewMaxSmthMonEpmTotal(n) = max(NewSmoothEpmTotal(Indt)); end n = n+1; end end [R P] = corrcoef(stdMonthlyMeanJunWindRotAng,NewSmthMonEpmTotal(Indxx),'Rows','complete'); scrsz = get(0, 'ScreenSize'); figure('position',[1,300,scrsz(3)*0.22,scrsz(4)*0.3]); scatter(stdMonthlyMeanJunWindRotAng,NewSmthMonEpmTotal(Indxx),60,'linewidth',2);hold on; myFit1 = NonLinearModel.fit(stdMonthlyMeanJunWindRotAng,NewSmthMonEpmTotal(Indxx), 'y ~ b0 + b1*x', [0,1]); Ind = ~isnan(NewSmthMonEpmTotal(Indxx)); temp = NewSmthMonEpmTotal(Indxx); yData = temp(Ind); pl = polyfit(stdMonthlyMeanJunWindRotAng(Ind)',yData,1); f = polyval(pl,stdMonthlyMeanJunWindRotAng(Ind)'); plot(stdMonthlyMeanJunWindRotAng(Ind)',f,'-','LineWidth',2);hold on; set(gca,'box','on');set(gca,'fontsize',18);xlim([7 29]); ylabel('Monthly Epm (J/kg)','fontsize',14);xlabel('WRA Standard Deviation (deg)','fontsize',14); title(['Corr. Coef = ' num2str(roundn(R(1,2),-2)) ', P = ' num2str(roundn(1-P(1,2),-2))]); %title({'Monthly Epm vs. WRA Standard Deviation'; ['coef = ' num2str(roundn(R(1,2),-2)) ' P = ' num2str(roundn(1-P(1,2),-2))]}); figure('position',[1,300,scrsz(3)*0.22,scrsz(4)*0.3]); [R2 P2] = corrcoef(NewSmthMonEpmTotal(Indxx),MonthlyMeanJunWindRotAng,'Rows','complete'); scatter(MonthlyMeanJunWindRotAng,NewSmthMonEpmTotal(Indxx),60,'linewidth',2);hold on; myFit2 = NonLinearModel.fit(MonthlyMeanJunWindRotAng,NewSmthMonEpmTotal(Indxx), 'y ~ b0 + b1*x', [0,1]); pl = polyfit(MonthlyMeanJunWindRotAng(Ind)',yData,1); f = polyval(pl,MonthlyMeanJunWindRotAng(Ind)'); plot(MonthlyMeanJunWindRotAng(Ind)',f,'-','LineWidth',2);hold on; set(gca,'box','on');set(gca,'fontsize',18);xlim([36 73]); ylabel('Monthly Epm (J/kg)','fontsize',14);xlabel('Wind Rotation Angle (deg)','fontsize',14); title(['Corr. Coef = ' num2str(roundn(R2(1,2),-2)) ', P = ' num2str(roundn(1-P2(1,2),-2))]); %title({'Monthly Epm vs. WRA'; ['coef = ' num2str(roundn(R2(1,2),-2)) ' P = ' num2str(roundn(1-P2(1,2),-2))]}); % polar vortex edge vs. std wind rot ang load('C:\Users\l1989\Desktop\PolarVortex\FigureofPolarVortex\PolarVortexEdge1980-2019.mat'); Alt = 45; [altInd,~] = find(PolarVortexEdge.ALTITUDE == Alt); for nYr = 2011:1:2019 Ind = find(PolarVortexEdge.DATE.Year == nYr & PolarVortexEdge.DATE.Month == 6); Ind2 = find(PolarVortexEdge.DATE.Year == nYr & PolarVortexEdge.DATE.Month == 7); MonthlyEdgeDisJun(nYr-2010) = nanmean(PolarVortexEdge.EDGE_DISTANCE(Ind,altInd)); MonthlyEdgeDisJul(nYr-2010) = nanmean(PolarVortexEdge.EDGE_DISTANCE(Ind2,altInd)); end [R1 P1] = corrcoef(MonthlyEdgeDisJun/1000,stdMonthlyMeanJunWindRotAng); figure('position',[1,300,scrsz(3)*0.22,scrsz(4)*0.3]); scatter(MonthlyEdgeDisJun/1000,stdMonthlyMeanJunWindRotAng,60,'linewidth',2);hold on; myFit3 = NonLinearModel.fit(MonthlyEdgeDisJun/1000,stdMonthlyMeanJunWindRotAng, 'y ~ b0 + b1*x', [0,1]); p = polyfit(MonthlyEdgeDisJun/1000,stdMonthlyMeanJunWindRotAng',1); f = polyval(p,MonthlyEdgeDisJun/1000); plot(MonthlyEdgeDisJun/1000,f,'-','LineWidth',2);hold on; set(gca,'box','on');set(gca,'fontsize',18);xlim([-4.6 -2.8]); xlabel('Vortex Edge Distance (10^3 km)','fontsize',14);ylabel('WRA Standard Deviation (deg)','fontsize',14); title(['Corr. Coeff = ' num2str(roundn(R1(1,2),-2)) ', P = ' num2str(roundn(1-P1(1,2),-2))]); %title({'Vortex Edge Distance vs. WRA Standard Deviation'; ['coef = ' num2str(roundn(R1(1,2),-2)) ' P = ' num2str(roundn(1-P1(1,2),-2))]}); %% contour figure; set(gcf,'outerposition',get(0,'screensize')); [h1, h2]=contourf(Time,lat(2:end-1),WinRotAng,'linestyle','none'); %shading interp; set(h2,'LevelStep',2); handlecolor=colorbar; colormap jet; % axis([start_time_day end_time_day 30 70]); set(gca,'fontsize',16); set(gca,'CLim',[0, 200]); set(gca,'YMinorTick','on');set(gca,'TickDir','out');set(gca,'TickLength',[0.006 0.0025]); set(get(handlecolor,'ylabel'),'String', 'deg');ylabel('latitude'); 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'),],'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'),-80.7,'2011','fontsize',20);text(datenum('20120101','yyyymmdd'),-80.7,'2012','fontsize',20); text(datenum('20130101','yyyymmdd'),-80.7,'2013','fontsize',20);text(datenum('20140101','yyyymmdd'),-80.7,'2014','fontsize',20); text(datenum('20150101','yyyymmdd'),-80.7,'2015','fontsize',20);text(datenum('20160101','yyyymmdd'),-80.7,'2016','fontsize',20); text(datenum('20170101','yyyymmdd'),-80.7,'2017','fontsize',20);text(datenum('20180101','yyyymmdd'),-80.7,'2018','fontsize',20); text(datenum('20190101','yyyymmdd'),-80.7,'2019','fontsize',20); title('WRA between 3 to 30 km'); %% [R P] = corrcoef([absMonEpmAnomaly(Indxx) absMonEpmAnomaly(Indyy)],[stdMonthlyMeanJunWindRotAng;stdMonthlyMeanJulWindRotAng]); figure; scatter([absMonEpmAnomaly(Indxx) absMonEpmAnomaly(Indyy)],[stdMonthlyMeanJunWindRotAng;stdMonthlyMeanJunWindRotAng],80,'linewidth',2);hold on; set(gca,'box','on');set(gca,'fontsize',20); xlabel('Epm Anomaly (J/kg)');ylabel('Wind Rotation Angle Std (deg.)'); title({'June Epm Anomaly vs. Wind Rotation Angle Std'; ['coef = ' num2str(roundn(R(1,2),-2)) ' ConfLev = ' num2str(roundn(1-P(1,2),-2))]}); figure; [R2 P2] = corrcoef([absMonEpmAnomaly(Indxx) absMonEpmAnomaly(Indyy)],[MonthlyMeanJunWindRotAng;MonthlyMeanJulWindRotAng]); scatter([absMonEpmAnomaly(Indxx) absMonEpmAnomaly(Indyy)],[MonthlyMeanJunWindRotAng;MonthlyMeanJulWindRotAng],80,'linewidth',2); set(gca,'box','on');set(gca,'fontsize',20); xlabel('Epm Anomaly (J/kg)');ylabel('Wind Rotation Angle (deg.)'); title({'June Epm Anomaly vs. Wind Rotation Angle'; ['coef = ' num2str(roundn(R2(1,2),-2)) ' ConfLev = ' num2str(roundn(1-P2(1,2),-2))]}); % polar vortex edge vs. std wind rot ang load('C:\Users\l1989\Desktop\PolarVortex\FigureofPolarVortex\PolarVortexEdge1980-2019.mat'); Alt = 45; [altInd,~] = find(PolarVortexEdge.ALTITUDE == Alt); for nYr = 2011:1:2019 Ind = find(PolarVortexEdge.DATE.Year == nYr & PolarVortexEdge.DATE.Month == 6); Ind2 = find(PolarVortexEdge.DATE.Year == nYr & PolarVortexEdge.DATE.Month == 7); MonthlyEdgeDisJun(nYr-2010) = nanmean(PolarVortexEdge.EDGE_DISTANCE(Ind,altInd)); MonthlyEdgeDisJul(nYr-2010) = nanmean(PolarVortexEdge.EDGE_DISTANCE(Ind2,altInd)); end [R1 P1] = corrcoef([MonthlyEdgeDisJun/1000 MonthlyEdgeDisJul/1000],[stdMonthlyMeanJunWindRotAng;stdMonthlyMeanJulWindRotAng]); figure; scatter([MonthlyEdgeDisJun/1000 MonthlyEdgeDisJul/1000],[stdMonthlyMeanJunWindRotAng;stdMonthlyMeanJulWindRotAng],80,'linewidth',2); set(gca,'box','on');set(gca,'fontsize',20); xlabel('Vortex Edge Distance from McMurdo (10^3 km)');ylabel('Wind Rotation Angle Std (deg.)'); title({'June Vortex Edge Distance vs. Wind Rotation Angle Std'; ['coef = ' num2str(roundn(R1(1,2),-2)) ' ConfLev = ' num2str(roundn(1-P1(1,2),-2))]}); %% plot Epm(Annual fitting) vs wind rotaion angle 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 SmthYdata = smooth(WinRotAng,14); p1 = plot(Time-dn_start,SmthYdata,'LineWidth',1.5);hold on;% Wind rotation angle is smoothed by 14 points(days) dtTime = datetime(Time,'convertfrom','datenum'); [~,IndJJ] = find(dtTime.Month == 6 | dtTime.Month == 7); % scatter(Time(IndJJ),SmthYdata(IndJJ),'filled'); ylim([0,360]);yticks([0 60 120 180 240 300 360]);ylabel('Wind Rotation Angle (deg.)','FontSize',30); yyaxis right % plot(timesortTotal,smooth(mean30to50EpmTotal,7),'LineWidth',2);% Epm is smoothed by 7 points for nYr = 2011:1:2019 hold on; dnAll = (datenum(nYr,1,1)-dn_start):1:(datenum(nYr,12,31)-dn_start); p2 = plot(dnAll,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))), '-r','LineWidth',1.5); end ylim([0,20]);xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); set(gca,'YMinorTick','on');set(gca,'TickDir','out'); 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',12); text(datenum('20110101','yyyymmdd')-dn_start,-3,'2011','fontsize',20);text(datenum('20120101','yyyymmdd')-dn_start,-3,'2012','fontsize',20); text(datenum('20130101','yyyymmdd')-dn_start,-3,'2013','fontsize',20);text(datenum('20140101','yyyymmdd')-dn_start,-3,'2014','fontsize',20); text(datenum('20150101','yyyymmdd')-dn_start,-3,'2015','fontsize',20);text(datenum('20160101','yyyymmdd')-dn_start,-3,'2016','fontsize',20); text(datenum('20170101','yyyymmdd')-dn_start,-3,'2017','fontsize',20);text(datenum('20180101','yyyymmdd')-dn_start,-3,'2018','fontsize',20); text(datenum('20190101','yyyymmdd')-dn_start,-3,'2019','fontsize',20); % legend('Wind Rotation Angle','E_{pm} (J/kg)','fontsize',20); % monthly mean wind rot angle IndJun = [];IndJul = [];TimeSortJun = [];TimeSortJul = []; MonthlyMeanJunWindRotAng = zeros(9,1);MonthlyMeanJulWindRotAng = zeros(9,1); stdMonthlyMeanJunWindRotAng = zeros(9,1);stdMonthlyMeanJulWindRotAng = zeros(9,1);stdJJMeanWindRotAng = zeros(9,1); stdMonthlyMeanWindRotAng = zeros(108,1); for nYr = 2011:1:2019 [~,IndJun] = find(dtTime.Year == nYr & dtTime.Month == 6); [~,IndJul] = find(dtTime.Year == nYr & dtTime.Month == 7); TimeSortJun(nYr-2010) = datenum(datetime(nYr,6,15)); TimeSortJul(nYr-2010) = datenum(datetime(nYr,7,15)); MonthlyMeanJunWindRotAng(nYr-2010) = mean(SmthYdata(IndJun)); MonthlyMeanJulWindRotAng(nYr-2010) = mean(SmthYdata(IndJul)); stdMonthlyMeanJunWindRotAng(nYr-2010) = std(SmthYdata(IndJun),0,'all'); stdMonthlyMeanJulWindRotAng(nYr-2010) = std(SmthYdata(IndJul),0,'all'); stdJJMeanWindRotAng(nYr-2010) = std([SmthYdata(IndJun);SmthYdata(IndJul)],0,'all'); end n = 1; for nYr = 2011:1:2019 for nMon = 1:1:12 [~,IndMon] = find(dtTime.Year == nYr & dtTime.Month == nMon); stdMonthlyMeanWindRotAng(n) = std(SmthYdata(IndMon),0,'all'); TimesortMon(n) = datenum(datetime(nYr,nMon,15)); n = n+1; end end yyaxis left p3 = scatter(TimeSortJun-dn_start,stdMonthlyMeanJunWindRotAng,'m','filled'); p4 = scatter(TimeSortJul-dn_start,stdMonthlyMeanJulWindRotAng.*8,'k','filled'); % p4 = scatter(TimeSortJun-dn_start,stdJJMeanWindRotAng,'g','filled'); p5 = plot(TimesortMon-dn_start,stdMonthlyMeanWindRotAng.*8,'-*b'); set(gca,'TickDir','out');set(gca,'TickLength',[0.006 0.0025]);set(gca,'fontsize',16); legend([p1 p2 p3 p4 p5],{'14 points smth Rot Ang'; 'Annual Fitting Epm'; 'Std of June'; 'Std of July'; 'Std of Monthly Wind Rot Ang'}); strlon1 = num2str(lon(lonRg(1))); strlon2 = num2str(lon(lonRg(end))); strlat1 = num2str(lat(latRg(1))); strlat2 = num2str(lat(latRg(end))); stralt1 = num2str(alt(altRg(1))); stralt2 = num2str(alt(altRg(end))); switch method case 1 title(['Wind Rot Angle: lon: ' strlon1 ' - ' strlon2 ' lat: ' strlat1 ' - ' strlat2 ' alt: ' stralt2 ' - ' stralt1 ' km']); case 2 title(['Wind Rot Angle(Mean for diff location): lon: ' strlon1 ' - ' strlon2 ' lat: ' strlat1 ' - ' strlat2 ' alt: ' stralt2 ' - ' stralt1 ' km']); case 3 title(['Wind Rot Angle(NewMethod): lon: ' strlon1 ' - ' strlon2 ' lat: ' strlat1 ' - ' strlat2 ' alt: ' stralt2 ' - ' stralt1 ' km'],'position',[(datenum('20150701','yyyymmdd')-dn_start),375],'FontSize',18); end %% plot 3km and 30km wind (set: altRg = [18 48]) scrsz = get(0, 'ScreenSize'); figure('position',[1,300,scrsz(3),scrsz(4)*0.8]); % 30km wind subplot('position',[0.05,0.11,0.9,0.35]); yyaxis left; Wind = sqrt(MeanU(1,:).^2+MeanV(1,:).^2); smthYData = smooth(Wind(end,:),14,'rlowess'); plot(Time-dn_start,smthYData,'linewidth',2);hold on; ylabel('Wind Speed ( m/s )');set(gca,'YMinorTick','on'); % [~,IndJJ] = find(dtTime.Month == 6 | dtTime.Month == 7); % scatter(Time(IndJJ)-dn_start,smthYData(IndJJ),'filled'); 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); ylim([0,20]);xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); set(gca,'YMinorTick','on');set(gca,'TickDir','out'); 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',12); text(datenum('20110101','yyyymmdd')-dn_start,-4,'2011','fontsize',20);text(datenum('20120101','yyyymmdd')-dn_start,-4,'2012','fontsize',20); text(datenum('20130101','yyyymmdd')-dn_start,-4,'2013','fontsize',20);text(datenum('20140101','yyyymmdd')-dn_start,-4,'2014','fontsize',20); text(datenum('20150101','yyyymmdd')-dn_start,-4,'2015','fontsize',20);text(datenum('20160101','yyyymmdd')-dn_start,-4,'2016','fontsize',20); text(datenum('20170101','yyyymmdd')-dn_start,-4,'2017','fontsize',20);text(datenum('20180101','yyyymmdd')-dn_start,-4,'2018','fontsize',20); text(datenum('20190101','yyyymmdd')-dn_start,-4,'2019','fontsize',20); strlon1 = num2str(lon(lonRg(1))); strlon2 = num2str(lon(lonRg(end))); strlat1 = num2str(lat(latRg(1))); strlat2 = num2str(lat(latRg(end))); stralt1 = num2str(alt(altRg(1))); stralt2 = num2str(alt(altRg(end))); set(gca,'TickDir','out');set(gca,'TickLength',[0.006 0.0025]);set(gca,'fontsize',16); title(['30km Wind(Mean for diff location): lon: ' strlon1 ' - ' strlon2 ' lat: ' strlat1 ' - ' strlat2],'position',[(datenum('20150701','yyyymmdd')-dn_start),105]); % 3km wind subplot('position',[0.05,0.59,0.9,0.35]); yyaxis left; Wind = sqrt(MeanU(end,:).^2+MeanV(end,:).^2); smthYData = smooth(Wind(end,:),14,'rlowess'); plot(Time-dn_start,smthYData,'linewidth',2);hold on; ylabel('Wind Speed ( m/s )');set(gca,'YMinorTick','on'); ylim([0 16]); % [~,IndJJ] = find(dtTime.Month == 6 | dtTime.Month == 7); % scatter(Time(IndJJ)-dn_start,smthYData(IndJJ),'filled'); 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); ylim([0,20]);xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); set(gca,'YMinorTick','on');set(gca,'TickDir','out'); 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',12); strlon1 = num2str(lon(lonRg(1))); strlon2 = num2str(lon(lonRg(end))); strlat1 = num2str(lat(latRg(1))); strlat2 = num2str(lat(latRg(end))); stralt1 = num2str(alt(altRg(1))); stralt2 = num2str(alt(altRg(end))); set(gca,'TickDir','out');set(gca,'TickLength',[0.006 0.0025]);set(gca,'fontsize',16); title(['3km Wind(Mean for diff location): lon: ' strlon1 ' - ' strlon2 ' lat: ' strlat1 ' - ' strlat2],'position',[(datenum('20150701','yyyymmdd')-dn_start),16.5]); %% monthly mean Epm and wind rotation angle corr n = 1; MonthlyEpm = [];MonthlyWRA = [];stdMonthlyWRA = []; WinRotAngSmth = smooth(WinRotAng,14,'rlowess'); for nYr = 2011:1:2019 for nMon = 1:1:12 IndMonEpm = find(timesort_datetime.Year == nYr & timesort_datetime.Month == nMon); IndMonWRA = find(dtSeries9yrs.Year == nYr & dtSeries9yrs.Month == nMon); MonthlyEpm(n) = nanmean(mean30to50EpmTotal(IndMonEpm)); MonthlyWRA(n) = nanmean(WinRotAngSmth(IndMonWRA)); stdMonthlyWRA(n) = std(WinRotAngSmth(IndMonWRA),0,'all'); n = n+1; end end IndNan = ~isnan(MonthlyEpm); [R1 P1] = corrcoef(MonthlyEpm(IndNan),MonthlyWRA(IndNan)); [R2 P2] = corrcoef(MonthlyEpm(IndNan),stdMonthlyWRA(IndNan)); figure; scatter(MonthlyEpm(IndNan),MonthlyWRA(IndNan)); n = 1; MonthlyJJEpm = [];MonthlyJJWRA = [];stdMonthlyJJWRA = []; for nYr = 2011:1:2019 for nMon = 7:1:7 IndMonJJEpm = find(timesort_datetime.Year == nYr & timesort_datetime.Month == nMon); IndMonJJWRA = find(dtSeries9yrs.Year == nYr & dtSeries9yrs.Month == nMon); MonthlyJJEpm(n) = nanmean(mean30to50EpmTotal(IndMonJJEpm)); MonthlyJJWRA(n) = nanmean(WinRotAngSmth(IndMonJJWRA)); stdMonthlyJJWRA(n) = std(WinRotAngSmth(IndMonJJWRA),0,'all'); n = n+1; end end IndNan = ~isnan(MonthlyJJEpm); [R3 P3] = corrcoef(MonthlyJJEpm(IndNan),MonthlyJJWRA(IndNan)); [R4 P4] = corrcoef(MonthlyJJEpm(IndNan),stdMonthlyJJWRA(IndNan)); figure; scatter(MonthlyJJEpm(IndNan),MonthlyJJWRA(IndNan)); %% monthly mean Epm and 3/30 km wind corr n = 1; Wind30 = sqrt(MeanU(1,:).^2+MeanV(1,:).^2); Wind3 = sqrt(MeanU(end,:).^2+MeanV(end,:).^2); Wind30 = smooth(Wind30,14,'rlowess'); Wind3 = smooth(Wind3,14,'rlowess'); MonthlyEpm = [];MonthlyWind30 = [];MonthlyWind3 = []; for nYr = 2011:1:2019 for nMon = 1:1:12 IndMonEpm = find(timesort_datetime.Year == nYr & timesort_datetime.Month == nMon); IndMonWRA = find(dtSeries9yrs.Year == nYr & dtSeries9yrs.Month == nMon); MonthlyEpm(n) = nanmean(mean30to50EpmTotal(IndMonEpm)); MonthlyWind30(n) = nanmean(Wind30(IndMonWRA)); MonthlyWind3(n) = nanmean(Wind3(IndMonWRA)); n = n+1; end end IndNan = ~isnan(MonthlyEpm); [R2 P2] = corrcoef(MonthlyEpm(IndNan),MonthlyWind30(IndNan)); [R3 P3] = corrcoef(MonthlyEpm(IndNan),MonthlyWind3(IndNan)); figure; scatter(MonthlyEpm(IndNan),MonthlyWind30(IndNan));hold on; scatter(MonthlyEpm(IndNan),MonthlyWind3(IndNan)); n = 1; Wind30 = sqrt(MeanU(1,:).^2+MeanV(1,:).^2); Wind3 = sqrt(MeanU(end,:).^2+MeanV(end,:).^2); MonthlyJJEpm = [];MonthlyJJWind30 = [];MonthlyJJWind3 = []; for nYr = 2011:1:2019 for nMon = 6:1:7 IndMonJJEpm = find(timesort_datetime.Year == nYr & timesort_datetime.Month == nMon); IndMonJJWRA = find(dtSeries9yrs.Year == nYr & dtSeries9yrs.Month == nMon); MonthlyJJEpm(n) = nanmean(mean30to50EpmTotal(IndMonJJEpm)); MonthlyJJWind30(n) = nanmean(Wind30(IndMonJJWRA)); MonthlyJJWind3(n) = nanmean(Wind3(IndMonJJWRA)); n = n+1; end end IndNan = ~isnan(MonthlyJJEpm); [R4 P4] = corrcoef(MonthlyJJEpm(IndNan),MonthlyJJWind30(IndNan)); [R5 P5] = corrcoef(MonthlyJJEpm(IndNan),MonthlyJJWind3(IndNan)); figure; scatter(MonthlyJJEpm(IndNan),MonthlyJJWind30(IndNan));hold on; scatter(MonthlyJJEpm(IndNan),MonthlyJJWind3(IndNan)); %% plot summer and winter respectively if ~exist('ft9yrs','var') load('C:\Users\l1989\Desktop\PolarVortex\Epm\paperFigures\FittingParam.mat'); end timesortWind = datetime(Time,'convertfrom','datenum'); clear WinterMonEpm; clear WinterMonWinRotAng; clear WinterEpm; clear dnTimesortWinterEpm; n = 1;nn = 1;nnn = 1; for nYr = 2011:1:2019 for nMon = 1:1:12 if nMon >=5 && nMon <= 8 IndEpm = find(timesort_datetime.Year == nYr & timesort_datetime.Month == nMon); if n == 1 dnTimesortWinterEpm = timesortTotal(IndEpm); WinterEpm = mean30to50EpmTotal(IndEpm); n = n+1; else dnTimesortWinterEpm = [dnTimesortWinterEpm;timesortTotal(IndEpm)]; WinterEpm = [WinterEpm mean30to50EpmTotal(IndEpm)]; end dnTimesortWinterMon(nnn) = datenum(nYr,nMon,15); IndWinterEpm = (timesort_datetime.Year == nYr & timesort_datetime.Month == nMon); WinterMonEpm(nnn) = mean(mean30to50EpmTotal(IndWinterEpm)); IndWinterWindRotAng = (timesortWind.Year == nYr & timesortWind.Month == nMon); WinterMonWinRotAng(nnn) = mean(WinRotAng(IndWinterWindRotAng)); dateVector = TimeSeries9yrs(IndWinterWindRotAng); WinterMonFitEpm(nnn) = mean(ft{nYr-2010}.Coefficients.Estimate(1) + ft{nYr-2010}.Coefficients.Estimate(2).*cos((2*pi/365).*(dateVector-ft{nYr-2010}.Coefficients.Estimate(3)))+ft{nYr-2010}.Coefficients.Estimate(4).*cos((4*pi/365).*(dateVector-ft{nYr-2010}.Coefficients.Estimate(5)))); nnn = nnn+1; end IndWinRot = (timesortWind.Year == nYr & timesortWind.Month == nMon); timesortWindMon(nn) = datenum(nYr,nMon,15); MonthlyMeanWindRotAng(nn) = mean(WinRotAng(IndWinRot)); IndEpmMon = (timesort_datetime.Year == nYr & timesort_datetime.Month == nMon); timesortEpmMon(nn) = datenum(nYr,nMon,15); MonthlyMeanEpm(nn) = mean(mean30to50EpmTotal(IndEpmMon)); nn = nn+1; end end f = figure; f.WindowState = 'maximized'; % scatter(timesortEpmMon, MonthlyMeanEpm,'filled');hold on; scatter(dnTimesortWinterMon,WinterMonEpm);hold on; scatter(dnTimesortWinterMon,WinterMonFitEpm,'filled'); ylim([0 20]);set(gca,'box','on');ylabel('E_p_m ( J/kg )'); yyaxis right plot(Time,smooth(WinRotAng(1,:),14,'rlowess'));% Wind rotation angle is smoothed by 14 points(days) scatter(dnTimesortWinterMon,WinterMonWinRotAng,'filled'); ylim([0 360]);yticks([0 60 120 180 240 300 360]); 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'),],'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',12); text(datenum('20110101','yyyymmdd'),-20,'2011','fontsize',20);text(datenum('20120101','yyyymmdd'),-20,'2012','fontsize',20); text(datenum('20130101','yyyymmdd'),-20,'2013','fontsize',20);text(datenum('20140101','yyyymmdd'),-20,'2014','fontsize',20); text(datenum('20150101','yyyymmdd'),-20,'2015','fontsize',20);text(datenum('20160101','yyyymmdd'),-20,'2016','fontsize',20); text(datenum('20170101','yyyymmdd'),-20,'2017','fontsize',20);text(datenum('20180101','yyyymmdd'),-20,'2018','fontsize',20); text(datenum('20190101','yyyymmdd'),-20,'2019','fontsize',20); ylabel('Wind roation angle ( deg. )'); legend('monthly mean Epm','fitting Epm','Wind Rotation Angle','fontsize',16); xlim([datenum('20110101','yyyymmdd'),datenum('20191231','yyyymmdd')]); % correlation bewteen Epm and wind rotaion angle figure; scatter(WinterMonWinRotAng,WinterMonEpm,'r');hold on; scatter(WinterMonWinRotAng,WinterMonFitEpm,'b'); legend({'monthly mean Epm' 'monthly fitting Epm'}); title([{['\color{red}coef: ' num2str(R(1,2)) 'conf level: ' num2str(P(1,2))]} {['\color{blue}coef:' num2str(Rfit(1,2)) 'conf level: ' num2str(Pfit(1,2))]}]); indNan = ~isnan(WinterMonEpm); [R P] = corrcoef(WinterMonEpm(indNan), WinterMonWinRotAng(indNan)); [Rfit Pfit] = corrcoef(WinterMonFitEpm, WinterMonWinRotAng); xlabel('Monthly mean wind rotation angle ( deg. )'); ylabel('Epm ( J/kg )');set(gca,'box','on'); %% wind contour ContourU = zeros(length(alt),length(U)); ContourSpeed = zeros(length(alt),length(U)); for nCnt = 1:1:length(U) ContourU(:,nCnt) = U{nCnt}(lonRg,latRg,:); ContourSpeed(:,nCnt) = sqrt(U{nCnt}(lonRg,latRg,:).^2+V{nCnt}(lonRg,latRg,:).^2); end %% dn_start = datenum(2010,12,31); f = figure; f.WindowState = 'maximized'; [h1, h2]=contourf(TimeSeries9yrs-dn_start,alt,ContourSpeed,'linestyle','none'); set(h2,'LevelStep',2); handlecolor=colorbar; colormap jet; set(gca,'fontsize',16) set(gca,'CLim',[0, 120]); % set(get(handlecolor,'ylabel'),'String', 'K'); set(get(handlecolor,'ylabel'),'fontsize', 16); 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',12); set(gca,'TickDir','out');ylabel('altitude ( km )');set(gca,'TickLength',[0.005 0.025]); text(datenum('20110101','yyyymmdd')-dn_start,-1,'2011','FontSize',18);text(datenum('20120101','yyyymmdd')-dn_start,-1,'2012','FontSize',18); text(datenum('20130101','yyyymmdd')-dn_start,-1,'2013','FontSize',18);text(datenum('20140101','yyyymmdd')-dn_start,-1,'2014','FontSize',18); text(datenum('20150101','yyyymmdd')-dn_start,-1,'2015','FontSize',18);text(datenum('20160101','yyyymmdd')-dn_start,-1,'2016','FontSize',18); text(datenum('20170101','yyyymmdd')-dn_start,-1,'2017','FontSize',18);text(datenum('20180101','yyyymmdd')-dn_start,-1,'2018','FontSize',18); text(datenum('20190101','yyyymmdd')-dn_start,-1,'2019','FontSize',18); title('McMurdo Wind Speed'); %% calculate and plot 3km wind and 30km wind - only use lat 25(-78) 3km - lev 56 30km - lev 26 WindSpeed = sqrt(U{1}.^2 + V{1}.^2); load('C:\Users\l1989\Desktop\PolarVortex\Epm\Epm2011-2019.mat'); f = figure(2); f.WindowState = 'maximized'; %plot wind @ 3km subplot(2,1,1); yyaxis left plot(Time,smooth(WindSpeed(56,:),14));% smoothing by 14 points(days) ylim([3,20]);ylabel('Wind @3km','FontSize',30); yyaxis right plot(timesortTotal,smooth(mean30to50EpmTotal,7));% smoothing by 7 points ylim([0,20]);xlim([datenum('20110101','yyyymmdd'),datenum('20191231','yyyymmdd')]); set(gca,'YMinorTick','on');set(gca,'TickDir','in'); 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'),],'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',12); text(datenum('20110101','yyyymmdd'),-2.5,'2011','fontsize',20);text(datenum('20120101','yyyymmdd'),-2.5,'2012','fontsize',20); text(datenum('20130101','yyyymmdd'),-2.5,'2013','fontsize',20);text(datenum('20140101','yyyymmdd'),-2.5,'2014','fontsize',20); text(datenum('20150101','yyyymmdd'),-2.5,'2015','fontsize',20);text(datenum('20160101','yyyymmdd'),-2.5,'2016','fontsize',20); text(datenum('20170101','yyyymmdd'),-2.5,'2017','fontsize',20);text(datenum('20180101','yyyymmdd'),-2.5,'2018','fontsize',20); text(datenum('20190101','yyyymmdd'),-2.5,'2019','fontsize',20); legend('Wind @3km','E_{pm} (J/kg)','location','northeast','fontsize',15); %plot wind at 30km subplot(2,1,2); yyaxis left plot(Time,smooth(WindSpeed(26,:),14));% smoothing by 14 points(days) ylim([0,120]);ylabel('Wind @30km','FontSize',30); yyaxis right plot(timesortTotal,smooth(mean30to50EpmTotal,7));% smoothing by 7 points ylim([0,20]);xlim([datenum('20110101','yyyymmdd'),datenum('20191231','yyyymmdd')]); set(gca,'YMinorTick','on');set(gca,'TickDir','in'); 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'),],'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',12); text(datenum('20110101','yyyymmdd'),-2.5,'2011','fontsize',20);text(datenum('20120101','yyyymmdd'),-2.5,'2012','fontsize',20); text(datenum('20130101','yyyymmdd'),-2.5,'2013','fontsize',20);text(datenum('20140101','yyyymmdd'),-2.5,'2014','fontsize',20); text(datenum('20150101','yyyymmdd'),-2.5,'2015','fontsize',20);text(datenum('20160101','yyyymmdd'),-2.5,'2016','fontsize',20); text(datenum('20170101','yyyymmdd'),-2.5,'2017','fontsize',20);text(datenum('20180101','yyyymmdd'),-2.5,'2018','fontsize',20); text(datenum('20190101','yyyymmdd'),-2.5,'2019','fontsize',20); legend('Wind @30km','E_{pm} (J/kg)','location','northeast','fontsize',15); %% calculate and plot correlation coefficient %{ In order to use func 'corrcoef/corr', the two arrays need to have the same size. So we are using monthly mean wind rotation angle and Epm. %} clc;clear;close all; load('test1.mat'); load('MERRA2_Diural_UVT_ModelLevel_20110101-20190331.mat');% wind rotation angle info load('C:\Users\l1989\Desktop\PolarVortex\Epm\Epm2011-2019.mat');% Epm info % take mean of all the Epm each month and Monthly mean wind rotation angle nCount = 1; for i = 2011:1:2015 indY1 = find(year(timesort_datetime) == i); indY2 = find(year(timesort) == i); for j = 1:1:12 indM1 = find(month(timesort_datetime) == j); ind1 = intersect(indY1,indM1); indM2 = find(month(timesort) == j); ind2 = intersect(indY2,indM2); if isempty(ind1) MonthlyEpm(nCount) = NaN; else MonthlyEpm(nCount) = mean(mean30to50EpmTotal(ind1)); end MonthlyWindRotAng(nCount) = mean(WinRotAng(1,ind2)); MonthlyWindSpeed3km(nCount) = mean(WindSpeed(56,ind2)); MonthlyWindSpeed30km(nCount) = mean(WindSpeed(26,ind2)); nCount = nCount+1; end end % the monthly wind rotation angle array is smaller than the Epm MonthlyEpm(isnan(MonthlyWindRotAng))=[]; MonthlyWindRotAng(isnan(MonthlyWindRotAng))=[]; MonthlyWindSpeed3km(isnan(MonthlyWindSpeed3km))=[]; MonthlyWindSpeed30km(isnan(MonthlyWindSpeed30km))=[]; % some month don't have Epm values % MonthlyWindRotAng(isnan(MonthlyEpm))=[]; % MonthlyEpm(isnan(MonthlyEpm))=[]; %% % calcucate correlation coefficient [R1 P1] = corrcoef(MonthlyEpm,MonthlyWindRotAng,'Rows','pairwise'); [R2 P2] = corrcoef(MonthlyEpm,MonthlyWindSpeed3km,'Rows','pairwise'); [R3 P3] = corrcoef(MonthlyEpm,MonthlyWindSpeed30km,'Rows','pairwise'); % do plotting len = size(MonthlyEpm,2); f = figure(3); f.WindowState = 'maximized'; yyaxis left plot(1:len,MonthlyWindRotAng,'-o'); ylim([0,360]);ylabel('Wind Rotation Angle','FontSize',30); yyaxis right plot(1:len,MonthlyEpm,'-o'); ylim([0,20]);ylabel('E_{pm} (J/kg)','FontSize',30); set(gca,'xtick',1:1:len,'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',12); text(1,-1,'2011','fontsize',20);text(13,-1,'2012','fontsize',20); text(25,-1,'2013','fontsize',20);text(37,-1,'2014','fontsize',20); text(49,-1,'2015','fontsize',20);text(61,-1,'2016','fontsize',20); text(73,-1,'2017','fontsize',20);text(85,-1,'2018','fontsize',20); text(97,-1,'2019','fontsize',20); title(['2011-2015 Correlation Coefficient: ' num2str(R1(1,2)) ' P: ' num2str(P1(1,2))]); f = figure(4); f.WindowState = 'maximized'; yyaxis left plot(1:len,MonthlyWindSpeed3km,'-o'); ylim([3,20]);ylabel('Wind @3km','FontSize',30); yyaxis right plot(1:len,MonthlyEpm,'-o'); ylim([0,20]);ylabel('E_{pm} (J/kg)','FontSize',30); set(gca,'xtick',1:1:len,'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',12); text(1,-1,'2011','fontsize',20);text(13,-1,'2012','fontsize',20); text(25,-1,'2013','fontsize',20);text(37,-1,'2014','fontsize',20); text(49,-1,'2015','fontsize',20);text(61,-1,'2016','fontsize',20); text(73,-1,'2017','fontsize',20);text(85,-1,'2018','fontsize',20); text(97,-1,'2019','fontsize',20); title(['2011-2015 Correlation Coefficient: ' num2str(R2(1,2)) ' P: ' num2str(P2(1,2))]); f = figure(5); f.WindowState = 'maximized'; yyaxis left plot(1:len,MonthlyWindSpeed30km,'-o'); ylim([0,120]);ylabel('Wind @30km','FontSize',30); yyaxis right plot(1:len,MonthlyEpm,'-o'); ylim([0,20]);ylabel('E_{pm} (J/kg)','FontSize',30); set(gca,'xtick',1:1:len,'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',12); text(1,-1,'2011','fontsize',20);text(13,-1,'2012','fontsize',20); text(25,-1,'2013','fontsize',20);text(37,-1,'2014','fontsize',20); text(49,-1,'2015','fontsize',20);text(61,-1,'2016','fontsize',20); text(73,-1,'2017','fontsize',20);text(85,-1,'2018','fontsize',20); text(97,-1,'2019','fontsize',20); title(['2011-2015 Correlation Coefficient: ' num2str(R3(1,2)) ' P: ' num2str(P3(1,2))]); %% testing wind rotation angle from McMurdo to equator clc;close all;clear;fclose all; mainpath = 'E:\New folder\'; cd(mainpath); path = 'MERRA2_20110101-20190630_Diurnal_UV_Temp_ModelLevel_-90to0Latitude\'; % ncdisp([path 'MERRA2_400.inst6_3d_ana_Nv.20110101.SUB.nc']); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Source: % E:\New folder\MERRA2_20110101-20190630_Diurnal_UV_Temp_ModelLevel_-90to0Latitude\MERRA2_400.inst6_3d_ana_Nv.20110101.SUB.nc % Format: % netcdf4 % Global Attributes: % CDI = 'Climate Data Interface version 1.6.9 (http://mpimet.mpg.de/cdi)' % history = 'Sun Jun 16 13:34:20 2019: cdo -s -L -f nc4 -timmean -sellonlatbox,-180.0,180.0,-80.0,0.0 -selname,T,U,V /tmpdata/regridder/s_1560692043_13241/MERRA2_400.inst6_3d_ana_Nv.20110101.nc4 /tmpdata/regridder/s_1560692043_13241/MERRA2_400.inst6_3d_ana_Nv.20110101.SUB.nc4' % Conventions = 'CF-1' % ks = 40 % nstep = 15760 % ptop = 1 % pint = 15039.3027 % ak = [1 2.00000024 3.27000046 4.75850105 6.60000134 8.93450165 11.9703016 15.9495029 21.134903 27.8526058 36.5041084 47.5806084 61.6779099 79.5134125 101.944023 130.051025 165.079025 208.49704 262.021057 327.643066 407.657104 504.680115 621.680115 761.984192 929.294189 1127.69019 1364.34021 1645.71033 1979.1604 2373.04053 2836.78052 3381.00073 4017.54102 4764.39111 5638.79102 6660.34131 7851.23145 9236.57227 10866.3018 12783.7031 15039.3027 17693.0039 20119.2012 21686.502 22436.3008 22389.8008 21877.5977 21214.998 20325.8984 19309.6953 18161.8965 16960.8965 15625.9961 14290.9951 12869.5938 11895.8623 10918.1709 9936.52148 8909.99219 7883.42188 7062.19824 6436.26367 5805.32129 5169.61084 4533.90088 3898.20093 3257.08081 2609.20068 1961.31055 1313.48035 659.375244 4.80482578 0] % bk = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 8.1754e-09 0.00696 0.02801 0.06372 0.1136 0.15622 0.20035 0.24674 0.2944 0.34338 0.39289 0.44374 0.49459 0.5463 0.58104 0.61582 0.65063 0.6859 0.72117 0.74938 0.77064 0.79195 0.8133 0.83466 0.85602 0.87743 0.89891 0.92039 0.94187 0.96341 0.98495 1] % History = 'Original file generated: Mon Feb 16 04:38:36 2015 GMT' % Comment = 'GMAO filename: d5124_m2_jan10.inst6_3d_ana_Nv.20110101_.nc4' % Filename = 'MERRA2_400.inst6_3d_ana_Nv.20110101.nc4' % Institution = 'NASA Global Modeling and Assimilation Office' % References = 'http://gmao.gsfc.nasa.gov' % Format = 'NetCDF-4/HDF-5' % SpatialCoverage = 'global' % VersionID = '5.12.4' % TemporalRange = '1980-01-01 -> 2016-12-31' % identifier_product_doi_authority = 'http://dx.doi.org/' % ShortName = 'M2I6NVANA' % GranuleID = 'MERRA2_400.inst6_3d_ana_Nv.20110101.nc4' % ProductionDateTime = 'Original file generated: Mon Feb 16 04:38:36 2015 GMT' % LongName = 'MERRA2 inst6_3d_ana_Nv: 3d,6-Hourly,Instantaneous,Model-Level,Analysis,Analyzed Meteorological Fields' % Title = 'MERRA2 inst6_3d_ana_Nv: 3d,6-Hourly,Instantaneous,Model-Level,Analysis,Analyzed Meteorological Fields' % SouthernmostLatitude = '-90.0' % NorthernmostLatitude = '90.0' % WesternmostLongitude = '-180.0' % EasternmostLongitude = '179.375' % LatitudeResolution = '0.5' % LongitudeResolution = '0.625' % DataResolution = '0.5 x 0.625 (72 native layers)' % Contact = 'http://gmao.gsfc.nasa.gov' % identifier_product_doi = '10.5067/IUUF4WB9FT4W' % RangeBeginningDate = '2011-01-01' % RangeBeginningTime = '00:00:00.000000' % RangeEndingDate = '2011-01-01' % RangeEndingTime = '18:00:00.000000' % CDO = 'Climate Data Operators version 1.6.9 (http://mpimet.mpg.de/cdo)' % Dimensions: % time = 1 (UNLIMITED) % lev = 72 % lat = 161 % lon = 576 % bnds = 2 % Variables: % T % Size: 576x161x72x1 % Dimensions: lon,lat,lev,time % Datatype: single % Attributes: % standard_name = 'air_temperature' % long_name = 'Air temperature' % units = 'K' % _FillValue = 999999986991104 % missing_value = 999999986991104 % fmissing_value = 999999986991104 % vmax = 999999986991104 % vmin = -999999986991104 % U % Size: 576x161x72x1 % 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 % V % Size: 576x161x72x1 % Dimensions: lon,lat,lev,time % Datatype: single % Attributes: % standard_name = 'northward_wind' % long_name = 'Northward wind component' % units = 'm/s' % _FillValue = 999999986991104 % missing_value = 999999986991104 % fmissing_value = 999999986991104 % vmax = 999999986991104 % vmin = -999999986991104 % lat % Size: 161x1 % Dimensions: lat % Datatype: double % Attributes: % standard_name = 'latitude' % long_name = 'latitude' % units = 'degrees_north' % axis = 'Y' % lev % Size: 72x1 % Dimensions: lev % Datatype: double % Attributes: % long_name = 'vertical level' % units = 'layer' % 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 2011-01-01 00:00:00' % calendar = 'standard' % time_bnds % Size: 2x1 % Dimensions: bnds,time % Datatype: double %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % load altitude(lev),lon,lat a=load([path 'MERRA2_ModelLevel 72.txt']); alt = struct('ModelLev',a(:,1),'Pressure',a(:,2)); alt.Altitude = atmospalt(alt.Pressure*100); lon = ncread([path 'MERRA2_400.inst6_3d_ana_Nv.20110101.SUB.nc'],'lon'); lat = ncread([path 'MERRA2_400.inst6_3d_ana_Nv.20110101.SUB.nc'],'lat'); % lat -80:0.5:0 lon -180:0.625:179.375 % read nc files - lat: 1 to 161 is -80 to 0 lon: 556 is 166.875 lev: 26 to 43 % are 11 to 30km lonNum = 558;% longtitude number DateStart = datetime([2011 01 01]); DateEnd = datetime([2019 03 31]); mon=months(datestr(DateStart),datestr(DateEnd))+1; day=days(DateEnd-DateStart); filename0 = [path 'MERRA2_400.inst6_3d_ana_Nv.'];filename1 = '.SUB.nc'; U = cell(161,1); V = cell(161,1); T = cell(161,1); Time = cell(1,1); for i=1:1:mon Days = eomday(year(DateStart),month(DateStart)); for j = 0:1:(Days-1) DateTemp = DateStart+days(j); str = datestr(DateTemp,'yyyymmdd'); Utemp = ncread([filename0 str filename1],'U'); Vtemp = ncread([filename0 str filename1],'V'); Ttemp = ncread([filename0 str filename1],'T'); Timetemp = datenum(str,'yyyymmdd'); for k = 1:1:161 if i == 1 && j == 0 U{k} = reshape(Utemp(lonNum,k,:,1),72,1); V{k} = reshape(Vtemp(lonNum,k,:,1),72,1); T{k} = reshape(Ttemp(lonNum,k,:,1),72,1); else U{k} = [U{k} reshape(Utemp(lonNum,k,:,1),72,1)]; V{k} = [V{k} reshape(Vtemp(lonNum,k,:,1),72,1)]; T{k} = [T{k} reshape(Ttemp(lonNum,k,:,1),72,1)]; end end if i == 1 && j == 0 Time = Timetemp; else Time = [Time;Timetemp]; end end DateStart = DateStart + days(Days); disp(i); end timesort = datetime(Time,'convertfrom','datenum'); if ~exist(['MERRA2_Diural_UVT_ModelLevel_lat(-80,0)_lon(' num2str(lon(lonNum)) ')_20110101-20190331.mat'],'file') save(['MERRA2_Diural_UVT_ModelLevel_lat(-80,0)_lon(' num2str(lon(lonNum)) ')_20110101-20190331.mat'],'lev','T','U','V','Time','timesort','lon','lat'); end %% load file clc;clear;close all; if exist('MERRA2_Diural_UVT_ModelLevel_lat(-80,0)_lon(166.875)_20110101-20190331.mat','file') load('MERRA2_Diural_UVT_ModelLevel_lat(-80,0)_lon(166.875)_20110101-20190331.mat'); else disp('Can not find database!'); return; end %% U including 161*0.5 in latitude(-80 to 0), 72 level in altitude and 3012 days(20110101 to 20190331) ang = cell(size(U,1),1); WinRotAng = cell(size(U,1),1); for i = 1:1:size(U,1) % deal with the outlier nanindex=(U{i}>9999 | U{i}<-9999); U{i}(nanindex) = NaN; nanindex=(V{i}>9999 | V{i}<-9999); V{i}(nanindex) = NaN; nanindex=(T{i}>9999 | T{i}<-9999); T{i}(nanindex) = NaN; %calc wind angle ang{i} = atan2(V{i},U{i}); %calc wind rotation angle, lev: 26 to 43 are 11 to 30km WinRotAng{i} = (max(ang{i}(26:43,:))-min(ang{i}(26:43,:)))/pi*180; %rotation angle > 360 equal to 360 index=WinRotAng{i}>360; WinRotAng{i}(index)=360; end %% plot Epm vs wind rotaion angle load('C:\Users\l1989\Desktop\PolarVortex\Epm\Epm2011-2019.mat'); f = figure(1); f.WindowState = 'maximized'; latNum = 5;% latitude number yyaxis left plot(Time,smooth(WinRotAng{latNum}(:),14));% Wind rotation angle is smoothed by 14 points(days) ylim([0,360]);ylabel(['Wind Rotation Angle @' num2str(lat(latNum)) 'degree'],'FontSize',30); yyaxis right plot(timesortTotal,smooth(mean30to50EpmTotal,7));% Epm is smoothed by 7 points ylim([0,20]);xlim([datenum('20110101','yyyymmdd'),datenum('20191231','yyyymmdd')]); set(gca,'YMinorTick','on');set(gca,'TickDir','in'); 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'),],'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',12); text(datenum('20110101','yyyymmdd'),-1,'2011','fontsize',20);text(datenum('20120101','yyyymmdd'),-1,'2012','fontsize',20); text(datenum('20130101','yyyymmdd'),-1,'2013','fontsize',20);text(datenum('20140101','yyyymmdd'),-1,'2014','fontsize',20); text(datenum('20150101','yyyymmdd'),-1,'2015','fontsize',20);text(datenum('20160101','yyyymmdd'),-1,'2016','fontsize',20); text(datenum('20170101','yyyymmdd'),-1,'2017','fontsize',20);text(datenum('20180101','yyyymmdd'),-1,'2018','fontsize',20); text(datenum('20190101','yyyymmdd'),-1,'2019','fontsize',20); legend('Wind Rotation Angle','E_{pm} (J/kg)','location','north','fontsize',20); %% try to see QBO in zonal wind % get zonal wind monthly mean ZonalWindMonthlyMean = cell(size(U,1),1); timeMonthlySort = NaT(1,108); for l = 1:1:size(U,1) nCount = 1; for i = 2011:1:2019 indY = find(year(timesort) == i); for j = 1:1:12 indM = find(month(timesort) == j); ind = intersect(indY,indM); if isnan(ind) return; else ZonalWindMonthlyMean{l}(1:72,nCount) = mean(U{l}(1:72,ind),2); timeMonthlySort(nCount) = datetime(i,j,15); nCount = nCount+1; end end end end %% clc;close all;clear; mainpath = 'E:\New folder\'; cd(mainpath); path = 'MERRA2_20110101-20190630_Diurnal_UV_Temp_ModelLevel_aroundMcMurdo\'; % load altitude(lev),lon,lat a=load('MERRA2_20110101-20190630_Diurnal_UV_Temp_ModelLevel_aroundMcMurdo\MERRA2_ModelLevel 72.txt'); alt = struct('ModelLev',a(:,1),'Pressure',a(:,2)); alt.Altitude = atmospalt(alt.Pressure*100); lon = ncread('MERRA2_20110101-20190630_Diurnal_UV_Temp_ModelLevel_aroundMcMurdo\MERRA2_400.inst6_3d_ana_Nv.20110101.SUB.nc','lon'); lat = ncread('MERRA2_20110101-20190630_Diurnal_UV_Temp_ModelLevel_aroundMcMurdo\MERRA2_400.inst6_3d_ana_Nv.20110101.SUB.nc','lat'); % read nc files - lat: 25 is -78 26 is -77.5 lon: 12 is 166.875 lev: 26 to 43 % are 11 to 30km DateStart = datetime([2011 01 01]); DateEnd = datetime([2019 03 31]); mon=months(datestr(DateStart),datestr(DateEnd))+1; day=days(DateEnd-DateStart); filename0 = 'MERRA2_20110101-20190630_Diurnal_UV_Temp_ModelLevel_aroundMcMurdo\MERRA2_400.inst6_3d_ana_Nv.';filename1 = '.SUB.nc'; U = cell(2,2); V = cell(2,2); T = cell(2,2); Time = cell(1,1); latNum1 = 39; latNum2 = 40; lonNum1 = 12; lonNum2 = 13; for i=1:1:mon Days = eomday(year(DateStart),month(DateStart)); for j = 0:1:(Days-1) DateTemp = DateStart+days(j); str = datestr(DateTemp,'yyyymmdd'); Utemp = ncread([filename0 str filename1],'U'); Vtemp = ncread([filename0 str filename1],'V'); Ttemp = ncread([filename0 str filename1],'T'); Timetemp = datenum(str,'yyyymmdd'); if i == 1 && j == 0 U{1,1} = reshape(Utemp(lonNum1,latNum1,:,1),72,1); U{1,2} = reshape(Utemp(lonNum1,latNum2,:,1),72,1); V{1,1} = reshape(Vtemp(lonNum1,latNum1,:,1),72,1); V{1,2} = reshape(Vtemp(lonNum1,latNum2,:,1),72,1); T{1,1} = reshape(Ttemp(lonNum1,latNum1,:,1),72,1); T{1,2} = reshape(Ttemp(lonNum1,latNum2,:,1),72,1); U{2,1} = reshape(Utemp(lonNum2,latNum1,:,1),72,1); U{2,2} = reshape(Utemp(lonNum2,latNum2,:,1),72,1); V{2,1} = reshape(Vtemp(lonNum2,latNum1,:,1),72,1); V{2,2} = reshape(Vtemp(lonNum2,latNum2,:,1),72,1); T{2,1} = reshape(Ttemp(lonNum2,latNum1,:,1),72,1); T{2,2} = reshape(Ttemp(lonNum2,latNum2,:,1),72,1); Time = Timetemp; else U{1,1} = [U{1,1} reshape(Utemp(lonNum1,latNum1,:,1),72,1)]; U{1,2} = [U{1,2} reshape(Utemp(lonNum1,latNum2,:,1),72,1)]; V{1,1} = [V{1,1} reshape(Vtemp(lonNum1,latNum1,:,1),72,1)]; V{1,2} = [V{1,2} reshape(Vtemp(lonNum1,latNum2,:,1),72,1)]; T{1,1} = [T{1,1} reshape(Ttemp(lonNum1,latNum1,:,1),72,1)]; T{1,2} = [T{1,2} reshape(Ttemp(lonNum1,latNum2,:,1),72,1)]; U{2,1} = [U{2,1} reshape(Utemp(lonNum2,latNum1,:,1),72,1)]; U{2,2} = [U{2,2} reshape(Utemp(lonNum2,latNum2,:,1),72,1)]; V{2,1} = [V{2,1} reshape(Vtemp(lonNum2,latNum1,:,1),72,1)]; V{2,2} = [V{2,2} reshape(Vtemp(lonNum2,latNum2,:,1),72,1)]; T{2,1} = [T{2,1} reshape(Ttemp(lonNum2,latNum1,:,1),72,1)]; T{2,2} = [T{2,2} reshape(Ttemp(lonNum2,latNum2,:,1),72,1)]; Time = [Time;Timetemp]; end end DateStart = DateStart + days(Days); disp(i); end timesort = datetime(Time,'convertfrom','datenum'); if ~exist(['MERRA2_Diural_UVT_ModelLevel72_lat(' num2str(lat(latNum1)) ',' num2str(lat(latNum2)) ')_lon(' num2str(lon(lonNum1)) ',' num2str(lon(lonNum2)) ')_20110101-20190331.mat'],'file') save(['MERRA2_Diural_UVT_ModelLevel72_lat(' num2str(lat(latNum1)) ',' num2str(lat(latNum2)) ')_lon(' num2str(lon(lonNum1)) ',' num2str(lon(lonNum2)) ')_20110101-20190331.mat'],'lev','T','U','V','Time','timesort'); end %% calculate Wind rotation angle = max angle(11 to 30km) - min angle(11 to 30km) %{ Wind rotaion angle is defined as the biggest wind angle difference between 11 to 30km. So first calculate the wind angle for all altitude, then using the maximum angle minus the minimum angle. Smoothing is hiding in the plotting region. %} clc;close all;clear; load('MERRA2_Diural_UVT_ModelLevel72_lat(-78,-77.5)_lon(166.875,167.5)_20110101-20190331.mat'); ang = cell(2,2); WinRotAng = cell(2,2); for i = 1:1:2 for j = 1:1:2 % deal with the outlier nanindex=(U{i,j}>9999 | U{i,j}<-9999); U{i,j}(nanindex) = NaN; nanindex=(V{i,j}>9999 | V{i,j}<-9999); V{i,j}(nanindex) = NaN; nanindex=(T{i,j}>9999 | T{i,j}<-9999); T{i,j}(nanindex) = NaN; % calc angle ang{i,j} = atan2(V{i,j}, U{i,j}); %calc wind rotation angle - lev: 26 to 43 are 11 to 30km WinRotAng{i,j} = (max(ang{i,j}(26:43,:))-min(ang{i,j}(26:43,:)))/pi*180; index=WinRotAng{i,j}>360; WinRotAng{i,j}(index)=360; end end %% plot Epm vs wind rotation angle(different latitude and longtitude) load('C:\Users\l1989\Desktop\PolarVortex\Epm\Epm2011-2019.mat'); f = figure(1); f.WindowState = 'maximized'; yyaxis left plot(Time,smooth(WinRotAng{2,1}(:),14));% Wind rotation angle is smoothed by 14 points(days) ylim([0,360]);ylabel('Wind Rotation Angle@(-77.5,166.875)','FontSize',30); yyaxis right plot(timesortTotal,smooth(mean30to50EpmTotal,7));% Epm is smoothed by 7 points ylim([0,20]);xlim([datenum('20110101','yyyymmdd'),datenum('20191231','yyyymmdd')]); set(gca,'YMinorTick','on');set(gca,'TickDir','in'); 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'),],'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',12); text(datenum('20110101','yyyymmdd'),-1,'2011','fontsize',20);text(datenum('20120101','yyyymmdd'),-1,'2012','fontsize',20); text(datenum('20130101','yyyymmdd'),-1,'2013','fontsize',20);text(datenum('20140101','yyyymmdd'),-1,'2014','fontsize',20); text(datenum('20150101','yyyymmdd'),-1,'2015','fontsize',20);text(datenum('20160101','yyyymmdd'),-1,'2016','fontsize',20); text(datenum('20170101','yyyymmdd'),-1,'2017','fontsize',20);text(datenum('20180101','yyyymmdd'),-1,'2018','fontsize',20); text(datenum('20190101','yyyymmdd'),-1,'2019','fontsize',20); legend('Wind Rotation Angle','E_{pm} (J/kg)','location','north','fontsize',20);