Narrowing region for tropical convections in the western North Pacific

Sanghyeon Yun, Namyoung Kang* and Chan Joo Jang

Abstract

Considering that the subtropical highs and tropical convections are observed as negative and positive vorticities respectively, the large-scale features of the atmospheric environment can be effectively represented using streamfunctions as defined by the Laplacian. By investigating the geographical patterns of streamfunctions from different modes of environmental variability, this study conceptualizes how the subtropical high expands and the region for tropical convections migrates in the western North Pacific. It is confirmed that, owing to the expansion of the subtropical high, the limited ocean area for tropical convections even bounded by the equator becomes narrower in the “La Niña mode” than that in the “El Niño mode”. This study finds that a warmer environment is likely to further expand the subtropical high to the west, and then the westernmost shift in the region for tropical convections appears in the “warmer La Niña mode”. A linear perspective suggests that every warmer La Niña environment could be one that people have scarcely experienced before.

Figure 1. Eight modes of environmental variability

library(png)
library(grid)
img = readPNG('Fig1.png')
grid.raster(img)

Fig. 1. Eight modes of environmental variability. A continuous variability space is formed by ENSO status (NSOI) and global ocean warmth (GMSST). PC1 and PC2 are the two principal components, which represent respectively the in-phase mode and out-of-phase mode of GMSST and NSOI. The eight modes of environmental variability are denoted as “Colder El Niño mode (CE)”, “El Niño mode (E)”, “Warmer El Niño mode (WE)”, “Colder mode (C)”, ’Warmer mode (W)”, “Colder La Niña mode (CL)”, “La Niña mode (L)”, and “Warmer La Niña mode (WL)” around the “Normal (N)” as each numbered from 1 to 8. The three circles indicate 0.5\(\sigma\), 1.0\(\sigma\), and 1.5\(\sigma\) outwards. CE and WE demonstrate the cold and warm anomalies around E, and therefore, 1 and 3 should have \(\sqrt{2}\) times larger values than E to be compared. The same for CL and WL as the anomalies around L, and 6 and 8 shows \(\sqrt{2}\) times larger values, too.

Figure 2. Funcionality of streamfunction

rm(list=ls())

intv=1
deg=seq(-180,180,intv)
wave=cos(deg*pi/180)
dwave=diff(wave) 
d2wave=diff(dwave) 

par(mar=c(6,8,3,3))
plot(-999,-999,xlim=c(-180,180),ylim=c(-1.1,1.4),axes=FALSE,xlab='',ylab='', xaxs='i')
abline(h=0,lty=2,col='grey')  ; #abline(v=0)
lines(deg,wave,col=1,lwd=4)
  deg=deg[1:(length(deg)-1)]+intv  ;  lines(deg,40*abs(dwave),col=4,lwd=2)
  deg=deg[1:(length(deg)-1)]+intv  ;  lines(deg,1500*d2wave,col=3,lwd=2)

axis(1, at = seq(-180,180,45),labels= c(rev(seq(0,180,45)),seq(45,180,45)), cex.axis=1.5)
axis(2, at = seq(-1,1,0.5), labels=c('-1.0','-0.5','0.0','0.5','1.0'), las=1, cex.axis=1.5)
mtext('Radial distance (conceptual)',side=1, line=4, cex=1.8)
mtext('Magnitude (conceptual)',side=2, line=5, cex=1.8)
legend( 'topleft', c(expression(paste( 'Streamfunction (',Psi,')')), expression(paste( 'Rotational Wind (','V'['r'],')')), expression(paste( 'Vorticity (',zeta,')')))
         ,cex=1.3, lwd=c(4,2,2), col=c(1,4,3), x.intersp=1, y.intersp=1)

box()

Fig. 2. Funcionality of streamfunction. The schematic diagram shows an interrelationship among streamfunction (\(\psi\)), rotational wind (\(V_{r}\)), and vorticity (\(\zeta\)) using a one-dimensional distribution of \(\psi\). \(\psi\)’s distribution follows a sine curve for an example. \(V_{r}\) is defined by the gradient of \(\psi\), and \(\zeta\) is formed by the gradient of \(V_{r}\). From the Laplacian, \(\psi\) is negatively proportional to \(\zeta\), and represents a large-scale pattern of \(\zeta\). The value on the x-axis conceptually represents the radial distance from the \(\psi\)’s maximum. For the purpose of comparison in the same plot, \(V_{r}\) and \(\zeta\) are 40 and 1500 times multiplied, respectively.

Figure 3. Geographical patterns of streamfunction

rm(list=ls())

year = 1985:2020 ; level = 850  # 200 / 500
laplace = "psi" # chi / psi
extraction = "r"    # d -> divergent wind / r -> rotational wind
fin = "vorticity"   # divergence / vorticity
#phase = "La-Niña"  # El-Niño / La-Niña / Normal

############################################
library(ncdf4)

## Warning: package 'ncdf4' was built under R version 4.1.1

ncfname_lap = paste(getwd(),"/",laplace,level,"hPa_JJASON_",year[1],"to",year[length(year)],".nc",sep="")

dname_lap = laplace

# 1. Lap. Regression
# Latitude of Chi&Psi are set to "degrees_north" (-90 ~ 90)
baselon = ncvar_get(nc_open(ncfname_lap),"lon") ; nlon = length(baselon)
baselat = ncvar_get(nc_open(ncfname_lap),"lat") ; nlat = length(baselat)
lon=rep(seq(0,357.5,2.5),73) ; lat=rep(seq(-90,90,2.5),each=144)
gmsst=read.table('gmsst_jjason_1854to2020.dat',T)[(year[1]-1853):(rev(year)[1]-1853),1]
nsoi=-1*read.table('soi_jjason_1951to2020.dat',T)[(year[1]-1950):(rev(year)[1]-1950),1]
############################################
lap = matrix(ncvar_get(nc_open(ncfname_lap),dname_lap),nrow=nlon*nlat)
lap = cbind(lon,lat,lap)
lap = lap[c(order(lap[,1],-lap[,2])),]

PC1 = (scale(gmsst)[,1] + scale(nsoi)[,1])/sqrt(2)
PC2 = (scale(gmsst)[,1] - scale(nsoi)[,1])/sqrt(2)

# El-Niño Phase
lap_ColdEl_pred = cbind(lap[,1],lap[,2],rep(NA,length(lon)))    # -PC2
lap_ColdEl_pval = cbind(lap[,1],lap[,2],rep(NA,length(lon)))    # -PC2
lap_El_pred = cbind(lap[,1],lap[,2],rep(NA,length(lon)))        # NSOI
lap_El_pval = cbind(lap[,1],lap[,2],rep(NA,length(lon)))        # NSOI
lap_WarmEl_pred = cbind(lap[,1],lap[,2],rep(NA,length(lon)))        # PC1
lap_WarmEl_pval = cbind(lap[,1],lap[,2],rep(NA,length(lon)))        # PC1

# Normal Phase
lap_Cold_pred = cbind(lap[,1],lap[,2],rep(NA,length(lon)))  # -GMSST
lap_Cold_pval = cbind(lap[,1],lap[,2],rep(NA,length(lon)))  # -GMSST
lap_Mean = cbind(lap[,1],lap[,2],rep(NA,length(lon)))   # Mean
lap_Warm_pred = cbind(lap[,1],lap[,2],rep(NA,length(lon)))  # GMSST
lap_Warm_pval = cbind(lap[,1],lap[,2],rep(NA,length(lon)))  # GMSST

# La-Niña Phase
lap_ColdLa_pred = cbind(lap[,1],lap[,2],rep(NA,length(lon)))    # -PC1
lap_ColdLa_pval = cbind(lap[,1],lap[,2],rep(NA,length(lon)))    # -PC1
lap_La_pred = cbind(lap[,1],lap[,2],rep(NA,length(lon)))    # -NSOI
lap_La_pval = cbind(lap[,1],lap[,2],rep(NA,length(lon)))    # -NSOI
lap_WarmLa_pred = cbind(lap[,1],lap[,2],rep(NA,length(lon)))    # PC2
lap_WarmLa_pval = cbind(lap[,1],lap[,2],rep(NA,length(lon)))    # PC2

# -PC2 -> NSOI -> PC1 -> -GMSST -> Mean -> GMSST -> -PC1 -> -NSOI -> PC2
# ColdEL -> EL -> WarmEL -> Cold -> Mean -> Warm -> ColdLA -> LA -> WarmLa

reglap = lap[,3:(length(year)+2)]
for(i in 1:length(lon))
{
    y = reglap[i,]
    #y = scale(reglap[i,])[,1]
    
    # 1. Cold pattern El-Niño phase (-PC2)
    x1 = -1*PC2
    model1 = lm(y~x1)
    lap_ColdEl_pred[i,3] = predict(model1,data.frame(x1=(1.5*sqrt(2)))) # 2.1-sigma
    lap_ColdEl_pval[i,3] = summary(model1)$coef[2,4]
    
    # 2. Normal pattern El-Niño phase (nsoi)
    x2 = scale(nsoi)[,1]
    model2 = lm(y~x2)
    lap_El_pred[i,3] = predict(model2,data.frame(x2=1.5))   # 1.5-sigma
    lap_El_pval[i,3] = summary(model2)$coef[2,4]
    
    # 3. Warm pattern El-Niño phase (PC1)
    x3 = PC1
    model3 = lm(y~x3)
    lap_WarmEl_pred[i,3] = predict(model3,data.frame(x3=(1.5*sqrt(2)))) # 2.1-sigma
    lap_WarmEl_pval[i,3] = summary(model3)$coef[2,4]
    
    # 4. Cold pattern Normal phase (-gmsst)
    x4 = scale(-1*gmsst)[,1]
    model4 = lm(y~x4)
    lap_Cold_pred[i,3] = predict(model4,data.frame(x4=1.5)) # 1.5-sigma
    lap_Cold_pval[i,3] = summary(model4)$coef[2,4]
    
    # 5. Normal phase (mean)
    lap_Mean[i,3] = mean(lap[i,3:(length(year)+2)])
    
    # 6. Warm pattern Normal phase (gmsst)
    x6 = scale(gmsst)[,1]
    model6 = lm(y~x6)
    lap_Warm_pred[i,3] = predict(model6,data.frame(x6=1.5)) # 1.5-sigma
    lap_Warm_pval[i,3] = summary(model6)$coef[2,4]
    
    # 7. Cold pattern La-Niña phase (-PC1)
    x7 = -1*PC1
    model7 = lm(y~x7)
    lap_ColdLa_pred[i,3] = predict(model7,data.frame(x7=(1.5*sqrt(2)))) # 2.1-sigma
    lap_ColdLa_pval[i,3] = summary(model7)$coef[2,4]
    
    # 8. Normal pattern La-Niña phase (-nsoi)
    x8 = scale(-1*nsoi)[,1]
    model8 = lm(y~x8)
    lap_La_pred[i,3] = predict(model8,data.frame(x8=1.5))   # 1.5-sigma
    lap_La_pval[i,3] = summary(model8)$coef[2,4]
    
    # 9. Warm pattern La-Niña phase (pc2)
    x9 = PC2
    model9 = lm(y~x9)
    lap_WarmLa_pred[i,3] = predict(model9,data.frame(x9=(1.5*sqrt(2)))) # 2.1-sigma
    lap_WarmLa_pval[i,3] = summary(model9)$coef[2,4]    
}

# El-Niño Phase
colnames(lap_ColdEl_pred)=c("lon","lat","pred")
colnames(lap_ColdEl_pval)=c("lon","lat","pval")
colnames(lap_El_pred)=c("lon","lat","pred")
colnames(lap_El_pval)=c("lon","lat","pval")
colnames(lap_WarmEl_pred)=c("lon","lat","pred")
colnames(lap_WarmEl_pval)=c("lon","lat","pval")

# Normal Phase
colnames(lap_Cold_pred)=c("lon","lat","pred")
colnames(lap_Cold_pval)=c("lon","lat","pval")
colnames(lap_Mean)=c("lon","lat","mean")
colnames(lap_Warm_pred)=c("lon","lat","pred")
colnames(lap_Warm_pval)=c("lon","lat","pred")

