A hypothesis of global warming influence on the rapid intensification of western North Pacific tropical cyclones
Nam-Young Kang and James Elsner
Abstract
In association with global warming, rapid intensification (RI) of tropical cyclones (TCs) is one of the rising concerns in our world. With reference to the previous studies on the TC intensification with fewer occurrences, so called the efficiency of intensity (EINT), this study poses a hypothesis that TC RI in the western North Pacific is related to the EINT environment favored by global warming. EINT environment represents the strongly anomalous high over most unstable tropical atmosphere, which supports the idea of blasting effect. Investigation during June to November over the 30-year period (2016–2015) reveals that EINT explains 93.0~% of the global warming influence on the the increasing proportion of RI-experiencing TCs. Consequently, global warming leads to larger proportion of RI events, fewers TCs, and thus neutral trend in the number of RI events after all.
Table 1. Conversion table of JMA LMIs equivalent to JTWC LMIs
# Sources
#--------------------------------------------------------------------------------
time=1986:2015
AB=read.table('Data/ABtrack2016.dat',T)
AB1=AB[(AB$bst=='JT' & AB$yr>=time[1] & AB$yr<=rev(time)[1] & AB$mo>=6 & AB$mo<=11),]
AB2=AB[(AB$bst=='JM' & AB$yr>=time[1] & AB$yr<=rev(time)[1] & AB$mo>=6 & AB$mo<=11),]
H1=hist(AB1$lmw,breaks=seq(32.5,182.5,5),prob=T,xlab='JTWC LMI')
H2=hist(AB2$lmw,breaks=seq(32.5,182.5,5),prob=T,xlab='JMA LMI')
Djt=cumsum(H1$counts)/sum(H1$count) #DensiTy for JTWC
Djm=cumsum(H2$counts)/sum(H2$count) #Density for JMA
# input
#--------------------------------------------------------------------------------
INTlev=seq(35,180,5)
n = length(seq(35,180,1))
# interpolation
#--------------------------------------------------------------------------------
spline1=spline(INTlev,Djm, n)
spline2=spline(INTlev,Djt, n)
X=spline1$x
Yjm=spline1$y
Yjt=spline2$y
Dimsi=NULL; for(i in INTlev) Dimsi=c(Dimsi,Yjm[X==i]) #get JMA densify at 5-kt intervals
Limsi=NULL; for(i in Dimsi) Limsi=c(Limsi,X[which.min(abs(Yjt-i))]) #find eq JT LMI at the nearest density
plot(-99,-99,xlim=c(35,180),ylim=c(0,1), xlab='LMI (kt)', ylab='Prob level')
points(H1$mids,Djt)
lines(H1$mids,Djt,col=4)
points(H2$mids,Djm,pch=16)
lines(H2$mids,Djm,col=1)
df=data.frame(JMALMI=INTlev,ProbLev=round(100*Dimsi,1), EqJMALMI=Limsi)[1:19,]
points(df$EqJMALMI,df$ProbLev,pch=16,cex=1,col=4) 
df## JMALMI ProbLev EqJMALMI
## 1 35 6.6 35
## 2 40 14.6 42
## 3 45 23.7 48
## 4 50 32.4 55
## 5 55 39.3 64
## 6 60 43.9 71
## 7 65 48.7 75
## 8 70 54.7 86
## 9 75 61.5 92
## 10 80 68.7 106
## 11 85 75.2 116
## 12 90 80.1 124
## 13 95 85.4 129
## 14 100 92.4 138
## 15 105 95.2 140
## 16 110 98.3 152
## 17 115 99.4 156
## 18 120 99.7 158
## 19 125 100.0 170Table. 1. Conversion table of JMA LMIs equivalent to JTWC LMIs. LMIs come from JMA and JTWC best-track data sets, recorded at 5-kt intervals. JMA and JTWC use 10-min and 1-min wind averaging periods, respectively. 1-min equivalent JMA LMIs refer to the values at the same probability level of JTWC LMIs from interpolated quantiles. The values are for JJASON over the 30-year period 1986–2015.
set all variables (vectors)
# Sources
#--------------------------------------------------------------------------------
time=1986:2015
AB=read.table('Data/ABtrack2016RIeq.dat',T)
AB1=AB[(AB$bst=='JT' & AB$yr>=1984 & AB$mo>=6 & AB$mo<=11),] #1:JT
AB2=AB[(AB$bst=='JM' & AB$yr>=1984 & AB$mo>=6 & AB$mo<=11),] #2:JM
AB2$LMI=AB2$LMIeq # temporary conversion
GMSST=read.table('Data/gmsst_jjason_1854to2015.dat',T)[(time[1]-1853):(rev(time)[1]-1853),1]
NSOI=-1*read.table('Data/SOIann_jjason_1951to2015.dat',T)[(time[1]-1950):(rev(time)[1]-1950),1]
ABnri1=AB1[(AB1$RI==0),] ; ABnri2=AB2[(AB2$RI==0),] #non-RI
ABri1=AB1[(AB1$RI==1),] ; ABri2=AB2[(AB2$RI==1),] #RI
#JT total
#---------------
FRQ1=NULL ; for(i in time) FRQ1=c(FRQ1,length(AB1$LMI[AB1$yr==i]))
INT1=NULL ; for(i in time) INT1=c(INT1,mean(AB1$LMI[AB1$yr==i]))
dur1=NULL ; for(i in time) dur1=c(dur1,mean(AB1$dur[AB1$yr==i]))
AT1=NULL ; for(i in time) AT1=c(AT1,mean(AB1$AT[AB1$yr==i]))
#JM total
#---------------
FRQ2=NULL ; for(i in time) FRQ2=c(FRQ2,length(AB2$LMI[AB2$yr==i]))
INT2=NULL ; for(i in time) INT2=c(INT2,mean(AB2$LMI[AB2$yr==i]))
dur2=NULL ; for(i in time) dur2=c(dur2,mean(AB2$dur[AB2$yr==i]))
AT2=NULL ; for(i in time) AT2=c(AT2,mean(AB2$AT[AB2$yr==i]))
#JT RI/non-RI
#---------------
FRQri1=NULL ; for(i in time) FRQri1=c(FRQri1,length(ABri1$LMI[ABri1$yr==i]))
INTri1=NULL ; for(i in time) INTri1=c(INTri1,mean(ABri1$LMI[ABri1$yr==i]))
durri1=NULL ; for(i in time) durri1=c(durri1,mean(ABri1$dur[ABri1$yr==i]))
ATri1=NULL ; for(i in time) ATri1=c(ATri1,mean(ABri1$AT[ABri1$yr==i]))
#JM RI/non-RI
#---------------
FRQri2=NULL ; for(i in time) FRQri2=c(FRQri2,length(ABri2$LMI[ABri2$yr==i]))
INTri2=NULL ; for(i in time) INTri2=c(INTri2,mean(ABri2$LMI[ABri2$yr==i]))
durri2=NULL ; for(i in time) durri2=c(durri2,mean(ABri2$dur[ABri2$yr==i]))
ATri2=NULL ; for(i in time) ATri2=c(ATri2,mean(ABri2$AT[ABri2$yr==i]))
#MERGE
#---------------
FRQs=scale((scale(FRQ1)+scale(FRQ2))/sqrt(2))
FRQ=FRQs*sd(FRQ1)+mean(FRQ1)
INT=(scale(INT1)+scale(INT2))/sqrt(2)
RIp1=FRQri1/FRQ1 ; RIp2=FRQri2/FRQ2
RIps=scale((scale(RIp1)+scale(RIp2))/sqrt(2))
RIp=RIps*sd(RIp1)+mean(RIp1)
ATri=(scale(ATri1)+scale(ATri2))/sqrt(2)
ACT=(scale(INT)+scale(FRQ))/sqrt(2)
EINT=(scale(INT)-scale(FRQ))/sqrt(2)
ATg1=(AT1-ATri1)/AT1
ATg2=(AT2-ATri2)/AT2
ATg=(scale(ATg1)+scale(ATg2))/sqrt(2)Figure 1
AB=read.table('Data/ABtrack2016RIeq.dat',T)
ABjt=AB[(AB$bst=='JT' & AB$yr>=1986 & AB$mo>=6 & AB$mo<=11),]
ABjm=AB[(AB$bst=='JM' & AB$yr>=1986 & AB$mo>=6 & AB$mo<=11),]
ABjm$LMI=ABjm$LMIeq # temporary conversion
ABrijt=AB[(AB$bst=='JT' & AB$yr>=1986 & AB$mo>=6 & AB$mo<=11 & AB$RI==1),]
ABrijm=AB[(AB$bst=='JM' & AB$yr>=1986 & AB$mo>=6 & AB$mo<=11 & AB$RI==1),]
ABrijm$LMI=ABrijm$LMIeq # temporary conversion
#Setting
library(fields) #for image.plot## Warning: package 'fields' was built under R version 3.3.2## Loading required package: spam## Warning: package 'spam' was built under R version 3.3.2## Loading required package: grid## Spam version 1.4-0 (2016-08-29) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.##
## Attaching package: 'spam'## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve## Loading required package: maps## Warning: package 'maps' was built under R version 3.3.2library(MASS) # in case it is not already loaded
require(RColorBrewer) # for Brewer.pal --> color## Loading required package: RColorBrewer## Warning: package 'RColorBrewer' was built under R version 3.3.2Collist =c('white',brewer.pal(11, "RdYlBu")[7:11])
#brk=c(0,seq(0.00008,0.00024,0.00004),0.00035)
#lab.brk=c('','8','12','16','20','24','')
brk=c(0,seq(0.00004,0.00020,0.00004),0.00032)
lab.brk=c('','4','8','12','16','20','')
par(mfrow=c(2,2),mar=c(4,4,3.5,4))
#Plot1
A=ABjt$LMI ; B=ABjt$AT
Den = kde2d(A,B, n=400)
image(Den,breaks=brk,xlim=c(35,170),ylim=c(0,170),col=Collist, axes=F)
axis(1,at=seq(0,200,50),labels=seq(0,200,50),cex=1.2)
axis(2,at=seq(0,250,50),labels=seq(0,250,50),las=1,cex=1.2)
axis(3,at=seq(0,200,50),labels=rep('',length(seq(0,200,50))),cex=1.2)
axis(4,at=seq(0,250,50),labels=rep('',length(seq(0,250,50))),cex=1.2)
mtext('Life-time maximum intensity (kt)',side=1,line=2.4,cex=1)
mtext('Arrival time (hr)',side=2,line=2.6,cex=1)
mtext('(a)',side=3,line=1,cex=1.2,adj=-0.2)
image.plot(Den,breaks=brk,lab.breaks=lab.brk,col=Collist
,legend.args=list(text=expression(paste('x10'^-4)),side=3,line=0.4,cex=0.6)
,axis.args = list(cex.axis =0.8)
,legend.only=T)
#Plot2
A=ABjm$LMI ; B=ABjm$AT
Den = kde2d(A,B, n=400)
image(Den,breaks=brk,xlim=c(35,170),ylim=c(0,170),col=Collist, axes=F)
axis(1,at=seq(0,200,50),labels=seq(0,200,50),cex=1.2)
axis(2,at=seq(0,250,50),labels=seq(0,250,50),las=1,cex=1.2)
axis(3,at=seq(0,200,50),labels=rep('',length(seq(0,200,50))),cex=1.2)
axis(4,at=seq(0,250,50),labels=rep('',length(seq(0,250,50))),cex=1.2)
mtext('Life-time maximum intensity (kt)',side=1,line=2.4,cex=1)
mtext('Arrival time (hr)',side=2,line=2.6,cex=1)
mtext('(b)',side=3,line=1,cex=1.2,adj=-0.2)
image.plot(Den,breaks=brk,lab.breaks=lab.brk,col=Collist
,legend.args=list(text=expression(paste('x10'^-4)),side=3,line=0.4,cex=0.6)
,axis.args = list(cex.axis =0.8)
,legend.only=T)
#Plot3
A=ABrijt$LMI ; B=ABrijt$AT
Den = kde2d(A,B, n=400)
image(Den,breaks=brk,xlim=c(35,170),ylim=c(0,170),col=Collist, axes=F)
axis(1,at=seq(0,200,50),labels=seq(0,200,50),cex=1.2)
axis(2,at=seq(0,250,50),labels=seq(0,250,50),las=1,cex=1.2)
axis(3,at=seq(0,200,50),labels=rep('',length(seq(0,200,50))),cex=1.2)
axis(4,at=seq(0,250,50),labels=rep('',length(seq(0,250,50))),cex=1.2)
mtext('Life-time maximum intensity (kt)',side=1,line=2.4,cex=1)
mtext('Arrival time (hr)',side=2,line=2.6,cex=1)
mtext('(c)',side=3,line=1,cex=1.2,adj=-0.2)
image.plot(Den,breaks=brk,lab.breaks=lab.brk,col=Collist
,legend.args=list(text=expression(paste('x10'^-4)),side=3,line=0.4,cex=0.6)
,axis.args = list(cex.axis =0.8)
,legend.only=T)
rect(0, 0, 65, 200, col = 'light grey', border = NA)
box()
#Plot4
A=ABrijm$LMI ; B=ABrijm$AT
Den = kde2d(A,B, n=400)
image(Den,breaks=brk,xlim=c(35,170),ylim=c(0,170),col=Collist, axes=F)
axis(1,at=seq(0,200,50),labels=seq(0,200,50),cex=1.2)
axis(2,at=seq(0,250,50),labels=seq(0,250,50),las=1,cex=1.2)
axis(3,at=seq(0,200,50),labels=rep('',length(seq(0,200,50))),cex=1.2)
axis(4,at=seq(0,250,50),labels=rep('',length(seq(0,250,50))),cex=1.2)
mtext('Life-time maximum intensity (kt)',side=1,line=2.4,cex=1)
mtext('Arrival time (hr)',side=2,line=2.6,cex=1)
mtext('(d)',side=3,line=1,cex=1.2,adj=-0.2)
image.plot(Den,breaks=brk,lab.breaks=lab.brk,col=Collist
,legend.args=list(text=expression(paste('x10'^-4)),side=3,line=0.4,cex=0.6)
,axis.args = list(cex.axis =0.8)
,legend.only=T)
rect(0, 0, 75, 200, col = 'light grey', border = NA)
box()
Fig. 1. 2-dimensional density distributions of LMI (kt) and its arrival time (hr). (a) and (b) are for all TCs (LMI \(\geq\) 34 kt), and (c) and (d) for RI-experiencing TCs for JJASON over 30 years (1986–2015). (a) and (c) are from JTWC, and (b) and (d) from JMA best-track data. 1-min equivalent LMIs are used for JMA through Table 1.
Figure 2.
# (using mereged best track data) Totalnum, RInum, RIp
par(mfrow=c(1,1),mar=c(4,4,4,5))
plot(-10, -10, xlim=c(1986, 2015), ylim=c(0, 35), axes=FALSE, xlab='', ylab='')
#abline(h=0, lty=3, col='grey')
lines(time, FRQ, lwd=5, col='#00339999') # number of total TCs
model1=lm(FRQ~time)
abline(model1,lwd=2,col='#00339999')
30*summary(model1)$coef[2,1:2]## Estimate Std. Error
## -7.711982 2.387832lines(time, FRQ*RIp, lwd=5, col='#00339955') # number of RI-experiencing TCs
model2=lm(FRQ*RIp~time)
abline(model2,lwd=2,col='#00339955')
30*summary(model2)$coef[2,1:2]## Estimate Std. Error
## -0.09971252 1.36208686lines(time,35*RIp, lwd=3, col="#ff7f00") # proportion
model3=lm(RIp~time)
abline(lm(35*RIp~time),lwd=2,col="#ff7f00")
100*30*summary(model3)$coef[2,1:2] # rescaled(100%) over 30 years## Estimate Std. Error
## 15.910155 4.858283# axis
axis(1, at = seq(1985, 2015, 5), labels=seq(1985, 2015, 5), cex.axis=1)
axis(2, at = seq(0, 35, 5), labels=seq(0,35,5), las=1, cex.axis=1)
axis(4, at = 35*seq(0,1,0.1), labels=100*seq(0,1,0.1), las=1, cex.axis=1, col="#ff7f00", col.axis="#ff330099")
mtext('Year', side=1, line=2.6, cex=1.2)
mtext('Number', side=2, line=2.5, cex=1.2)
text(par("usr")[2] +4.1, 0.8, srt=-90, adj =0.94, labels = 'Proportion (%) of RI-experiencing TCs'
, cex=1.2, col="#ff7f00", xpd = TRUE)
legend(1997,35, bty='n',
c('Total TCs','RI-experiencing TCs', 'Proportion'),
text.col=c('#00339999','#00339955','#ff7f00'),
col=c('#00339999','#00339955','#ff7f00'),seg.len=3,
cex=0.8 ,lwd=c(5,5,3), ncol=1)
box()
Fig. 2. Time series of the number of RI-experiencing TCs, and the proportion (%) among the total TCs. The number of RI-experiencing TCs shows neutral trend (-0.10 \(\pm\) 1.36 (s.e.) /30 yr), while the proportion among the total TCs significantly increases (+15.9 \(\pm\) 4.86 (s.e.) %/30 yr). For comparison, annual variation of the total TCs is also presented, showing significantly decreasing trend (-7.7 \(\pm\) 2.39 (s.e.) /30 yr).
Figure 3.
#(using mereged best track data) RIp~GMSST+EINT, RIp~NSOI+EINT, RIp~non_AT+EINT
#RIp~EINT+GMSST
BEST=predict(lm(RIp~EINT+GMSST))
#A=acos(cor(BEST,EINT))*180/pi ; B=acos(cor(BEST,GMSST))*180/pi
r.sq=cor(RIp,BEST)^2 #variance portion of BEST
Va=r.sq*cor(BEST,GMSST)^2 #contribution of A
Vb=r.sq*cor(BEST,EINT)^2 #contribution of B
Int=Va+Vb-r.sq #Intersection
round(100*c(Va-Int, Int, Vb-Int),1)## [1] 3.7 47.7 16.1round(100*Va,1) #GMSST portion## [,1]
## [1,] 51.3round(100*Vb,1) #EINTportion## [,1]
## [1,] 63.8round(100*round(100*Int,1)/round(100*Va,1),1) #portion of EINT in GMSST influence## [,1]
## [1,] 93#RIp~EINT+NSOI
BEST=predict(lm(RIp~EINT+NSOI))
#A=acos(cor(BEST,EINT))*180/pi ; B=acos(cor(BEST,NSOI))*180/pi
r.sq=cor(RIp,BEST)^2 #variance portion of BEST
Va=r.sq*cor(BEST,NSOI)^2 #contribution of A
Vb=r.sq*cor(BEST,EINT)^2 #contribution of B
Int=Va+Vb-r.sq #Intersection
round(100*c(Va-Int, Int, Vb-Int),1)## [1] 4.2 11.9 51.9round(100*Va,1) #NSOI portion## [,1]
## [1,] 16.1round(100*Vb,1) #EINTportion## [,1]
## [1,] 63.8round(100*round(100*Int,1)/round(100*Va,1),1) #portion of EINT in NSOI influence## [,1]
## [1,] 73.9#RIp~EINT+ATg
BEST=predict(lm(RIp~EINT+ATg))
#A=acos(cor(BEST,EINT))*180/pi ; B=acos(cor(BEST,ATg))*180/pi
r.sq=cor(RIp,BEST)^2 #variance portion of BEST
Va=r.sq*cor(BEST,ATg)^2 #contribution of A
Vb=r.sq*cor(BEST,EINT)^2 #contribution of B
Int=Va+Vb-r.sq #Intersection
round(100*c(Va-Int, Int, Vb-Int),1)## [1] 11.1 29.4 34.4round(100*Va,1) #ATg portion## [,1]
## [1,] 40.5round(100*Vb,1) #EINTportion## [,1]
## [1,] 63.8round(100*round(100*Int,1)/round(100*Va,1),1) #portion of EINT in ATg influence## [,1]
## [1,] 72.6par(mar=c(0,0,0,0))
require(png)## Loading required package: png## Warning: package 'png' was built under R version 3.3.2require(grid)
img <- readPNG('Figure3_VariancePartitions.png')
grid.raster(img)
Fig. 3. Venn diagram for variance partitions of RIp. Contribution of (a) GMSST and EINT, (b) NSOI and EINT, and (c) ATg and EINT to the variance of RIp. The overlapping area means the variance portion how much the RIp variance can be explained. All values are rounded up to the first digit after the decimal point.
Figure 4.
library(RColorBrewer)
# 1. set option
#------------------------------------------------
Year=1986:2015
part1=data.frame(EINT=EINT,GMSST=GMSST,RIp=RIp)
# 2. read files
# (t, q, z, rh)
#------------------------------------------------
t=read.table('Data/temp.1984to2015.jjason.wnp.txt',FALSE)[-c(1:34),]# temperature (K)-273.15
q=read.table('Data/shum.1984to2015.jjason.wnp.txt',FALSE)[-c(1:34),]# specific humidity (g/kg)
z=read.table('Data/hgt.1984to2015.jjason.wnp.txt',FALSE)[-c(1:34),] # geopotential height (m)
# 3.calculate annual MSE (Moist Static Energy)
#------------------------------------------------
Cp= 1004 # specific heat capacity of air (J/(kg.K))
L= LHV=2.5*10^6 # latent heat of vaporization (J/kg)
g= 9.81 # accerelation due to gravity (m/s2)
SHann=Cp*(t+273.15) # sensible heating energy (J/kg)
LHann=L*q/1000 # latent heating (J/kg) # q over 300hPa is ignored to 0
PEann=g*z # potential energy (m2/s2)=(J/kg)
MSEann=SHann+PEann+LHann #annual Moist static energy
# 4.make each profile (30S~30N area-weighted mean values)
#------------------------------------------------
Lat=seq(-90,90,2.5)
n1=which(Lat==0) ; n2=which(Lat==30) # tropics
var=t
profileset=NULL #dim[time,lev]
for (lev in 1:17) # for each lev (1000 hPa to 10 hPa)
{tmsrs=NULL
annualnumset= seq(lev,length(Year)*17,17)[1:length(Year)]
for (ann in annualnumset) # time-series of area weight mean (0~30N)
{area_weighted_value=var[ann,(n1:n2)]*cos(Lat[n1:n2]*pi/180)/sum(cos(Lat[n1:n2]*pi/180))
tmsrs=c(tmsrs, sum(as.numeric(area_weighted_value)))}
profileset=rbind(profileset,tmsrs)
}
profileset_t=profileset
var=q
profileset=NULL #dim[time,lev]
for (lev in 1:8) # for each lev (1000 hPa to 10 hPa)
{tmsrs=NULL
annualnumset= seq(lev,length(Year)*17,17)[1:length(Year)]
for (ann in annualnumset) # time-series of area weight mean (0~30N)
{area_weighted_value=var[ann,(n1:n2)]*cos(Lat[n1:n2]*pi/180)/sum(cos(Lat[n1:n2]*pi/180))
tmsrs=c(tmsrs, sum(as.numeric(area_weighted_value)))}
profileset=rbind(profileset,tmsrs)
}
profileset_q=profileset
var=z
profileset=NULL #dim[time,lev]
for (lev in 1:17) # for each lev (1000 hPa to 10 hPa)
{tmsrs=NULL
annualnumset= seq(lev,length(Year)*17,17)[1:length(Year)]
for (ann in annualnumset) # time-series of area weight mean (0~30N)
{area_weighted_value=var[ann,(n1:n2)]*cos(Lat[n1:n2]*pi/180)/sum(cos(Lat[n1:n2]*pi/180))
tmsrs=c(tmsrs, sum(as.numeric(area_weighted_value)))}
profileset=rbind(profileset,tmsrs)
}
profileset_z=profileset
var=MSEann
profileset=NULL #dim[time,lev]
for (lev in 1:17) # for each lev (1000 hPa to 10 hPa)
{tmsrs=NULL
annualnumset= seq(lev,length(Year)*17,17)[1:length(Year)]
for (ann in annualnumset) # time-series of area weight mean (0~30N)
{area_weighted_value=var[ann,(n1:n2)]*cos(Lat[n1:n2]*pi/180)/sum(cos(Lat[n1:n2]*pi/180))
tmsrs=c(tmsrs, sum(as.numeric(area_weighted_value)))}
profileset=rbind(profileset,tmsrs)
}
profileset_MSE=profileset
# 5.calculate correlations
#------------------------------------------------
prof=profileset_z ; corr=NULL
for (lev in 1:17)
{part2=prof[lev,]
corr=rbind(corr,as.numeric(c(cor(part1[,1],part2),cor(part1[,2],part2),cor(part1[,3],part2))))}
corr_z=corr
prof=profileset_MSE ; corr=NULL
for (lev in 1:17)
{part2=prof[lev,]
corr=rbind(corr,as.numeric(c(cor(part1[,1],part2),cor(part1[,2],part2),cor(part1[,3],part2))))}
corr_MSE=corr
# 6. plot
#------------------------------------------------
p=log(c(1000,925,850,700,600,500,400,300,250,200,150,100,70,50,30,20,10))
p2=c(1000,925,850,700,600,500,400,300,250,200,150,100,70,50,30,20,10)
collist1=brewer.pal(n = 9, name = 'PuBu')[c(8,6,4)] #color for Geopotential height
collist2=brewer.pal(n = 9, name = 'OrRd')[c(8,6,4)] #color for MSE
par(mar=c(7,6,5,3))
plot(-999,-999,xlim=c(-1,1),ylim=c(log(1000),log(150)),xlab='',ylab='',axes=FALSE)
lines(corr_MSE[,3],p,col=collist2[3],lty=1,lwd=3) #EINT
lines(corr_MSE[,2],p,col=collist2[2],lty=1,lwd=3) #GMSST
lines(corr_MSE[,1],p,col=collist2[1],lty=1,lwd=7) #RIp
lines(corr_z[,3],p,col=collist1[3],lty=1,lwd=3) #EINT
lines(corr_z[,2],p,col=collist1[2],lty=1,lwd=3) #GMSST
lines(corr_z[,1],p,col=collist1[1],lty=1,lwd=7) #RIp
abline(v=0,col=1,lty=2)
axis(1,at=seq(-1,1,0.2),labels=seq(-1,1,0.2),cex.axis=1)
axis(2,at=p,labels=p2,cex.axis=1,las=1)
axis(3,at =seq(-1,1,0.2),labels=rep('',length=length(seq(-1,1,0.2))))
axis(4, at=p,labels=rep('',length(p)))
mtext('Correlation Coefficient',side=1,line=3, cex=1.4)
mtext('Pressure (hPa)',side=2,line=3.5, cex=1.4)
legend(-1.1,log(200), bty='n',
c('Geopotential Height by EINT','MSE by EINT',NA,
'Geopotential Height by GMSST','MSE by GMSST',NA,
'Geopotential Height by RIp','MSE by RIp'),
text.col=c(collist1[1],collist2[1],NA,collist1[2],collist2[2],NA,collist1[3],collist2[3]),
col=c(collist1[1],collist2[1],NA,collist1[2],collist2[2],NA,collist1[3],collist2[3]),seg.len=5,
cex=0.8 ,lty=c(1,1,NA,1,1,NA,1,1),lwd=c(7,7,NA,3,3,NA,3,3), ncol=1)
box()
Fig. 4. Correlation profiles of moist static energy (MSE) and geopotential height with GMSST, EINT, and RIp. Values are calculated over the tropics ((0\(^\circ\)–30\(^\circ\)N) in the western North Pacific (100\(^\circ\)E–180\(^\circ\)), for JJASON over 30 years (1986–2015). The larger correlation coefficients can be understood as the larger anomalies.