# La-Niña Phase
colnames(lap_ColdLa_pred)=c("lon","lat","pred")
colnames(lap_ColdLa_pval)=c("lon","lat","pval")
colnames(lap_La_pred)=c("lon","lat","pred")
colnames(lap_La_pval)=c("lon","lat","pval")
colnames(lap_WarmLa_pred)=c("lon","lat","pred")
colnames(lap_WarmLa_pval)=c("lon","lat","pval")


# El-Niño Phase
lap_ColdEl_pred = lap_ColdEl_pred[c(order(-lap_ColdEl_pred[,2])),]
lap_ColdEl_pval = lap_ColdEl_pval[c(order(-lap_ColdEl_pval[,2])),]
lap_El_pred = lap_El_pred[c(order(-lap_El_pred[,2])),]
lap_El_pval = lap_El_pval[c(order(-lap_El_pval[,2])),]
lap_WarmEl_pred = lap_WarmEl_pred[c(order(-lap_WarmEl_pred[,2])),]
lap_WarmEl_pval = lap_WarmEl_pval[c(order(-lap_WarmEl_pval[,2])),]

# Normal Phase
lap_Cold_pred = lap_Cold_pred[c(order(-lap_Cold_pred[,2])),]
lap_Cold_pval = lap_Cold_pval[c(order(-lap_Cold_pval[,2])),]
lap_Mean = lap_Mean[c(order(-lap_Mean[,2])),]
lap_Warm_pred = lap_Warm_pred[c(order(-lap_Warm_pred[,2])),]
lap_Warm_pval = lap_Warm_pval[c(order(-lap_Warm_pval[,2])),]

# La-Niña Phase
lap_ColdLa_pred = lap_ColdLa_pred[c(order(-lap_ColdLa_pred[,2])),]
lap_ColdLa_pval = lap_ColdLa_pval[c(order(-lap_ColdLa_pval[,2])),]
lap_La_pred = lap_La_pred[c(order(-lap_La_pred[,2])),]
lap_La_pval = lap_La_pval[c(order(-lap_La_pval[,2])),]
lap_WarmLa_pred = lap_WarmLa_pred[c(order(-lap_WarmLa_pred[,2])),]
lap_WarmLa_pval = lap_WarmLa_pval[c(order(-lap_WarmLa_pval[,2])),]

############################################
############################################
library(geosphere)

## Warning: package 'geosphere' was built under R version 4.1.1

# Cal. partial x
p1=cbind(rep(0,73),seq(90,-90,-2.5)) ; p2=cbind(rep(2.5,73),seq(90,-90,-2.5))
dx=distVincentyEllipsoid(p1,p2,a=6378137, b=6356752.3142, f=1/298.257223563)

# Cal. partial y
imsi=c(90,seq(90,-90,-2.5),-90) ; n=length(imsi)
p1=cbind(rep(0,(n-2)),imsi[1:(n-2)]) ; p2= cbind(rep(0,(n-2)),imsi[3:n])
dy=(1/2)*distVincentyEllipsoid(p1,p2,a=6378137, b=6356752.3142, f=1/298.257223563)

# 1-1. Calculate New Divergent/Rotational Wind from Lap~ColdEl.
imsiU = NULL ; imsiV = NULL
imsiU = lap_ColdEl_pred ; imsiV = lap_ColdEl_pred
imsiU = imsiU[c(order(imsiU[,1])),] ; imsiV = imsiV[c(order(imsiV[,1])),]

U = as.matrix(imsiU[,3]) ; V = as.matrix(imsiV[,3])

Lon = seq(0,357.5,2.5) ; ilen=length(Lon)
Lat = seq(90,-90,-2.5) ; jlen=length(Lat)

dUdx=NULL ; dVdx=NULL ; dUdy=NULL ; dVdy=NULL
for(k in 1)   #year 
{
    dUdx_ann=NULL ; dVdx_ann=NULL ; dUdy_ann=NULL ; dVdy_ann=NULL
    U_ann=matrix(U[,k],ncol=ilen) ; V_ann=matrix(V[,k],ncol=ilen)
    for(j in 1:jlen) {d=diff(c(U_ann[j,ilen], U_ann[j,], U_ann[j,1])) ; dUdx_ann=c(dUdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
    for(i in 1:ilen) {d=diff(c(U_ann[1,i], U_ann[,i], U_ann[jlen,i])) ; dUdy_ann=c(dUdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
    for(j in 1:jlen) {d=diff(c(V_ann[j,ilen], V_ann[j,], V_ann[j,1])) ; dVdx_ann=c(dVdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
    for(i in 1:ilen) {d=diff(c(V_ann[1,i], V_ann[,i], V_ann[jlen,i])) ; dVdy_ann=c(dVdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
    
    dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)
}

#dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)

dUdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dUdx)
dVdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dVdx)
dUdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dUdy)
dVdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dVdy)

dUdy = dUdy[c(order(-dUdy[,2])),] ; dVdy = dVdy[c(order(-dVdy[,2])),]
dUdy[,3] = -1*dUdy[,3]
dVdy[,3] = -1*dVdy[,3]

dUdy[,3] = -1*dUdy[,3]

ud_ColdEl = dUdx ; vd_ColdEl = dVdy
ur_ColdEl = dUdy ; vr_ColdEl = dVdx

rm(imsiU,imsiV,U,V,dUdx,dUdx_ann,dVdx,dVdx_ann,dUdy,dUdy_ann,dVdy,dVdy_ann)

############################################
# 1-2. Calculate New Divergent/Rotational Wind from Lap~El.
imsiU = NULL ; imsiV = NULL
imsiU = lap_El_pred ; imsiV = lap_El_pred
imsiU = imsiU[c(order(imsiU[,1])),] ; imsiV = imsiV[c(order(imsiV[,1])),]

U = as.matrix(imsiU[,3]) ; V = as.matrix(imsiV[,3])

Lon = seq(0,357.5,2.5) ; ilen=length(Lon)
Lat = seq(90,-90,-2.5) ; jlen=length(Lat)

dUdx=NULL ; dVdx=NULL ; dUdy=NULL ; dVdy=NULL
dUdx_ann=NULL ; dVdx_ann=NULL ; dUdy_ann=NULL ; dVdy_ann=NULL
U_ann=matrix(U[,1],ncol=ilen) ; V_ann=matrix(V[,1],ncol=ilen)
for(j in 1:jlen) {d=diff(c(U_ann[j,ilen], U_ann[j,], U_ann[j,1])) ; dUdx_ann=c(dUdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(U_ann[1,i], U_ann[,i], U_ann[jlen,i])) ; dUdy_ann=c(dUdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
for(j in 1:jlen) {d=diff(c(V_ann[j,ilen], V_ann[j,], V_ann[j,1])) ; dVdx_ann=c(dVdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(V_ann[1,i], V_ann[,i], V_ann[jlen,i])) ; dVdy_ann=c(dVdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}

dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)

dUdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dUdx)
dVdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dVdx)
dUdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dUdy)
dVdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dVdy)

dUdy = dUdy[c(order(-dUdy[,2])),] ; dVdy = dVdy[c(order(-dVdy[,2])),]
dUdy[,3] = -1*dUdy[,3]
dVdy[,3] = -1*dVdy[,3]

dUdy[,3] = -1*dUdy[,3]

ud_El = dUdx ; vd_El = dVdy
ur_El = dUdy ; vr_El = dVdx

rm(imsiU,imsiV,U,V,dUdx,dUdx_ann,dVdx,dVdx_ann,dUdy,dUdy_ann,dVdy,dVdy_ann)

############################################
# 1-3. Calculate New Divergent/Rotational Wind from Lap~WarmEl.
imsiU = NULL ; imsiV = NULL
imsiU = lap_WarmEl_pred ; imsiV = lap_WarmEl_pred
imsiU = imsiU[c(order(imsiU[,1])),] ; imsiV = imsiV[c(order(imsiV[,1])),]

U = as.matrix(imsiU[,3]) ; V = as.matrix(imsiV[,3])

Lon = seq(0,357.5,2.5) ; ilen=length(Lon)
Lat = seq(90,-90,-2.5) ; jlen=length(Lat)

dUdx=NULL ; dVdx=NULL ; dUdy=NULL ; dVdy=NULL
for(k in 1)   #year 
{
    dUdx_ann=NULL ; dVdx_ann=NULL ; dUdy_ann=NULL ; dVdy_ann=NULL
    U_ann=matrix(U[,k],ncol=ilen) ; V_ann=matrix(V[,k],ncol=ilen)
    for(j in 1:jlen) {d=diff(c(U_ann[j,ilen], U_ann[j,], U_ann[j,1])) ; dUdx_ann=c(dUdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
    for(i in 1:ilen) {d=diff(c(U_ann[1,i], U_ann[,i], U_ann[jlen,i])) ; dUdy_ann=c(dUdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
    for(j in 1:jlen) {d=diff(c(V_ann[j,ilen], V_ann[j,], V_ann[j,1])) ; dVdx_ann=c(dVdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
    for(i in 1:ilen) {d=diff(c(V_ann[1,i], V_ann[,i], V_ann[jlen,i])) ; dVdy_ann=c(dVdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
    
    dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)
}

#dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)

dUdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dUdx)
dVdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dVdx)
dUdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dUdy)
dVdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dVdy)

dUdy = dUdy[c(order(-dUdy[,2])),] ; dVdy = dVdy[c(order(-dVdy[,2])),]
dUdy[,3] = -1*dUdy[,3]
dVdy[,3] = -1*dVdy[,3]

dUdy[,3] = -1*dUdy[,3]

ud_WarmEl = dUdx ; vd_WarmEl = dVdy
ur_WarmEl = dUdy ; vr_WarmEl = dVdx

rm(imsiU,imsiV,U,V,dUdx,dUdx_ann,dVdx,dVdx_ann,dUdy,dUdy_ann,dVdy,dVdy_ann)

############################################
# 1-4. Calculate New Divergent/Rotational Wind from Lap~Cold.
imsiU = NULL ; imsiV = NULL
imsiU = lap_Cold_pred ; imsiV = lap_Cold_pred
imsiU = imsiU[c(order(imsiU[,1])),] ; imsiV = imsiV[c(order(imsiV[,1])),]

U = as.matrix(imsiU[,3]) ; V = as.matrix(imsiV[,3])

Lon = seq(0,357.5,2.5) ; ilen=length(Lon)
Lat = seq(90,-90,-2.5) ; jlen=length(Lat)

dUdx=NULL ; dVdx=NULL ; dUdy=NULL ; dVdy=NULL
for(k in 1)   #year 
{
    dUdx_ann=NULL ; dVdx_ann=NULL ; dUdy_ann=NULL ; dVdy_ann=NULL
    U_ann=matrix(U[,k],ncol=ilen) ; V_ann=matrix(V[,k],ncol=ilen)
    for(j in 1:jlen) {d=diff(c(U_ann[j,ilen], U_ann[j,], U_ann[j,1])) ; dUdx_ann=c(dUdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
    for(i in 1:ilen) {d=diff(c(U_ann[1,i], U_ann[,i], U_ann[jlen,i])) ; dUdy_ann=c(dUdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
    for(j in 1:jlen) {d=diff(c(V_ann[j,ilen], V_ann[j,], V_ann[j,1])) ; dVdx_ann=c(dVdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
    for(i in 1:ilen) {d=diff(c(V_ann[1,i], V_ann[,i], V_ann[jlen,i])) ; dVdy_ann=c(dVdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
    
    dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)
}

#dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)

dUdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dUdx)
dVdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dVdx)
dUdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dUdy)
dVdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dVdy)

dUdy = dUdy[c(order(-dUdy[,2])),] ; dVdy = dVdy[c(order(-dVdy[,2])),]
dUdy[,3] = -1*dUdy[,3]
dVdy[,3] = -1*dVdy[,3]

dUdy[,3] = -1*dUdy[,3]

ud_Cold = dUdx ; vd_Cold = dVdy
ur_Cold = dUdy ; vr_Cold = dVdx

rm(imsiU,imsiV,U,V,dUdx,dUdx_ann,dVdx,dVdx_ann,dUdy,dUdy_ann,dVdy,dVdy_ann)

############################################
# 1-5. Calculate New Divergent/Rotational Wind from Lap~Mean.
imsiU = NULL ; imsiV = NULL
imsiU = lap_Mean ; imsiV = lap_Mean
imsiU = imsiU[c(order(imsiU[,1])),] ; imsiV = imsiV[c(order(imsiV[,1])),]

U = as.matrix(imsiU[,3]) ; V = as.matrix(imsiV[,3])

Lon = seq(0,357.5,2.5) ; ilen=length(Lon)
Lat = seq(90,-90,-2.5) ; jlen=length(Lat)

dUdx=NULL ; dVdx=NULL ; dUdy=NULL ; dVdy=NULL
for(k in 1)   #year 
{
    dUdx_ann=NULL ; dVdx_ann=NULL ; dUdy_ann=NULL ; dVdy_ann=NULL
    U_ann=matrix(U[,k],ncol=ilen) ; V_ann=matrix(V[,k],ncol=ilen)
    for(j in 1:jlen) {d=diff(c(U_ann[j,ilen], U_ann[j,], U_ann[j,1])) ; dUdx_ann=c(dUdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
    for(i in 1:ilen) {d=diff(c(U_ann[1,i], U_ann[,i], U_ann[jlen,i])) ; dUdy_ann=c(dUdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
    for(j in 1:jlen) {d=diff(c(V_ann[j,ilen], V_ann[j,], V_ann[j,1])) ; dVdx_ann=c(dVdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
    for(i in 1:ilen) {d=diff(c(V_ann[1,i], V_ann[,i], V_ann[jlen,i])) ; dVdy_ann=c(dVdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
    
    dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)
}

#dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)

dUdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dUdx)
dVdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dVdx)
dUdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dUdy)
dVdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dVdy)

dUdy = dUdy[c(order(-dUdy[,2])),] ; dVdy = dVdy[c(order(-dVdy[,2])),]
dUdy[,3] = -1*dUdy[,3]
dVdy[,3] = -1*dVdy[,3]

dUdy[,3] = -1*dUdy[,3]

ud_Mean = dUdx ; vd_Mean = dVdy
ur_Mean = dUdy ; vr_Mean = dVdx

rm(imsiU,imsiV,U,V,dUdx,dUdx_ann,dVdx,dVdx_ann,dUdy,dUdy_ann,dVdy,dVdy_ann)

############################################
# 1-6. Calculate New Divergent/Rotational Wind from Lap~Warm.
imsiU = NULL ; imsiV = NULL
imsiU = lap_Warm_pred ; imsiV = lap_Warm_pred
imsiU = imsiU[c(order(imsiU[,1])),] ; imsiV = imsiV[c(order(imsiV[,1])),]

U = as.matrix(imsiU[,3]) ; V = as.matrix(imsiV[,3])

Lon = seq(0,357.5,2.5) ; ilen=length(Lon)
Lat = seq(90,-90,-2.5) ; jlen=length(Lat)

dUdx=NULL ; dVdx=NULL ; dUdy=NULL ; dVdy=NULL
for(k in 1)   #year 
{
    dUdx_ann=NULL ; dVdx_ann=NULL ; dUdy_ann=NULL ; dVdy_ann=NULL
    U_ann=matrix(U[,k],ncol=ilen) ; V_ann=matrix(V[,k],ncol=ilen)
    for(j in 1:jlen) {d=diff(c(U_ann[j,ilen], U_ann[j,], U_ann[j,1])) ; dUdx_ann=c(dUdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
    for(i in 1:ilen) {d=diff(c(U_ann[1,i], U_ann[,i], U_ann[jlen,i])) ; dUdy_ann=c(dUdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
    for(j in 1:jlen) {d=diff(c(V_ann[j,ilen], V_ann[j,], V_ann[j,1])) ; dVdx_ann=c(dVdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
    for(i in 1:ilen) {d=diff(c(V_ann[1,i], V_ann[,i], V_ann[jlen,i])) ; dVdy_ann=c(dVdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
    
    dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)
}

#dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)

dUdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dUdx)
dVdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dVdx)
dUdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dUdy)
dVdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dVdy)

dUdy = dUdy[c(order(-dUdy[,2])),] ; dVdy = dVdy[c(order(-dVdy[,2])),]
dUdy[,3] = -1*dUdy[,3]
dVdy[,3] = -1*dVdy[,3]

dUdy[,3] = -1*dUdy[,3]

ud_Warm = dUdx ; vd_Warm = dVdy
ur_Warm = dUdy ; vr_Warm = dVdx

rm(imsiU,imsiV,U,V,dUdx,dUdx_ann,dVdx,dVdx_ann,dUdy,dUdy_ann,dVdy,dVdy_ann)

############################################
# 1-7. Calculate New Divergent/Rotational Wind from Lap~ColdLa.
imsiU = NULL ; imsiV = NULL
imsiU = lap_ColdLa_pred ; imsiV = lap_ColdLa_pred
imsiU = imsiU[c(order(imsiU[,1])),] ; imsiV = imsiV[c(order(imsiV[,1])),]

U = as.matrix(imsiU[,3]) ; V = as.matrix(imsiV[,3])

Lon = seq(0,357.5,2.5) ; ilen=length(Lon)
Lat = seq(90,-90,-2.5) ; jlen=length(Lat)

dUdx=NULL ; dVdx=NULL ; dUdy=NULL ; dVdy=NULL
for(k in 1)   #year 
{
    dUdx_ann=NULL ; dVdx_ann=NULL ; dUdy_ann=NULL ; dVdy_ann=NULL
    U_ann=matrix(U[,k],ncol=ilen) ; V_ann=matrix(V[,k],ncol=ilen)
    for(j in 1:jlen) {d=diff(c(U_ann[j,ilen], U_ann[j,], U_ann[j,1])) ; dUdx_ann=c(dUdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
    for(i in 1:ilen) {d=diff(c(U_ann[1,i], U_ann[,i], U_ann[jlen,i])) ; dUdy_ann=c(dUdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
    for(j in 1:jlen) {d=diff(c(V_ann[j,ilen], V_ann[j,], V_ann[j,1])) ; dVdx_ann=c(dVdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
    for(i in 1:ilen) {d=diff(c(V_ann[1,i], V_ann[,i], V_ann[jlen,i])) ; dVdy_ann=c(dVdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
    
    dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)
}

#dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)

dUdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dUdx)
dVdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dVdx)
dUdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dUdy)
dVdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dVdy)

dUdy = dUdy[c(order(-dUdy[,2])),] ; dVdy = dVdy[c(order(-dVdy[,2])),]
dUdy[,3] = -1*dUdy[,3]
dVdy[,3] = -1*dVdy[,3]

dUdy[,3] = -1*dUdy[,3]

ud_ColdLa = dUdx ; vd_ColdLa = dVdy
ur_ColdLa = dUdy ; vr_ColdLa = dVdx

rm(imsiU,imsiV,U,V,dUdx,dUdx_ann,dVdx,dVdx_ann,dUdy,dUdy_ann,dVdy,dVdy_ann)

############################################
# 1-8. Calculate New Divergent/Rotational Wind from Lap~La.
imsiU = NULL ; imsiV = NULL
imsiU = lap_La_pred ; imsiV = lap_La_pred
imsiU = imsiU[c(order(imsiU[,1])),] ; imsiV = imsiV[c(order(imsiV[,1])),]

U = as.matrix(imsiU[,3]) ; V = as.matrix(imsiV[,3])

Lon = seq(0,357.5,2.5) ; ilen=length(Lon)
Lat = seq(90,-90,-2.5) ; jlen=length(Lat)

dUdx=NULL ; dVdx=NULL ; dUdy=NULL ; dVdy=NULL
dUdx_ann=NULL ; dVdx_ann=NULL ; dUdy_ann=NULL ; dVdy_ann=NULL
U_ann=matrix(U[,1],ncol=ilen) ; V_ann=matrix(V[,1],ncol=ilen)
for(j in 1:jlen) {d=diff(c(U_ann[j,ilen], U_ann[j,], U_ann[j,1])) ; dUdx_ann=c(dUdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(U_ann[1,i], U_ann[,i], U_ann[jlen,i])) ; dUdy_ann=c(dUdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
for(j in 1:jlen) {d=diff(c(V_ann[j,ilen], V_ann[j,], V_ann[j,1])) ; dVdx_ann=c(dVdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(V_ann[1,i], V_ann[,i], V_ann[jlen,i])) ; dVdy_ann=c(dVdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}

dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)

dUdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dUdx)
dVdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dVdx)
dUdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dUdy)
dVdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dVdy)

dUdy = dUdy[c(order(-dUdy[,2])),] ; dVdy = dVdy[c(order(-dVdy[,2])),]
dUdy[,3] = -1*dUdy[,3]
dVdy[,3] = -1*dVdy[,3]

dUdy[,3] = -1*dUdy[,3]

ud_La = dUdx ; vd_La = dVdy
ur_La = dUdy ; vr_La = dVdx

rm(imsiU,imsiV,U,V,dUdx,dUdx_ann,dVdx,dVdx_ann,dUdy,dUdy_ann,dVdy,dVdy_ann)

############################################
# 1-9. Calculate New Divergent/Rotational Wind from Lap~WarmLa.
imsiU = NULL ; imsiV = NULL
imsiU = lap_WarmLa_pred ; imsiV = lap_WarmLa_pred
imsiU = imsiU[c(order(imsiU[,1])),] ; imsiV = imsiV[c(order(imsiV[,1])),]

U = as.matrix(imsiU[,3]) ; V = as.matrix(imsiV[,3])

Lon = seq(0,357.5,2.5) ; ilen=length(Lon)
Lat = seq(90,-90,-2.5) ; jlen=length(Lat)

dUdx=NULL ; dVdx=NULL ; dUdy=NULL ; dVdy=NULL
for(k in 1)   #year 
{
    dUdx_ann=NULL ; dVdx_ann=NULL ; dUdy_ann=NULL ; dVdy_ann=NULL
    U_ann=matrix(U[,k],ncol=ilen) ; V_ann=matrix(V[,k],ncol=ilen)
    for(j in 1:jlen) {d=diff(c(U_ann[j,ilen], U_ann[j,], U_ann[j,1])) ; dUdx_ann=c(dUdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
    for(i in 1:ilen) {d=diff(c(U_ann[1,i], U_ann[,i], U_ann[jlen,i])) ; dUdy_ann=c(dUdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
    for(j in 1:jlen) {d=diff(c(V_ann[j,ilen], V_ann[j,], V_ann[j,1])) ; dVdx_ann=c(dVdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
    for(i in 1:ilen) {d=diff(c(V_ann[1,i], V_ann[,i], V_ann[jlen,i])) ; dVdy_ann=c(dVdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
    
    dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)
}

#dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)

dUdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dUdx)
dVdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dVdx)
dUdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dUdy)
dVdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dVdy)

dUdy = dUdy[c(order(-dUdy[,2])),] ; dVdy = dVdy[c(order(-dVdy[,2])),]
dUdy[,3] = -1*dUdy[,3]
dVdy[,3] = -1*dVdy[,3]

dUdy[,3] = -1*dUdy[,3]

ud_WarmLa = dUdx ; vd_WarmLa = dVdy
ur_WarmLa = dUdy ; vr_WarmLa = dVdx

rm(imsiU,imsiV,U,V,dUdx,dUdx_ann,dVdx,dVdx_ann,dUdy,dUdy_ann,dVdy,dVdy_ann)

############################################

exlon=rep(seq(0,357.5,2.5),73) ; exlat=rep(seq(90,-90,-2.5),each=144)
nexlon = length(exlon) ; nexlat = length(exlat)

############################################
############################################
# 2-1. Scalar Wind (Lap~ColdEl)
exu_ColdEl = get(paste("u",extraction,"_ColdEl",sep=""))
exv_ColdEl = get(paste("v",extraction,"_ColdEl",sep=""))

exu_ColdEl = exu_ColdEl[,3] ; exv_ColdEl = exv_ColdEl[,3]

exwind_ColdEl = cbind(exlon,exlat,rep(NA,length(exlon)))
for(i in 1:nexlon)
{
    exwind_ColdEl[i,3] = sqrt((exu_ColdEl[i]^2)+(exv_ColdEl[i]^2))
}

############################################
# 2-2. Scalar Wind (Lap~El)
exu_El = get(paste("u",extraction,"_El",sep=""))
exv_El = get(paste("v",extraction,"_El",sep=""))

exu_El = exu_El[,3] ; exv_El = exv_El[,3]

exwind_El = cbind(exlon,exlat,rep(NA,length(exlon)))
for(i in 1:nexlon)
{
    exwind_El[i,3] = sqrt((exu_El[i]^2)+(exv_El[i]^2))
}

############################################
# 2-3. Scalar Wind (Lap~WarmEl)
exu_WarmEl = get(paste("u",extraction,"_WarmEl",sep=""))
exv_WarmEl = get(paste("v",extraction,"_WarmEl",sep=""))

exu_WarmEl = exu_WarmEl[,3] ; exv_WarmEl = exv_WarmEl[,3]

exwind_WarmEl = cbind(exlon,exlat,rep(NA,length(exlon)))
for(i in 1:nexlon)
{
    exwind_WarmEl[i,3] = sqrt((exu_WarmEl[i]^2)+(exv_WarmEl[i]^2))
}

############################################
# 2-4. Scalar Wind (Lap~Cold)
exu_Cold = get(paste("u",extraction,"_Cold",sep=""))
exv_Cold = get(paste("v",extraction,"_Cold",sep=""))

exu_Cold = exu_Cold[,3] ; exv_Cold = exv_Cold[,3]

exwind_Cold = cbind(exlon,exlat,rep(NA,length(exlon)))
for(i in 1:nexlon)
{
    exwind_Cold[i,3] = sqrt((exu_Cold[i]^2)+(exv_Cold[i]^2))
}

############################################
# 2-5. Scalar Wind (Lap~Mean)
exu_Mean = get(paste("u",extraction,"_Mean",sep=""))
exv_Mean = get(paste("v",extraction,"_Mean",sep=""))

exu_Mean = exu_Mean[,3] ; exv_Mean = exv_Mean[,3]

exwind_Mean = cbind(exlon,exlat,rep(NA,length(exlon)))
for(i in 1:nexlon)
{
    exwind_Mean[i,3] = sqrt((exu_Mean[i]^2)+(exv_Mean[i]^2))
}

############################################
# 2-6. Scalar Wind (Lap~Warm)
exu_Warm = get(paste("u",extraction,"_Warm",sep=""))
exv_Warm = get(paste("v",extraction,"_Warm",sep=""))

exu_Warm = exu_Warm[,3] ; exv_Warm = exv_Warm[,3]

exwind_Warm = cbind(exlon,exlat,rep(NA,length(exlon)))
for(i in 1:nexlon)
{
    exwind_Warm[i,3] = sqrt((exu_Warm[i]^2)+(exv_Warm[i]^2))
}

############################################
# 2-7. Scalar Wind (Lap~ColdLa)
exu_ColdLa = get(paste("u",extraction,"_ColdLa",sep=""))
exv_ColdLa = get(paste("v",extraction,"_ColdLa",sep=""))

exu_ColdLa = exu_ColdLa[,3] ; exv_ColdLa = exv_ColdLa[,3]

exwind_ColdLa = cbind(exlon,exlat,rep(NA,length(exlon)))
for(i in 1:nexlon)
{
    exwind_ColdLa[i,3] = sqrt((exu_ColdLa[i]^2)+(exv_ColdLa[i]^2))
}

############################################
# 2-8. Scalar Wind (Lap~La)
exu_La = get(paste("u",extraction,"_La",sep=""))
exv_La = get(paste("v",extraction,"_La",sep=""))

exu_La = exu_La[,3] ; exv_La = exv_La[,3]

exwind_La = cbind(exlon,exlat,rep(NA,length(exlon)))
for(i in 1:nexlon)
{
    exwind_La[i,3] = sqrt((exu_La[i]^2)+(exv_La[i]^2))
}

############################################
# 2-9. Scalar Wind (Lap~WarmLa)
exu_WarmLa = get(paste("u",extraction,"_WarmLa",sep=""))
exv_WarmLa = get(paste("v",extraction,"_WarmLa",sep=""))

exu_WarmLa = exu_WarmLa[,3] ; exv_WarmLa = exv_WarmLa[,3]

exwind_WarmLa = cbind(exlon,exlat,rep(NA,length(exlon)))
for(i in 1:nexlon)
{
    exwind_WarmLa[i,3] = sqrt((exu_WarmLa[i]^2)+(exv_WarmLa[i]^2))
}

############################################
############################################
library(geosphere)
# Cal. partial x
p1=cbind(rep(0,73),seq(90,-90,-2.5)) ; p2=cbind(rep(2.5,73),seq(90,-90,-2.5))
dx=distVincentyEllipsoid(p1,p2,a=6378137, b=6356752.3142, f=1/298.257223563)

# Cal. partial y
imsi=c(90,seq(90,-90,-2.5),-90) ; n=length(imsi)
p1=cbind(rep(0,(n-2)),imsi[1:(n-2)]) ; p2= cbind(rep(0,(n-2)),imsi[3:n])
dy=(1/2)*distVincentyEllipsoid(p1,p2,a=6378137, b=6356752.3142, f=1/298.257223563)


# 3-1. Fin.value from Extraction wind (Lap~ColdEl)
imsiU = NULL ; imsiV = NULL
imsiU = cbind(exlon,exlat,exu_ColdEl) ; imsiV = cbind(exlon,exlat,exv_ColdEl)
imsiU = imsiU[c(order(imsiU[,1])),] ; imsiV = imsiV[c(order(imsiV[,1])),]

U = as.matrix(imsiU[,3]) ; V = as.matrix(imsiV[,3])

Lon = seq(0,357.5,2.5) ; ilen=length(Lon)
Lat = seq(90,-90,-2.5) ; jlen=length(Lat)

dUdx=NULL ; dVdx=NULL ; dUdy=NULL ; dVdy=NULL
dUdx_ann=NULL ; dVdx_ann=NULL ; dUdy_ann=NULL ; dVdy_ann=NULL
U_ann=matrix(U[,1],ncol=ilen) ; V_ann=matrix(V[,1],ncol=ilen)
for(j in 1:jlen) {d=diff(c(U_ann[j,ilen], U_ann[j,], U_ann[j,1])) ; dUdx_ann=c(dUdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(U_ann[1,i], U_ann[,i], U_ann[jlen,i])) ; dUdy_ann=c(dUdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
for(j in 1:jlen) {d=diff(c(V_ann[j,ilen], V_ann[j,], V_ann[j,1])) ; dVdx_ann=c(dVdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(V_ann[1,i], V_ann[,i], V_ann[jlen,i])) ; dVdy_ann=c(dVdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}

dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)
dUdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dUdx)
dVdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dVdx)
dUdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dUdy)
dVdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dVdy)

dUdy = dUdy[c(order(-dUdy[,2])),] ; dVdy = dVdy[c(order(-dVdy[,2])),]
dUdy[,3] = -1*dUdy[,3]
dVdy[,3] = -1*dVdy[,3]

divergence = dUdx[,3] + dVdy[,3]
#vorticity = (-1*dUdy[,3:(length(year)+2)]) - dVdx[,3:(length(year)+2)]
vorticity = dVdx[,3] - dUdy[,3]

exlon=rep(seq(0,357.5,2.5),73) ; exlat=rep(seq(90,-90,-2.5),each=144)
nexlon = length(exlon) ; nexlat = length(exlat)

divergence_ColdEl = cbind(exlon,exlat,divergence)
vorticity_ColdEl = cbind(exlon,exlat,vorticity)

rm(imsiU,imsiV,U,V,dUdx,dUdx_ann,dVdx,dVdx_ann,dUdy,dUdy_ann,dVdy,dVdy_ann)

############################################
# 3-2. Fin.value from Extraction wind (Lap~El)
imsiU = NULL ; imsiV = NULL
imsiU = cbind(exlon,exlat,exu_El) ; imsiV = cbind(exlon,exlat,exv_El)
imsiU = imsiU[c(order(imsiU[,1])),] ; imsiV = imsiV[c(order(imsiV[,1])),]

U = as.matrix(imsiU[,3]) ; V = as.matrix(imsiV[,3])

Lon = seq(0,357.5,2.5) ; ilen=length(Lon)
Lat = seq(90,-90,-2.5) ; jlen=length(Lat)

dUdx=NULL ; dVdx=NULL ; dUdy=NULL ; dVdy=NULL
dUdx_ann=NULL ; dVdx_ann=NULL ; dUdy_ann=NULL ; dVdy_ann=NULL
U_ann=matrix(U[,1],ncol=ilen) ; V_ann=matrix(V[,1],ncol=ilen)
for(j in 1:jlen) {d=diff(c(U_ann[j,ilen], U_ann[j,], U_ann[j,1])) ; dUdx_ann=c(dUdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(U_ann[1,i], U_ann[,i], U_ann[jlen,i])) ; dUdy_ann=c(dUdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
for(j in 1:jlen) {d=diff(c(V_ann[j,ilen], V_ann[j,], V_ann[j,1])) ; dVdx_ann=c(dVdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(V_ann[1,i], V_ann[,i], V_ann[jlen,i])) ; dVdy_ann=c(dVdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}

dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)
dUdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dUdx)
dVdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dVdx)
dUdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dUdy)
dVdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dVdy)

dUdy = dUdy[c(order(-dUdy[,2])),] ; dVdy = dVdy[c(order(-dVdy[,2])),]
dUdy[,3] = -1*dUdy[,3]
dVdy[,3] = -1*dVdy[,3]

divergence = dUdx[,3] + dVdy[,3]
#vorticity = (-1*dUdy[,3:(length(year)+2)]) - dVdx[,3:(length(year)+2)]
vorticity = dVdx[,3] - dUdy[,3]

exlon=rep(seq(0,357.5,2.5),73) ; exlat=rep(seq(90,-90,-2.5),each=144)
nexlon = length(exlon) ; nexlat = length(exlat)

divergence_El = cbind(exlon,exlat,divergence)
vorticity_El = cbind(exlon,exlat,vorticity)

rm(imsiU,imsiV,U,V,dUdx,dUdx_ann,dVdx,dVdx_ann,dUdy,dUdy_ann,dVdy,dVdy_ann)

############################################
# 3-3. Fin.value from Extraction wind (Lap~WarmEl)
imsiU = NULL ; imsiV = NULL
imsiU = cbind(exlon,exlat,exu_WarmEl) ; imsiV = cbind(exlon,exlat,exv_WarmEl)
imsiU = imsiU[c(order(imsiU[,1])),] ; imsiV = imsiV[c(order(imsiV[,1])),]

U = as.matrix(imsiU[,3]) ; V = as.matrix(imsiV[,3])

Lon = seq(0,357.5,2.5) ; ilen=length(Lon)
Lat = seq(90,-90,-2.5) ; jlen=length(Lat)

dUdx=NULL ; dVdx=NULL ; dUdy=NULL ; dVdy=NULL
dUdx_ann=NULL ; dVdx_ann=NULL ; dUdy_ann=NULL ; dVdy_ann=NULL
U_ann=matrix(U[,1],ncol=ilen) ; V_ann=matrix(V[,1],ncol=ilen)
for(j in 1:jlen) {d=diff(c(U_ann[j,ilen], U_ann[j,], U_ann[j,1])) ; dUdx_ann=c(dUdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(U_ann[1,i], U_ann[,i], U_ann[jlen,i])) ; dUdy_ann=c(dUdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
for(j in 1:jlen) {d=diff(c(V_ann[j,ilen], V_ann[j,], V_ann[j,1])) ; dVdx_ann=c(dVdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(V_ann[1,i], V_ann[,i], V_ann[jlen,i])) ; dVdy_ann=c(dVdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}

dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)
dUdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dUdx)
dVdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dVdx)
dUdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dUdy)
dVdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dVdy)

dUdy = dUdy[c(order(-dUdy[,2])),] ; dVdy = dVdy[c(order(-dVdy[,2])),]
dUdy[,3] = -1*dUdy[,3]
dVdy[,3] = -1*dVdy[,3]

divergence = dUdx[,3] + dVdy[,3]
#vorticity = (-1*dUdy[,3:(length(year)+2)]) - dVdx[,3:(length(year)+2)]
vorticity = dVdx[,3] - dUdy[,3]

exlon=rep(seq(0,357.5,2.5),73) ; exlat=rep(seq(90,-90,-2.5),each=144)
nexlon = length(exlon) ; nexlat = length(exlat)

divergence_WarmEl = cbind(exlon,exlat,divergence)
vorticity_WarmEl = cbind(exlon,exlat,vorticity)

fin_ColdEl = paste(fin,"_ColdEl",sep="")
fin_El = paste(fin,"_El",sep="")
fin_WarmEl = paste(fin,"_WarmEl",sep="")

rm(imsiU,imsiV,U,V,dUdx,dUdx_ann,dVdx,dVdx_ann,dUdy,dUdy_ann,dVdy,dVdy_ann)

############################################
# 3-4. Fin.value from Extraction wind (Lap~Cold)
imsiU = NULL ; imsiV = NULL
imsiU = cbind(exlon,exlat,exu_Cold) ; imsiV = cbind(exlon,exlat,exv_Cold)
imsiU = imsiU[c(order(imsiU[,1])),] ; imsiV = imsiV[c(order(imsiV[,1])),]

U = as.matrix(imsiU[,3]) ; V = as.matrix(imsiV[,3])

Lon = seq(0,357.5,2.5) ; ilen=length(Lon)
Lat = seq(90,-90,-2.5) ; jlen=length(Lat)

dUdx=NULL ; dVdx=NULL ; dUdy=NULL ; dVdy=NULL
dUdx_ann=NULL ; dVdx_ann=NULL ; dUdy_ann=NULL ; dVdy_ann=NULL
U_ann=matrix(U[,1],ncol=ilen) ; V_ann=matrix(V[,1],ncol=ilen)
for(j in 1:jlen) {d=diff(c(U_ann[j,ilen], U_ann[j,], U_ann[j,1])) ; dUdx_ann=c(dUdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(U_ann[1,i], U_ann[,i], U_ann[jlen,i])) ; dUdy_ann=c(dUdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
for(j in 1:jlen) {d=diff(c(V_ann[j,ilen], V_ann[j,], V_ann[j,1])) ; dVdx_ann=c(dVdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(V_ann[1,i], V_ann[,i], V_ann[jlen,i])) ; dVdy_ann=c(dVdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}

dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)
dUdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dUdx)
dVdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dVdx)
dUdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dUdy)
dVdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dVdy)

dUdy = dUdy[c(order(-dUdy[,2])),] ; dVdy = dVdy[c(order(-dVdy[,2])),]
dUdy[,3] = -1*dUdy[,3]
dVdy[,3] = -1*dVdy[,3]

divergence = dUdx[,3] + dVdy[,3]
#vorticity = (-1*dUdy[,3:(length(year)+2)]) - dVdx[,3:(length(year)+2)]
vorticity = dVdx[,3] - dUdy[,3]

exlon=rep(seq(0,357.5,2.5),73) ; exlat=rep(seq(90,-90,-2.5),each=144)
nexlon = length(exlon) ; nexlat = length(exlat)

divergence_Cold = cbind(exlon,exlat,divergence)
vorticity_Cold = cbind(exlon,exlat,vorticity)

rm(imsiU,imsiV,U,V,dUdx,dUdx_ann,dVdx,dVdx_ann,dUdy,dUdy_ann,dVdy,dVdy_ann)

############################################
# 3-5. Fin.value from Extraction wind (Lap~Mean)
imsiU = NULL ; imsiV = NULL
imsiU = cbind(exlon,exlat,exu_Mean) ; imsiV = cbind(exlon,exlat,exv_Mean)
imsiU = imsiU[c(order(imsiU[,1])),] ; imsiV = imsiV[c(order(imsiV[,1])),]

U = as.matrix(imsiU[,3]) ; V = as.matrix(imsiV[,3])

Lon = seq(0,357.5,2.5) ; ilen=length(Lon)
Lat = seq(90,-90,-2.5) ; jlen=length(Lat)

dUdx=NULL ; dVdx=NULL ; dUdy=NULL ; dVdy=NULL
dUdx_ann=NULL ; dVdx_ann=NULL ; dUdy_ann=NULL ; dVdy_ann=NULL
U_ann=matrix(U[,1],ncol=ilen) ; V_ann=matrix(V[,1],ncol=ilen)
for(j in 1:jlen) {d=diff(c(U_ann[j,ilen], U_ann[j,], U_ann[j,1])) ; dUdx_ann=c(dUdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(U_ann[1,i], U_ann[,i], U_ann[jlen,i])) ; dUdy_ann=c(dUdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
for(j in 1:jlen) {d=diff(c(V_ann[j,ilen], V_ann[j,], V_ann[j,1])) ; dVdx_ann=c(dVdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(V_ann[1,i], V_ann[,i], V_ann[jlen,i])) ; dVdy_ann=c(dVdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}

dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)
dUdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dUdx)
dVdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dVdx)
dUdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dUdy)
dVdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dVdy)

dUdy = dUdy[c(order(-dUdy[,2])),] ; dVdy = dVdy[c(order(-dVdy[,2])),]
dUdy[,3] = -1*dUdy[,3]
dVdy[,3] = -1*dVdy[,3]

divergence = dUdx[,3] + dVdy[,3]
#vorticity = (-1*dUdy[,3:(length(year)+2)]) - dVdx[,3:(length(year)+2)]
vorticity = dVdx[,3] - dUdy[,3]

exlon=rep(seq(0,357.5,2.5),73) ; exlat=rep(seq(90,-90,-2.5),each=144)
nexlon = length(exlon) ; nexlat = length(exlat)

divergence_Mean = cbind(exlon,exlat,divergence)
vorticity_Mean = cbind(exlon,exlat,vorticity)

rm(imsiU,imsiV,U,V,dUdx,dUdx_ann,dVdx,dVdx_ann,dUdy,dUdy_ann,dVdy,dVdy_ann)

############################################
# 3-6. Fin.value from Extraction wind (Lap~Warm)
imsiU = NULL ; imsiV = NULL
imsiU = cbind(exlon,exlat,exu_Warm) ; imsiV = cbind(exlon,exlat,exv_Warm)
imsiU = imsiU[c(order(imsiU[,1])),] ; imsiV = imsiV[c(order(imsiV[,1])),]

U = as.matrix(imsiU[,3]) ; V = as.matrix(imsiV[,3])

Lon = seq(0,357.5,2.5) ; ilen=length(Lon)
Lat = seq(90,-90,-2.5) ; jlen=length(Lat)

dUdx=NULL ; dVdx=NULL ; dUdy=NULL ; dVdy=NULL
dUdx_ann=NULL ; dVdx_ann=NULL ; dUdy_ann=NULL ; dVdy_ann=NULL
U_ann=matrix(U[,1],ncol=ilen) ; V_ann=matrix(V[,1],ncol=ilen)
for(j in 1:jlen) {d=diff(c(U_ann[j,ilen], U_ann[j,], U_ann[j,1])) ; dUdx_ann=c(dUdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(U_ann[1,i], U_ann[,i], U_ann[jlen,i])) ; dUdy_ann=c(dUdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
for(j in 1:jlen) {d=diff(c(V_ann[j,ilen], V_ann[j,], V_ann[j,1])) ; dVdx_ann=c(dVdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(V_ann[1,i], V_ann[,i], V_ann[jlen,i])) ; dVdy_ann=c(dVdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}

dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)
dUdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dUdx)
dVdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dVdx)
dUdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dUdy)
dVdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dVdy)

dUdy = dUdy[c(order(-dUdy[,2])),] ; dVdy = dVdy[c(order(-dVdy[,2])),]
dUdy[,3] = -1*dUdy[,3]
dVdy[,3] = -1*dVdy[,3]

divergence = dUdx[,3] + dVdy[,3]
#vorticity = (-1*dUdy[,3:(length(year)+2)]) - dVdx[,3:(length(year)+2)]
vorticity = dVdx[,3] - dUdy[,3]

exlon=rep(seq(0,357.5,2.5),73) ; exlat=rep(seq(90,-90,-2.5),each=144)
nexlon = length(exlon) ; nexlat = length(exlat)

divergence_Warm = cbind(exlon,exlat,divergence)
vorticity_Warm = cbind(exlon,exlat,vorticity)

fin_Cold = paste(fin,"_Cold",sep="")
fin_Mean = paste(fin,"_Mean",sep="")
fin_Warm = paste(fin,"_Warm",sep="")

rm(imsiU,imsiV,U,V,dUdx,dUdx_ann,dVdx,dVdx_ann,dUdy,dUdy_ann,dVdy,dVdy_ann)

############################################
# 3-7. Fin.value from Extraction wind (Lap~ColdLa)
imsiU = NULL ; imsiV = NULL
imsiU = cbind(exlon,exlat,exu_ColdLa) ; imsiV = cbind(exlon,exlat,exv_ColdLa)
imsiU = imsiU[c(order(imsiU[,1])),] ; imsiV = imsiV[c(order(imsiV[,1])),]

U = as.matrix(imsiU[,3]) ; V = as.matrix(imsiV[,3])

Lon = seq(0,357.5,2.5) ; ilen=length(Lon)
Lat = seq(90,-90,-2.5) ; jlen=length(Lat)

dUdx=NULL ; dVdx=NULL ; dUdy=NULL ; dVdy=NULL
dUdx_ann=NULL ; dVdx_ann=NULL ; dUdy_ann=NULL ; dVdy_ann=NULL
U_ann=matrix(U[,1],ncol=ilen) ; V_ann=matrix(V[,1],ncol=ilen)
for(j in 1:jlen) {d=diff(c(U_ann[j,ilen], U_ann[j,], U_ann[j,1])) ; dUdx_ann=c(dUdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(U_ann[1,i], U_ann[,i], U_ann[jlen,i])) ; dUdy_ann=c(dUdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
for(j in 1:jlen) {d=diff(c(V_ann[j,ilen], V_ann[j,], V_ann[j,1])) ; dVdx_ann=c(dVdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(V_ann[1,i], V_ann[,i], V_ann[jlen,i])) ; dVdy_ann=c(dVdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}

dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)
dUdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dUdx)
dVdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dVdx)
dUdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dUdy)
dVdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dVdy)

dUdy = dUdy[c(order(-dUdy[,2])),] ; dVdy = dVdy[c(order(-dVdy[,2])),]
dUdy[,3] = -1*dUdy[,3]
dVdy[,3] = -1*dVdy[,3]

divergence = dUdx[,3] + dVdy[,3]
#vorticity = (-1*dUdy[,3:(length(year)+2)]) - dVdx[,3:(length(year)+2)]
vorticity = dVdx[,3] - dUdy[,3]

exlon=rep(seq(0,357.5,2.5),73) ; exlat=rep(seq(90,-90,-2.5),each=144)
nexlon = length(exlon) ; nexlat = length(exlat)

divergence_ColdLa = cbind(exlon,exlat,divergence)
vorticity_ColdLa = cbind(exlon,exlat,vorticity)

rm(imsiU,imsiV,U,V,dUdx,dUdx_ann,dVdx,dVdx_ann,dUdy,dUdy_ann,dVdy,dVdy_ann)

############################################
# 3-8. Fin.value from Extraction wind (Lap~La)
imsiU = NULL ; imsiV = NULL
imsiU = cbind(exlon,exlat,exu_La) ; imsiV = cbind(exlon,exlat,exv_La)
imsiU = imsiU[c(order(imsiU[,1])),] ; imsiV = imsiV[c(order(imsiV[,1])),]

U = as.matrix(imsiU[,3]) ; V = as.matrix(imsiV[,3])

Lon = seq(0,357.5,2.5) ; ilen=length(Lon)
Lat = seq(90,-90,-2.5) ; jlen=length(Lat)

dUdx=NULL ; dVdx=NULL ; dUdy=NULL ; dVdy=NULL
dUdx_ann=NULL ; dVdx_ann=NULL ; dUdy_ann=NULL ; dVdy_ann=NULL
U_ann=matrix(U[,1],ncol=ilen) ; V_ann=matrix(V[,1],ncol=ilen)
for(j in 1:jlen) {d=diff(c(U_ann[j,ilen], U_ann[j,], U_ann[j,1])) ; dUdx_ann=c(dUdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(U_ann[1,i], U_ann[,i], U_ann[jlen,i])) ; dUdy_ann=c(dUdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
for(j in 1:jlen) {d=diff(c(V_ann[j,ilen], V_ann[j,], V_ann[j,1])) ; dVdx_ann=c(dVdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(V_ann[1,i], V_ann[,i], V_ann[jlen,i])) ; dVdy_ann=c(dVdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}

dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)
dUdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dUdx)
dVdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dVdx)
dUdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dUdy)
dVdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dVdy)

dUdy = dUdy[c(order(-dUdy[,2])),] ; dVdy = dVdy[c(order(-dVdy[,2])),]
dUdy[,3] = -1*dUdy[,3]
dVdy[,3] = -1*dVdy[,3]

divergence = dUdx[,3] + dVdy[,3]
#vorticity = (-1*dUdy[,3:(length(year)+2)]) - dVdx[,3:(length(year)+2)]
vorticity = dVdx[,3] - dUdy[,3]

exlon=rep(seq(0,357.5,2.5),73) ; exlat=rep(seq(90,-90,-2.5),each=144)
nexlon = length(exlon) ; nexlat = length(exlat)

divergence_La = cbind(exlon,exlat,divergence)
vorticity_La = cbind(exlon,exlat,vorticity)

rm(imsiU,imsiV,U,V,dUdx,dUdx_ann,dVdx,dVdx_ann,dUdy,dUdy_ann,dVdy,dVdy_ann)

############################################
# 3-9. Fin.value from Extraction wind (Lap~WarmLa)
imsiU = NULL ; imsiV = NULL
imsiU = cbind(exlon,exlat,exu_WarmLa) ; imsiV = cbind(exlon,exlat,exv_WarmLa)
imsiU = imsiU[c(order(imsiU[,1])),] ; imsiV = imsiV[c(order(imsiV[,1])),]

U = as.matrix(imsiU[,3]) ; V = as.matrix(imsiV[,3])

Lon = seq(0,357.5,2.5) ; ilen=length(Lon)
Lat = seq(90,-90,-2.5) ; jlen=length(Lat)

dUdx=NULL ; dVdx=NULL ; dUdy=NULL ; dVdy=NULL
dUdx_ann=NULL ; dVdx_ann=NULL ; dUdy_ann=NULL ; dVdy_ann=NULL
U_ann=matrix(U[,1],ncol=ilen) ; V_ann=matrix(V[,1],ncol=ilen)
for(j in 1:jlen) {d=diff(c(U_ann[j,ilen], U_ann[j,], U_ann[j,1])) ; dUdx_ann=c(dUdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(U_ann[1,i], U_ann[,i], U_ann[jlen,i])) ; dUdy_ann=c(dUdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}
for(j in 1:jlen) {d=diff(c(V_ann[j,ilen], V_ann[j,], V_ann[j,1])) ; dVdx_ann=c(dVdx_ann,(d[2:(ilen+1)]+d[1:ilen])/2*(1/dx[j]))}
for(i in 1:ilen) {d=diff(c(V_ann[1,i], V_ann[,i], V_ann[jlen,i])) ; dVdy_ann=c(dVdy_ann,(d[2:(jlen+1)]+d[1:jlen])/2*(1/dy))}

dUdx=cbind(dUdx,dUdx_ann) ; dVdx=cbind(dVdx,dVdx_ann) ; dUdy=cbind(dUdy,dUdy_ann) ; dVdy=cbind(dVdy,dVdy_ann)
dUdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dUdx)
dVdx = cbind(rep(seq(0,357.5,2.5),73),rep(seq(90,-90,-2.5),each=144),dVdx)
dUdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dUdy)
dVdy = cbind(rep(seq(0,357.5,2.5),each=73),rep(seq(90,-90,-2.5),144),dVdy)

dUdy = dUdy[c(order(-dUdy[,2])),] ; dVdy = dVdy[c(order(-dVdy[,2])),]
dUdy[,3] = -1*dUdy[,3]
dVdy[,3] = -1*dVdy[,3]

divergence = dUdx[,3] + dVdy[,3]
#vorticity = (-1*dUdy[,3:(length(year)+2)]) - dVdx[,3:(length(year)+2)]
vorticity = dVdx[,3] - dUdy[,3]

exlon=rep(seq(0,357.5,2.5),73) ; exlat=rep(seq(90,-90,-2.5),each=144)
nexlon = length(exlon) ; nexlat = length(exlat)

divergence_WarmLa = cbind(exlon,exlat,divergence)
vorticity_WarmLa = cbind(exlon,exlat,vorticity)

fin_ColdLa = paste(fin,"_ColdLa",sep="")
fin_La = paste(fin,"_La",sep="")
fin_WarmLa = paste(fin,"_WarmLa",sep="")

rm(imsiU,imsiV,U,V,dUdx,dUdx_ann,dVdx,dVdx_ann,dUdy,dUdy_ann,dVdy,dVdy_ann)

library(maps)
library(akima)

## Warning: package 'akima' was built under R version 4.1.1

library(fields)

## Warning: package 'fields' was built under R version 4.1.1

## Loading required package: spam

## Loading required package: dotCall64

## Spam version 2.6-0 (2020-12-14) 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: viridis

## Loading required package: viridisLite

## 
## Attaching package: 'viridis'

## The following object is masked from 'package:maps':
## 
##     unemp

## 
## Try help(fields) to get started.

library(RColorBrewer)

colist_negative = colorRampPalette(brewer.pal(9,'BuPu'))(10)
colist_positive = colorRampPalette(brewer.pal(9,'YlOrRd'))(10)
colist = c(rev(colist_negative),colist_positive)

if(laplace == "chi"){lapname = "Chi"}
if(laplace == "psi"){lapname = "Psi"}
if(extraction == "d"){windname = "Divergent wind"}
if(extraction == "r"){windname = "Rotational wind"}
if(fin == "divergence"){finname = "Divergence"}
if(fin == "vorticity"){finname = "Vorticity"}
if(year[length(year)] == 2012){yearname = "1991to2012"}
if(year[length(year)] == 2020){yearname = "1991to2020"}
############################################
# 1. Laplace Map (Chi or Psi)
lap_ColdEl_pred = lap_ColdEl_pred[145:10368,]
lap_El_pred = lap_El_pred[145:10368,]
lap_WarmEl_pred = lap_WarmEl_pred[145:10368,]

lap_Cold_pred = lap_Cold_pred[145:10368,]
lap_Mean = lap_Mean[145:10368,]
lap_Warm_pred = lap_Warm_pred[145:10368,]

lap_ColdLa_pred = lap_ColdLa_pred[145:10368,]
lap_La_pred = lap_La_pred[145:10368,]
lap_WarmLa_pred = lap_WarmLa_pred[145:10368,]

# 2. Extraction wind Map (wind_d or wind_r)
exwind_ColdEl = exwind_ColdEl[145:10368,]
exwind_El = exwind_El[145:10368,]
exwind_WarmEl = exwind_WarmEl[145:10368,]

exwind_Cold = exwind_Cold[145:10368,]
exwind_Mean = exwind_Mean[145:10368,]
exwind_Warm = exwind_Warm[145:10368,]

exwind_ColdLa = exwind_ColdLa[145:10368,]
exwind_La = exwind_La[145:10368,]
exwind_WarmLa = exwind_WarmLa[145:10368,]

# 3. Fin val. Map (Chi or Psi)
finvar_ColdEl = get(fin_ColdEl) ; finvar_ColdEl = finvar_ColdEl[145:10368,]
finvar_El = get(fin_El) ; finvar_El = finvar_El[145:10368,]
finvar_WarmEl = get(fin_WarmEl) ; finvar_WarmEl = finvar_WarmEl[145:10368,]

finvar_Cold = get(fin_Cold) ; finvar_Cold = finvar_Cold[145:10368,]
finvar_Mean = get(fin_Mean) ; finvar_Mean = finvar_Mean[145:10368,]
finvar_Warm = get(fin_Warm) ; finvar_Warm = finvar_Warm[145:10368,]

finvar_ColdLa = get(fin_ColdLa) ; finvar_ColdLa = finvar_ColdLa[145:10368,]
finvar_La = get(fin_La) ; finvar_La = finvar_La[145:10368,]
finvar_WarmLa = get(fin_WarmLa) ; finvar_WarmLa = finvar_WarmLa[145:10368,]
############################################
# Chi or Psi : lap_(mode name)_pred
# Wind (D or R) : exwind_(mode name)
# Fin val. (delta or zeta) : finvar_(mode name)

# Fig.3 a
A1 = finvar_Mean[,1]    # Longitude
A2 = finvar_Mean[,2]    # Latitude
A3 = finvar_Mean[,3]    # Vorticity
A4 = lap_Mean[,3]   # Streamfunction

par(mar=c(6,7,4,5))

plot(-999,-999,xlim=c(100,185),ylim=c(0,50),xlab='',ylab='',axes=F,xaxs='i',yaxs='i')
image(interp(A1,A2,A3,xo=seq(min(A1),max(A1),length=144*15),yo=seq(min(A2),max(A2),length=73*15),duplicate='mean'),axes=F,xlab='',ylab='',breaks=c(min(A3),seq(-2e-05,2e-05,length=19),max(A3)),col=colist,main='',xlim=c(100,180),ylim=c(0,50),add=T)

m=map('world',col=colors()[213],interior=F,add=T,lty=0,xlim=c(100,180),ylim=c(0,50))
map(m,boundary=T,interior=F,add=T,col=colors()[213],lwd=2,xlim=c(100,180),ylim=c(0,50))

abline(v=seq(0,360,20),lty=1,col=grey(0.8))
abline(h=seq(-90,90,10),lty=1,col=grey(0.8))

contour(interp(A1,A2,A4,xo=seq(min(A1),max(A1),length=144*15),yo=seq(min(A2),max(A2),length=73*15),duplicate='mean'),levels=c(-2e06,0),drawlabels=T,labels=c(-2,0),lwd=1,add=T,col='Dark Blue',labcex=0.9)
#grey(0.6)
contour(interp(A1,A2,A4,xo=seq(min(A1),max(A1),length=144*15),yo=seq(min(A2),max(A2),length=73*15),duplicate='mean'),levels=0,drawlabels=T,labels=0,lwd=2,add=T,col='Dark Blue',labcex=0.9,method='flattest')
contour(interp(A1,A2,A4,xo=seq(min(A1),max(A1),length=144*15),yo=seq(min(A2),max(A2),length=73*15),duplicate='mean'),levels=c(2e06,4e06,6e06,8e06),drawlabels=T,labels=c(2,4,6,8),lwd=1,add=T,col='Dark Blue',labcex=0.9,method='edge')

axis(1,at=seq(100,180,20),labels=expression(100~degree~E,120~degree~E,140~degree~E,160~degree~E,180~degree),cex.axis=1.2)
axis(2,at=seq(-90,90,10),labels=expression(90~degree~S,80~degree~S,70~degree~S,60~degree~S,50~degree~S,40~degree~S,30~degree~S,20~degree~S,10~degree~S,EQ.,10~degree~N,20~degree~N,30~degree~N,40~degree~N,50~degree~N,60~degree~N,70~degree~N,80~degree~N,90~degree~N),cex.axis=1.2,las=1)
axis(3,at=seq(100,180,20),labels=rep('',length(seq(100,180,20))))

mtext('Longitude',at=c(140,-10),side=1,line=3.5,cex=1.5)
mtext('Latitude',side=2,line=4.8,cex=1.5)

text(85,60,'a',cex=2.8,font=2,xpd=T)

rect(180,0,185,50,col='white',border='white',xpd=T)

imsicol = c("#8C82BC", "#A9C3DE", "#E2EDF5", "#FFFFCC", "#FEDD7F", "#FD9D43")
image.plot(zlim=c(-12,12),legend.only=T,nlevel=12,col=imsicol,breaks=seq(-12,12,length=7),horizontal=F,legend.args=list(text=expression(paste(zeta,' (x',10^-6,')',sep='')),side=3,line=0.4,cex=0.8),legend.cex=0.6,legend.width=1.2,legend.shrink=1,lab.breaks=c('',-8,-4,0,4,8,''))

rect(100,0,180,50,lty=1,lwd=1,,xpd=T,lend='square',ljoin='mitre')
rect(100,20,180,30,lty=1,lwd=4,border=grey(0.3),xpd=T,lend='square',ljoin='mitre')

# Fig.3 b

# 1. Lat. mean
# extraction for Lat.mean
lap_CE_pred = lap_ColdEl_pred[3313:4032,]
lap_E_pred = lap_El_pred[3313:4032,]
lap_WE_pred = lap_WarmEl_pred[3313:4032,]

lap_C_pred = lap_Cold_pred[3313:4032,]
lap_N_pred = lap_Mean[3313:4032,]   # it's not predict value (setting for loop)
lap_W_pred = lap_Warm_pred[3313:4032,]

lap_CL_pred = lap_ColdLa_pred[3313:4032,]
lap_L_pred = lap_La_pred[3313:4032,]
lap_WL_pred = lap_WarmLa_pred[3313:4032,]

EVmode = c('CE','E','WE','C','N','W','CL','L','WL')
lapval = paste('lap_',EVmode,'_pred',sep="")
EVlap = array(rep(NA,144*9),dim=c(144,9))
for(i in 1:length(EVmode))
{
    LM = get(lapval[i]) # Lat. Mean
    lapmat = matrix(get(lapval[i])[,3],ncol=5)
    for(j in 1:144)
    {
        EVlap[j,i]=mean(lapmat[j,])
    }
}
colnames(EVlap) = EVmode ; lon = seq(0,357.5,2.5) ; EVlap = cbind(lon,EVlap)
EVlap = EVlap[41:73,]   # range : -3.5e06 ~ 1.4e07

# Shading
Warmcol = adjustcolor('Dark Orange',alpha=0.5)
Coldcol = adjustcolor('#0070c0',alpha=0.5)
BaseMode = c(3,6,9) # Number of Base Phase (El Niño, Mean, La Niña)

############################################
par(mar=c(6,7,4,5))

plot(-2e07,-2e07,xlim=c(100,185),ylim=c(-4e06,15e06),xlab='',ylab='',axes=F,xaxs='i',yaxs='i')
xline=c(100:180,rep(0,19)) ; yline=rep(0,100)
abline(v=seq(100,180,20),lty=1,col=grey(0.8))
lines(xline,yline,lwd=2,col='Dark Blue',lend='square',ljoin='mitre')

for(i in BaseMode)
{
polygon(c(EVlap[,1],rev(EVlap[,1])),c(EVlap[,i+1],rev(EVlap[,i])),border=NA,col=Warmcol)
polygon(c(EVlap[,1],rev(EVlap[,1])),c(EVlap[,i-1],rev(EVlap[,i])),border=NA,col=Coldcol)
    if(i == 6)
    {lines(EVlap[,1],EVlap[,i],col='Black',lwd=3,lend='square')}
    else
    {lines(EVlap[,1],EVlap[,i],col='Black',lwd=1,lend='square')}
    
    lines(EVlap[,1],EVlap[,i-1],col='#0070c0',lwd=1,lend='square')
    lines(EVlap[,1],EVlap[,i+1],col='Dark Orange',lwd=1,lend='square')
}

rect(180,4e06,185,15e06,col='white',border='white',xpd=T)

axis(1,at=seq(100,180,20),labels=expression(100~degree~E,120~degree~E,140~degree~E,160~degree~E,180~degree),cex.axis=1.2)
axis(2,at=seq(-18e06,18e06,3e06),labels=seq(-18,18,3),cex.axis=1.2,las=1)
axis(3,at=seq(100,180,20),labels=rep('',length(seq(100,180,20))))
mtext('Longitude',at=c(140,-6e06),side=1,line=3.5,cex=1.5)
mtext(expression(paste('Streamfunction (',10^6,' ',m^2,s^-1,')')),side=2,line=3.5,cex=1.5)

############################################
xloc=182 ; cirsize = 1.7 ; texsize = 0.5
# 2.1 sigma CE mode
points(xloc,EVlap[33,2]-5e05,pch=16,cex=cirsize,col='#0070c0',xpd=T)
points(xloc,EVlap[33,2]-5e05,pch='1',cex=texsize,col='white',xpd=T)

# 1.5 sigma E mode
points(xloc,EVlap[33,3]-2e05,pch=16,cex=cirsize,col='Black',xpd=T)
points(xloc,EVlap[33,3]-2e05,pch='2',cex=texsize,col='white',xpd=T)

# 2.1 sigma WE mode
points(xloc,EVlap[33,4]+2e05,pch=16,cex=cirsize,col='Dark Orange',xpd=T)
points(xloc,EVlap[33,4]+2e05,pch='3',cex=texsize,col='white',xpd=T)

############################################
# 1.5 sigma C mode
points(xloc,EVlap[33,5]-6e05,pch=16,cex=cirsize,col='#0070c0',xpd=T)
points(xloc,EVlap[33,5]-6e05,pch='4',cex=texsize,col='white',xpd=T)

# N phase (This is not mode)
points(xloc,EVlap[33,6],pch=16,cex=cirsize,col='Black',xpd=T)
points(xloc,EVlap[33,6],pch='N',cex=texsize,col='white',xpd=T)

# 1.5 sigma W mode
points(xloc,EVlap[33,7]+6e05,pch=16,cex=cirsize,col='Dark Orange',xpd=T)
points(xloc,EVlap[33,7]+6e05,pch='5',cex=texsize,col='white',xpd=T)

############################################
# 2.1 sigma CL mode
points(xloc,EVlap[33,8]-2.5e05,pch=16,cex=cirsize,col='#0070c0',xpd=T)
points(xloc,EVlap[33,8]-2.5e05,pch='6',cex=texsize,col='white',xpd=T)

# 1.5 sigma L mode
points(xloc,EVlap[33,9]+1.5e05,pch=16,cex=cirsize,col='Black',xpd=T)
points(xloc,EVlap[33,9]+1.5e05,pch='7',cex=texsize,col='white',xpd=T)

# 2.1 sigma WL mode
points(xloc,EVlap[33,10]+4.5e05,pch=16,cex=cirsize,col='Dark Orange',xpd=T)
points(xloc,EVlap[33,10]+4.5e05,pch='8',cex=texsize,col='white',xpd=T)

############################################
# Legend
cirsize = 2.2 ; texsize = 1
legpo = 13.8e06
points(104, legpo,pch=16,cex=cirsize,col='#0070c0',xpd=T)
points(104, legpo,pch='1',cex=texsize,col='white',xpd=T)
text(108, legpo,expression(paste('Colder El Ni',tilde(n),'o',sep='')),cex=0.7,xpd=T,adj=0)

points(104, legpo-1.2e06,pch=16,cex=cirsize,col='Black',xpd=T)
points(104, legpo-1.2e06,pch='2',cex=texsize,col='white',xpd=T)
text(108, legpo-1.2e06,expression(paste('El Ni',tilde(n),'o',sep='')),cex=0.7,xpd=T,adj=0)

points(104, legpo-2.4e06,pch=16,cex=cirsize,col='Dark Orange',xpd=T)
points(104, legpo-2.4e06,pch='3',cex=texsize,col='white',xpd=T)
text(108, legpo-2.4e06,expression(paste('Warmer El Ni',tilde(n),'o',sep='')),cex=0.7,xpd=T,adj=0)

points(104, legpo-3.6e06,pch=16,cex=cirsize,col='#0070c0',xpd=T)
points(104, legpo-3.6e06,pch='4',cex=texsize,col='white',xpd=T)
text(108, legpo-3.6e06,'Warmer',cex=0.7,xpd=T,adj=0)

points(104, legpo-4.8e06,pch=16,cex=cirsize,col='Black',xpd=T)
points(104, legpo-4.8e06,pch='N',cex=texsize,col='white',xpd=T)
text(108, legpo-4.8e06,'Normal',cex=0.7,xpd=T,adj=0)

points(104, legpo-6e06,pch=16,cex=cirsize,col='Dark Orange',xpd=T)
points(104, legpo-6e06,pch='5',cex=texsize,col='white',xpd=T)
text(108, legpo-6e06,'Colder',cex=0.7,xpd=T,adj=0)

points(104, legpo-7.2e06,pch=16,cex=cirsize,col='#0070c0',xpd=T)
points(104, legpo-7.2e06,pch='6',cex=texsize,col='white',xpd=T)
text(108, legpo-7.2e06,expression(paste('Colder La Ni',tilde(n),'a',sep='')),cex=0.7,xpd=T,adj=0)

points(104, legpo-8.4e06,pch=16,cex=cirsize,col='Black',xpd=T)
points(104, legpo-8.4e06,pch='7',cex=texsize,col='white',xpd=T)
text(108, legpo-8.4e06,expression(paste('La Ni',tilde(n),'a',sep='')),cex=0.7,xpd=T,adj=0)

points(104, legpo-9.6e06,pch=16,cex=cirsize,col='Dark Orange',xpd=T)
points(104, legpo-9.6e06,pch='8',cex=texsize,col='white',xpd=T)
text(108, legpo-9.6e06,expression(paste('Warmer La Ni',tilde(n),'a',sep='')),cex=0.7,xpd=T,adj=0)

text(85,18.5e06,'b',cex=2.8,font=2,xpd=T)
rect(100,-4e06,180,15e06,lty=1,lwd=1,,xpd=T,lend='square',ljoin='mitre')

Fig. 3. Geographical patterns of streamfunction. (a) Mean fields of streamfunction (\(\psi\)) and vorticity (\(\zeta\)), and (b) the sequential patterns of \(\psi\) corresponding to the eight modes of variability. All values are originally calculated in a global domain, and zoom in for the western North Pacific. In (a), local vorticities are shaded following the color legend on the right side, and streamlines are contoured at 10\(^6\) m\(^2\) s\(^{-1}\) intervals. A thick dark blue lines indicates a streamline having zero values. The lines in (b) show each longitudinal distribution of \(\psi\) averaged over the latitudes between 20\(^\circ\)–30\(^\circ\)N as denoted by a black rectange in (a). The thick black line represents the zonal distribution of the mean field as a reference, and the two thin black lines show the ENSO anomalies. Orange and blue shades represent the warm and cold anomalies, respectively. It is found that the patterns are sequentially continuous showing a geographical spectrum, which implies the narrowing region for the tropical convections as the chain number moves from 1 to 8.

Figure 4. Longitudinal distribution of vorticity

# 1. Lat. mean
# extraction for Lat.mean
finvar_CE = finvar_ColdEl[3889:4896,]
finvar_E = finvar_El[3889:4896,]
finvar_WE = finvar_WarmEl[3889:4896,]

finvar_C = finvar_Cold[3889:4896,]
finvar_N = finvar_Mean[3889:4896,]
finvar_W = finvar_Warm[3889:4896,]

finvar_CL = finvar_ColdLa[3889:4896,]
finvar_L = finvar_La[3889:4896,]
finvar_WL = finvar_WarmLa[3889:4896,]

EVmode = c('CE','E','WE','C','N','W','CL','L','WL')
finval = paste('finvar_',EVmode,sep="")
EVfin = array(rep(NA,144*9),dim=c(144,9))
for(i in 1:length(EVmode))
{
    LM = get(finval[i]) # Lat. Mean
    finmat = matrix(get(finval[i])[,3],ncol=7)
    for(j in 1:144)
    {
        EVfin[j,i]=mean(finmat[j,])
    }
}
colnames(EVfin) = EVmode ; lon = seq(0,357.5,2.5) ; EVfin = cbind(lon,EVfin)
EVfin = EVfin[41:73,]   # range : -3.5e06 ~ 1.4e07

# Shading
Warmcol = adjustcolor('Dark Orange',alpha=0.4)
Coldcol = adjustcolor('Dark Blue',alpha=0.4)
Meancol = adjustcolor('Black',alpha=0.4)
BaseMode = c(3,6,9) # Number of Base Phase (El Niño, Mean, La Niña)

# Fig.4 a
#######################################################################################
# 1. El Niño phase
par(mfrow=c(1,1),mar=c(6,7,4,5))

plot(-2e07,-2e07,xlim=c(100,185),ylim=c(-3e-06,6e-06),xlab='',ylab='',axes=F,xaxs='i',yaxs='i')
xline=c(100:180,rep(0,19)) ; yline=rep(0,100)
lines(xline,yline,lwd=1,col=1,lend='square',ljoin='mitre')
abline(v=seq(100,180,20),lty=1,col=grey(0.8))
for(i in BaseMode)
{
  if(i == 6)
  {lines(EVfin[,1],EVfin[,i],col=grey(0.9),lwd=1.5,lend='square')}
  else
  {lines(EVfin[,1],EVfin[,i],col=grey(0.9),lwd=1,lend='square')}
  
  lines(EVfin[,1],EVfin[,i-1],col=grey(0.9),lwd=1,lend='square')
  lines(EVfin[,1],EVfin[,i+1],col=grey(0.9),lwd=1,lend='square')
}
lines(EVfin[,1],EVfin[,2],col='Dark Blue',lwd=1,lend='square')
lines(EVfin[,1],EVfin[,3],col='Black',lwd=1,lend='square')
lines(EVfin[,1],EVfin[,4],col='Dark Orange',lwd=1,lend='square')

rect(180,3e-06,185,6e-06,col='white',border='white',xpd=T)

axis(1,at=seq(100,180,20),labels=expression(100~degree~E,120~degree~E,140~degree~E,160~degree~E,180~degree),cex.axis=1.2)
axis(2,at=seq(-18e-06,18e-06,3e-06),labels=seq(-18,18,3),cex.axis=1.2,las=1)
axis(3,at=seq(100,180,20),labels=rep('',length(seq(100,180,20))))
mtext('Longitude',at=c(140,-6e06),side=1,line=3,cex=1)
mtext(expression(paste('Vorticity (',10^-6,' ',s^-1,')')),side=2,line=3.5,cex=1)

text(85,7.5e-06,'a',cex=2.8,font=2,xpd=T)
text(146.5,4.3e-06,'CE',cex=1.2,xpd=T,adj=0,col='Dark Blue')
text(156.5,4.3e-06,'E',cex=1.2,xpd=T,adj=0,col='Black')
text(162.5,4.3e-06,'WE',cex=1.2,xpd=T,adj=0,col='Dark Orange')
rect(100,-3e-06,180,6e-06,lty=1,lwd=1,,xpd=T,lend='square',ljoin='mitre')

# Fig.4 b
#######################################################################################
# 2. Normal phase
par(mfrow=c(1,1),mar=c(6,7,4,5))

plot(-2e07,-2e07,xlim=c(100,185),ylim=c(-3e-06,6e-06),xlab='',ylab='',axes=F,xaxs='i',yaxs='i')
xline=c(100:180,rep(0,19)) ; yline=rep(0,100)
lines(xline,yline,lwd=1,col=1,lend='square',ljoin='mitre')
abline(v=seq(100,180,20),lty=1,col=grey(0.8))
for(i in BaseMode)
{
  if(i == 6)
  {lines(EVfin[,1],EVfin[,i],col=grey(0.9),lwd=1.5,lend='square')}
  else
  {lines(EVfin[,1],EVfin[,i],col=grey(0.9),lwd=1,lend='square')}
  
  lines(EVfin[,1],EVfin[,i-1],col=grey(0.9),lwd=1,lend='square')
  lines(EVfin[,1],EVfin[,i+1],col=grey(0.9),lwd=1,lend='square')
}
lines(EVfin[,1],EVfin[,5],col='Dark Blue',lwd=1,lend='square')
lines(EVfin[,1],EVfin[,6],col='Black',lwd=1.5,lend='square')
lines(EVfin[,1],EVfin[,7],col='Dark Orange',lwd=1,lend='square')

rect(180,3e-06,185,6e-06,col='white',border='white',xpd=T)

axis(1,at=seq(100,180,20),labels=expression(100~degree~E,120~degree~E,140~degree~E,160~degree~E,180~degree),cex.axis=1.2)
axis(2,at=seq(-18e-06,18e-06,3e-06),labels=seq(-18,18,3),cex.axis=1.2,las=1)
axis(3,at=seq(100,180,20),labels=rep('',length(seq(100,180,20))))
mtext('Longitude',at=c(140,-6e06),side=1,line=3,cex=1)
mtext(expression(paste('Vorticity (',10^-6,' ',s^-1,')')),side=2,line=3.5,cex=1)

text(85,7.5e-06,'b',cex=2.8,font=2,xpd=T)
text(146.5,2.2e-06,'C',cex=1.2,xpd=T,adj=0,col='Dark Blue')
text(154.5,2.2e-06,'N',cex=1.2,xpd=T,adj=0,col='Black')
text(162.5,2.2e-06,'W',cex=1.2,xpd=T,adj=0,col='Dark Orange')
rect(100,-3e-06,180,6e-06,lty=1,lwd=1,,xpd=T,lend='square',ljoin='mitre')

# Fig.4 c
#######################################################################################
# 3. La Niña phase
par(mfrow=c(1,1),mar=c(6,7,4,5))

plot(-2e07,-2e07,xlim=c(100,185),ylim=c(-3e-06,6e-06),xlab='',ylab='',axes=F,xaxs='i',yaxs='i')
xline=c(100:180,rep(0,19)) ; yline=rep(0,100)
lines(xline,yline,lwd=1,col=1,lend='square',ljoin='mitre')
abline(v=seq(100,180,20),lty=1,col=grey(0.8))
for(i in BaseMode)
{
  if(i == 6)
  {lines(EVfin[,1],EVfin[,i],col=grey(0.9),lwd=1.5,lend='square')}
  else
  {lines(EVfin[,1],EVfin[,i],col=grey(0.9),lwd=1,lend='square')}
  
  lines(EVfin[,1],EVfin[,i-1],col=grey(0.9),lwd=1,lend='square')
  lines(EVfin[,1],EVfin[,i+1],col=grey(0.9),lwd=1,lend='square')
}
lines(EVfin[,1],EVfin[,8],col='Dark Blue',lwd=1,lend='square')
lines(EVfin[,1],EVfin[,9],col='Black',lwd=1,lend='square')
lines(EVfin[,1],EVfin[,10],col='Dark Orange',lwd=1,lend='square')

rect(180,3e-06,185,6e-06,col='white',border='white',xpd=T)

axis(1,at=seq(100,180,20),labels=expression(100~degree~E,120~degree~E,140~degree~E,160~degree~E,180~degree),cex.axis=1.2)
axis(2,at=seq(-18e-06,18e-06,3e-06),labels=seq(-18,18,3),cex.axis=1.2,las=1)
axis(3,at=seq(100,180,20),labels=rep('',length(seq(100,180,20))))
mtext('Longitude',at=c(140,-6e06),side=1,line=3,cex=1)
mtext(expression(paste('Vorticity (',10^-6,' ',s^-1,')')),side=2,line=3.5,cex=1)

text(85,7.5e-06,'c',cex=2.8,font=2,xpd=T)
text(146.5,0.5e-06,'CL',cex=1.2,xpd=T,adj=0,col='Dark Blue')
text(156.5,0.5e-06,'L',cex=1.2,xpd=T,adj=0,col='Black')
text(162.5,0.5e-06,'WL',cex=1.2,xpd=T,adj=0,col='Dark Orange')
rect(100,-3e-06,180,6e-06,lty=1,lwd=1,,xpd=T,lend='square',ljoin='mitre')

Fig. 4. Longitudinal distribution of vorticity. The local vorticities between 5\(^\circ\)–20\(^\circ\)N are averaged along 100\(^\circ\)E–180\(^\circ\) for (a) “El Niño mode (E) and its anomaly modes (CE, WE)”, (b) “Colder mode and Warmer mode (C, W)” around Normal (N), and (c) “La Niña mode and its anomaly modes (CL, WL)”. The black lines from top to bottom represent E, N, and L, and the blue and red lines in each panel show the colder and ther warmer anomalies, respectively.

Figure 5. Chronology of the modes of environmental variability

rm(list=ls())

yearN=1:36  #1985~2020 (36 years)
GMSST=read.table('GMSST_6.txt',T)[yearN,1]
SOI=read.table('SOI_6.txt',T)[yearN,1] ; NSOI=-SOI
PC1=(scale(GMSST)+scale(NSOI))/sqrt(2)
PC2=(scale(GMSST)-scale(NSOI))/sqrt(2)

#Assigning the climate condition to each year
CE   = -PC2                 # Cold El Niño
E     = scale(NSOI)         # El Niño
WE  = PC1                   # Warm El Niño
C    = scale(-GMSST)      # Cold Normal
W   = scale(GMSST)       # Warm Normal
CL   = -PC1                 # Cold La Niña
L     = scale(-NSOI)        # La Niña
WL  =  PC2                  # Warm La Niña

ClimCond=cbind(WL,L,CL,W,C,WE,E,CE)
patternN=NULL ; level=NULL 
for(i in 1:36) {imsi=ClimCond[i,] ; max=max(imsi) ; level=c(level,max) ; patternN=c(patternN,mean(which.max(imsi)))}

par(mar=c(6,8,4,6))
plot(-999,-999,xlim=c(0.5,8.5),ylim=c(2020,1985), axes=FALSE, xlab='', ylab='') 
abline(v=1:8,lty=2,col='grey')
lines(patternN,1985:2020,lwd=1,lty=1,col='black')
points(patternN,1985:2020,pch=16,cex=2*level, col='dark orange')
points(patternN,1985:2020,pch=21,cex=2*level, col=1)
legend( 'topleft', c(expression(paste('1.0', sigma)), expression(paste( '2.0',sigma)))
         ,cex=1.1, pt.cex=c(2,4), pch=21, col=1, pt.bg='dark orange', x.intersp=1.2, y.intersp=1.3)

axis(1, at = 1:8, labels= c('WL','L','CL','W','C','WE','E','CE'), cex.axis=1.1)
axis(2, at = seq(1950,2030,5), labels= seq(1950,2030,5),las=1, cex.axis=1.1)
axis(3, at = 1:8, labels= c('8','7','6','5','4','3','2','1'), cex.axis=1.1)
axis(4, at = seq(1950,2030,5), labels= seq(1950,2030,5),las=1, cex.axis=1.1)
mtext('Mode',side=1,line=3, cex=1.5)
mtext('Year',side=2,line=4, cex=1.5)

box()

Fig. 5. Chronology of the modes of environmental variability. For each year, the closest variability direction among the eight modes is shown with its magnitude. The smaller circle implies that the variability is less anomalous and closer to N. El Niño and La Niña trends repeat over time and tilt toward each warmer mode such as WE and WL. Considering the westernmost expansion of the atmospheric pattern in WL, a linear perspective suggests that every warmer La Niña environment could be one which people have scarcely experienced before.

END