一、前言
原代码在NCl官网MJO板块可以找到：NCL: Madden Julian Oscillation Climate Variability
分别是mjoclivar2，mjoclivar14，mjoclivar16.
所需要的数据为olr，u850，v850和u200.我这里用的是ERA5的观测数据
先上图片
二、mjoclivar2部分
分别对olr，u850和v850求取异常场，需要跑三次以下代码，把文件名和变量名改一改就行，得到输出的三个异常场的nc文件。
originmjoclivar_2.ncl
;-------------------------------------------------------------
; These files are loaded by default in NCL V6.2.0 and newer
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;******************************************************
begin
;******************************************************
; User specifications
;******************************************************
NC = True ; create netCDF?
PLOT = True ; sample plots?
ymdStrt = 20000101 ; start yyyymmdd
ymdLast = 20191231 ; last
yrStrt = ymdStrt/10000
yrLast = ymdLast/10000
nhar = 4 ; number of fourier comp
var = "olr" ; name of file
vName = "OLR" ; name for plots
;diri = "/project/cas/shea/CDC/" ; input directory
;diri = "/Users/shea/Data/AMWG/" ; input directory
diri = "/mnt/c/project/MJO/data/" ; new input directory
fili = "noaa_daily_olr_1974_2019.nc" ; input file
if (NC) then
diro= "/mnt/c/project/MJO/data/" ; output dir
;filo= var+".day.anomalies.nc" ; output file
filo= vName+".day.anomalies."+yrStrt+"-"+yrLast+".nc" ; output file
end if
if (PLOT) then
wksType = "png" ; send graphics to PNG file
wksName = "mjoclivar" ; "mjo."+yrStrt+"_"+yrLast
end if
;***********************************************************
; Read user specified time and create required yyyyddd
;***********************************************************
f = addfile (diri+fili , "r")
time = f->time ; all times on file
ymd = cd_calendar(time, -2) ; yyyymmdd
iStrt = ind(ymd.eq.ymdStrt) ; index start
iLast = ind(ymd.eq.ymdLast) ; index last
delete(time)
delete(ymd)
;***********************************************************
; Read user specified time and create required yyyyddd
;***********************************************************
time = f->time(iStrt:iLast) ; time:units = "hours since"
TIME = cd_calendar(time, 0) ; type float
year = floattointeger( TIME(:,0) )
month = floattointeger( TIME(:,1) )
day = floattointeger( TIME(:,2) )
ddd = day_of_year(year, month, day)
yyyyddd = year*1000 + ddd ; needed for input
;***********************************************************
; Read data: short2flt
;***********************************************************
x = short2flt( f->$var$(iStrt:iLast,:,:) ) ; convert to float
printVarSummary( x )
;***********************************************************
; Compute daily climatology: raw and then 'smoothed'
;***********************************************************
xClmDay = clmDayTLL(x, yyyyddd) ; daily climatology at each grid point
printVarSummary(xClmDay)
;***********************************************************
; Compute smoothed daily climatology using 'nhar' harmonics
;***********************************************************
xClmDay_sm = smthClmDayTLL(xClmDay, nhar)
printVarSummary(xClmDay_sm)
;***********************************************************
; Compute daily anomalies using raw and smoothed climatologies
;***********************************************************
xAnom = calcDayAnomTLL (x, yyyyddd, xClmDay)
printVarSummary(xAnom)
printMinMax(xAnom, True)
xAnom_sm = calcDayAnomTLL (x, yyyyddd, xClmDay_sm)
xAnom_sm@long_name = "Anomalies from Smooth Daily Climatology"
printVarSummary(xAnom_sm)
printMinMax(xAnom_sm, True)
delete( x ) ; no longer needed
;***********************************************************
; Create netCDF: convenience use 'simple' method
;***********************************************************
dimx = dimsizes(xAnom)
ntim = dimx(0)
nlat = dimx(1)
mlon = dimx(2)
if (NC) then
lat = f->lat
lon = f->lon
system("/bin/rm -f "+diro+filo) ; rm any pre-exist file, if any
fnc = addfile (diro+filo, "c")
filAtt = 0
filAtt@title = vName+": Daily Anomalies: "+yrStrt+"-"+yrLast
filAtt@source_file = fili
filAtt@creation_date = systemfunc("date")
fileattdef( fnc, filAtt ) ; copy file attributes
setfileoption(fnc,"DefineMode",True)
varNC = vName+"_anom"
varNC_sm = vName+"_anom_sm"
dimNames = (/"time", "lat", "lon"/)
dimSizes = (/ -1 , nlat, mlon/)
dimUnlim = (/ True , False, False/)
filedimdef(fnc,dimNames,dimSizes,dimUnlim)
filevardef(fnc, "time" ,typeof(time),getvardims(time))
filevardef(fnc, "lat" ,typeof(lat) ,getvardims(lat))
filevardef(fnc, "lon" ,typeof(lon) ,getvardims(lon))
filevardef(fnc, varNC ,typeof(xAnom) ,getvardims(xAnom))
filevardef(fnc, varNC_sm,typeof(xAnom) ,getvardims(xAnom))
filevarattdef(fnc,"time" ,time) ; copy time attributes
filevarattdef(fnc,"lat" ,lat) ; copy lat attributes
filevarattdef(fnc,"lon" ,lon) ; copy lon attributes
filevarattdef(fnc,varNC ,xAnom)
filevarattdef(fnc,varNC_sm,xAnom_sm)
fnc->time = (/time/)
fnc->lat = (/lat/)
fnc->lon = (/lon/)
fnc->$varNC$ = (/xAnom /)
fnc->$varNC_sm$ = (/xAnom_sm/)
end if
;************************************************
; plotting parameters
;************************************************
if (PLOT) then
LAT = (/ 60, 45, 5, -5, -45, 60 /)
LON = (/270, 30, 90, 90, 180, 0 /)
nPts = dimsizes( LAT )
plot = new ( nPts, graphic)
data = new ( (/2,366/), typeof(xClmDay), getFillValue(xClmDay))
wks = gsn_open_wks (wksType,wksName)
res = True ; plot mods desired
res@gsnDraw = False
res@gsnFrame = False
res@trXMinF = 1
res@trXMaxF = 366
;res@tiMainString = ""
res@xyLineThicknesses = (/1.0, 2.0/) ; make 2nd lines thicker
res@xyLineColors = (/"blue","red"/) ; change line color
res@xyMonoDashPattern = True ; all solid
do np=0,nPts-1
data(0,:) = xClmDay(:,{LAT(np)},{LON(np)})
data(1,:) = xClmDay_sm(:,{LAT(np)},{LON(np)})
res@gsnCenterString = "lat="+LAT(np)+" lon="+LON(np)
plot(np) = gsn_csm_y (wks,data,res) ; create plot
end do
resP = True ; modify the panel plot
resP@gsnPanelMainString = vName+": Raw/Smooth Daily Climatology: "+yrStrt+"-"+yrLast
resP@gsnMaximize = True
resP@gsnPaperOrientation = "portrait"
gsn_panel(wks,plot,(/(nPts/2),2/),resP) ; now draw as one plot
;==========
; Plot anomalies for an arbitrarily selected near equatorial location
; Time: Oct 1, 1996 to April 1,1997 [arbitrary selection]
;==========
LATX = 0
LONX = 90
yyyymmdd = cd_calendar(time, -2)
;;yrfrac = yyyymmdd_to_yyyyfrac (yyyymmdd, 0)
;;delete(yrfrac@long_name)
xAnom@long_name = "Anomalies from Raw" ; short labels for plot
xAnom_sm@long_name = "Anomalies from Smooth"
ntBegin = ind(yyyymmdd.eq.19961001)
ntEnd = ind(yyyymmdd.eq.19970401)
monthLabels = (/1,4,7,10/)
monNam = (/"Jan","Feb","Mar","Apr","May","Jun" \
,"Jul","Aug","Sep","Oct","Nov","Dec" /)
xVal = new(ntim, typeof(xAnom&time) , "No_FillValue") ; bigger than
xLab = new(ntim, "string", "No_FillValue") ; needed
xValm = new(ntim, typeof(xAnom&time) , "No_FillValue") ; bigger than
ntm = -1
cr = inttochar(10) ; carriage return
do nt=ntBegin,ntEnd
if (day(nt).eq.1) then
ntm = ntm + 1
xVal(ntm) = xAnom&time(nt)
xLab(ntm) = monNam(month(nt)-1)
if (month(nt).eq.1) then
xLab(ntm) = xLab(ntm) + cr +sprinti("%0.4i", year(nt))
end if
end if
end do
rxy = True
rxy@gsnDraw = False
rxy@gsnFrame = False
rxy@gsnYRefLine = 0.0 ; create a reference line
rxy@gsnAboveYRefLineColor = "red" ; above ref line fill red
rxy@gsnBelowYRefLineColor = "blue" ; below ref line fill blue
rxy@xyLineThicknessF = 2.0
rxy@vpHeightF = 0.4 ; resize
rxy@vpWidthF = 0.8
rxy@tmXBMode = "Explicit"
rxy@tmXBValues = xVal(0:ntm)
rxy@tmXBLabels = xLab(0:ntm)
plot(0) = gsn_csm_xy (wks,time(ntBegin:ntEnd) \
,xAnom(ntBegin:ntEnd,{0},{90}),rxy)
plot(1) = gsn_csm_xy (wks,time(ntBegin:ntEnd) \
,xAnom_sm(ntBegin:ntEnd,{0},{90}),rxy)
resP@gsnPanelMainString = vName+": Anomalies: (0,90E)"
gsn_panel(wks,plot(0:1),(/2,1/),resP)
end if
end三、将nc文件进行插值
由于我的uv文件lon lat精度是360*181，代码里需要的精度是144*181，所以需要对uv文件进行插值。我这里用的是ubuntu，已经装好了cdo，可以进行数据切片等操作。
先将nc文件网格格式generic改为lonlat：
1、创建一个新的网格文件，将“generic”替换为“lonlat”
cdo griddes 你的文件名称.nc > mygrid
sed -i "s/generic/lonlat/g" mygrid
2、使用cdo命令重新设置网格
cdo setgrid,mygrid 你的文件名.nc 输出文件名.nc
也可解决 cdo转换经度-180~180 为0~360时出现报错：cdo sellonlatbox: Unsupported grid type: generic的问题
——————————
更改分辨率 cdo remapbil,r100x100 input.nc output.nc
第一个100表示lon，第二个表示lat，这里改成r144x181
四、mjoclivar14
进行eof分解，所需数据为olr异常场，插值后的u850和u200异常场，结果会输出一个新的nc文件，需要在下个部分用上
;******************************************************
;
; mjoclivar_14.ncl
;
;***********************************************************
; Combined EOFs
; Latest Update: July, 2016: Eun-Pa Lim; Bureau of Meteorology, Australia
;***********************************************************
;;
;; The following are automatically loaded from 6.2.0 onward
;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
undef("read_rename")
function read_rename(f[1]:file, varName[1]:string \
,iStrt[1]:integer, iLast[1]:integer \
,latS[1]:numeric , latN[1]:numeric )
; Utility to force specific named dimensions
; This is done for historical reasons (convenience)
begin
work = f->$varName$(iStrt:iLast,{latS:latN},:) ; (time,lat,lon)
work!0 = "time" ; CAM model names
work!1 = "lat"
work!2 = "lon"
return(work)
end
; =========================> MAIN <==============================
begin
neof = 2
latS = -15
latN = 15
ymdStrt = 20000101 ; start yyyymmdd
ymdLast = 20191231 ; last
yrStrt = ymdStrt/10000
yrLast = ymdLast/10000
pltDir = "/mnt/c/project/MJO/mjo14/"; plot directory
pltType = "png"
pltName = "mjoclivar"
diri = "/mnt/c/project/MJO/data/"; input directory
filolr = "OLR.day.anomalies.2000-2019.nc"
filu200 = "U200.day.anomalies.2000-2019.144.nc"
filu850 = "U850.day.anomalies.2000-2019.144.nc"
;************************************************
; create BandPass Filter
;************************************************
ihp = 2 ; bpf=>band pass filter
nWgt = 201
sigma = 1.0 ; Lanczos sigma
fca = 1./100.
fcb = 1./20.
wgt = filwgts_lanczos (nWgt, ihp, fca, fcb, sigma );计算lanczos滤波器的权重系数
;***********************************************************
; Find the indices corresponding to the start/end times
;***********************************************************
f = addfile (diri+filolr , "r")
TIME = f->time ; days since ...
YMD = cd_calendar(TIME, -2) ; entire (time,6) ；转换为公历日期
iStrt = ind(YMD.eq.ymdStrt) ; index start
iLast = ind(YMD.eq.ymdLast) ; index last
delete([/ TIME, YMD /])
;***********************************************************
; Read anomalies
;***********************************************************
work0 = read_rename(f,"OLR_anom",iStrt,iLast,latS,latN) ; (time,lat,lon)
OLR = dim_avg_n_Wrap(work0, 1) ; (time,lon)；计算数组work0第1维的平均值
f = addfile (diri+filu850 , "r")
work1 = read_rename(f,"U850_anom",iStrt,iLast,latS,latN) ; (time,lat,lon)
U850 = dim_avg_n_Wrap(work1, 1) ; (time,lon)
;计算数组work1第1维的平均值
f = addfile (diri+filu200 , "r")
work2 = read_rename(f,"U200_anom",iStrt,iLast,latS,latN) ; (time,lat,lon)
U200 = dim_avg_n_Wrap(work2, 1) ; (time,lon)
dimw = dimsizes( work0 ) ;返回维数
ntim = dimw(0) ;第0维的维数 1826
nlat = dimw(1) ;第1维的维数 13
mlon = dimw(2) ;第2维的维数 144
;mlon=144;mlon和olr维数需要一致，手动设置为144
;print(mlon)
;printVarSummary(work2)
;printVarSummary(work0)
delete(work0)
;printVarSummary(ntim)
lon = OLR&lon
time = OLR&time
date = cd_calendar(time, -2) ; yyyymmdd
;printVarSummary(lon)
;************************************************
; Apply the band pass filter to the original anomalies
;************************************************
olr = wgt_runave_n_Wrap ( OLR, wgt, 0, 0) ; (time,lon)；对数组olr第0维以wgt为权重进行滑动平均
u850 = wgt_runave_n_Wrap (U850, wgt, 0, 0)
u200 = wgt_runave_n_Wrap (U200, wgt, 0, 0)
;************************************************
; remove temporal means of band pass series: *not* necessary ；删除带通序列的时间平均值
;************************************************
olr = dim_rmvmean_n( olr, 0) ; (time,lon)；计算olr数组第0维的距平
u850 = dim_rmvmean_n(u850, 0) ; u850
u200 = dim_rmvmean_n(u200, 0) ; u200
;************************************************
; Compute the temporal variance at each lon；计算每个lon的时间方差
;************************************************
var_olr = dim_variance_n_Wrap( olr, 0) ; (lon)
var_u850 = dim_variance_n_Wrap(u850, 0)
var_u200 = dim_variance_n_Wrap(u200, 0)
;************************************************
; Compute the zonal mean of the temporal variance
;************************************************
zavg_var_olr = dim_avg_n_Wrap( var_olr , 0) ;计算时间方差的平均值
zavg_var_u850 = dim_avg_n_Wrap( var_u850, 0)
zavg_var_u200 = dim_avg_n_Wrap( var_u200, 0)
;************************************************
; Normalize by sqrt(avg_var*)；规范化：开平方
;************************************************
olr = olr/sqrt(zavg_var_olr ) ; (time,lon)
u850 = u850/sqrt(zavg_var_u850)
u200 = u200/sqrt(zavg_var_u200)
printVarSummary(olr)
;************************************************
; Combine the normalized data into one variable；将规范化的数据合并为一个变量
;************************************************
cdata = new ( (/3*mlon,ntim/), typeof(olr), getFillValue(olr));new（维数，数值类型，缺测值的数值）
do ml=0,mlon-1
cdata(ml ,:) = (/ olr(:,ml) /)
cdata(ml+ mlon,:) = (/ u850(:,ml) /)
cdata(ml+2*mlon,:) = (/ u200(:,ml) /)
end do
;************************************************
; Compute **combined** EOF; Sign of EOF is arbitrary；计算**组合**EOF；EOF的符号是任意的
;************************************************
eof_cdata = eofunc(cdata , neof, False) ; (neof,3*mlon)；经验函数正交分解？neof之前赋值为2；计算前2个主要模态
print("==============")
printVarSummary(eof_cdata)
printMinMax(eof_cdata, True)
eof_ts_cdata = eofunc_ts(cdata,eof_cdata,False) ; (neof,3*ntim)；再计算前两个模态所对应的时间序列
print("==============")
printVarSummary(eof_ts_cdata)
printMinMax(eof_ts_cdata, True)
;************************************************
; For clarity, explicitly extract each variable. Create time series ；显式提取每个变量。创建时间序列
;************************************************
nvar = 3 ; "olr", "u850", "u200"
ceof = new( (/nvar,neof,mlon/), typeof(cdata), getFillValue(cdata))
do n=0,neof-1
ceof(0,n,:) = eof_cdata(n,0:mlon-1) ; olr
ceof(1,n,:) = eof_cdata(n,mlon:2*mlon-1) ; u850
ceof(2,n,:) = eof_cdata(n,2*mlon:) ; u200
end do
ceof!0 = "var" ;ceof的第零维命名为var
ceof!1 = "eof"
ceof!2 = "lon"
;dimc=dimsizes(ceof)
;dierwei=dimc(2)
;print(dierwei)
;print(lon)
;printVarSummary(olr)
ceof&lon = olr&lon ;ceof的lon坐标变量
ceof_ts = new( (/nvar,neof,ntim/), typeof(cdata), getFillValue(cdata))
ceof_ts(0,:,:) = eofunc_ts_Wrap( olr(lon|:,time|:),ceof(0,:,:),False) ; (0,neof,ntim)
printVarSummary(olr)
printVarSummary(ceof)
printVarSummary(ceof_ts)
printVarSummary(u850)
printVarSummary(u200)
ceof_ts(1,:,:) = eofunc_ts_Wrap(u850(lon|:,time|:),ceof(1,:,:),False) ; (1,neof,ntim)
ceof_ts(2,:,:) = eofunc_ts_Wrap(u200(lon|:,time|:),ceof(2,:,:),False) ; (2,neof,ntim)
;计算模态的时间序列
;**********************************************t*
; Add code contributed by Marcus N. Morgan, Florida Institute of Technology; Feb 2015
; Calculate % variance (pcv_ )accounted for by OLR, U850 and U200 ；计算olr u850 u200引起的方差百分比
;************************************************
pcv_eof_olr = new(neof,typeof(ceof))
pcv_eof_u850 = new(neof,typeof(ceof))
pcv_eof_u200 = new(neof,typeof(ceof))
do n=0,neof-1
pcv_eof_olr(n) = avg((ceof(0,n,:)*sqrt(ceof@eval(n)))^2)*100
pcv_eof_u850(n) = avg((ceof(1,n,:)*sqrt(ceof@eval(n)))^2)*100
pcv_eof_u200(n) = avg((ceof(2,n,:)*sqrt(ceof@eval(n)))^2)*100
;;print("pcv: neof="+(n+1)+": "+pcv_eof_olr(n)+" "+pcv_eof_u850(n)+" "+pcv_eof_u200(n))
end do
;************************************************
; Change sign of EOFs for spatial structures of PC1 and PC2
; to represent convection over the tropical Indian Ocean and the tropical western Pacific Ocean, respectively
; (Ad hoc approach)
;************************************************
imax_olr_eof1 = maxind(ceof(0,0,:)) ;返回数组ceof中第一个最大值的下标变量 (?)
imax_olr_eof2 = maxind(ceof(0,1,:))
lonmax_eof1 = ceof&lon(imax_olr_eof1) ; longitude of max value (i.e. suppressed convection)
lonmax_eof2 = ceof&lon(imax_olr_eof2)
if (lonmax_eof1.ge.60 .and. lonmax_eof1.lt.180) then ; Change the sign of EOF1
ceof(:,0,:) = -ceof(:,0,:) ; if OLR is positive
ceof_ts(:,0,:) = -ceof_ts(:,0,:) ; over the tropical Indian Ocean
eof_cdata(0,:) = -eof_cdata(0,:)
eof_ts_cdata(0,:) = -eof_ts_cdata(0,:)
end if
if (lonmax_eof2.ge.120 .and. lonmax_eof2.lt.180) then ; Change the sign of EOF2
ceof(:,1,:) = -ceof(:,1,:) ; if OLR is positive
ceof_ts(:,1,:) = -ceof_ts(:,1,:) ; over the tropical western Pacific Ocean
eof_cdata(1,:) = -eof_cdata(1,:)
eof_ts_cdata(1,:) = -eof_ts_cdata(1,:)
end if
print("==============")
printVarSummary(eof_cdata)
printMinMax(eof_cdata, True)
;************************************************
; Compute cross correlation of each variable's EOF time series at zero-lag :计算零滞后条件下每个变量EOF时间序列的互相关
;************************************************
r_olr_u850 = escorc(ceof_ts(0,:,:) , ceof_ts(1,:,:) ) ; (neof)
r_olr_u200 = escorc(ceof_ts(0,:,:) , ceof_ts(2,:,:) ) ;计算x最右边维和y最右边维的同时线性相关系数
r_u850_u200 = escorc(ceof_ts(1,:,:) , ceof_ts(2,:,:) )
print("==============")
do n=0,neof-1
print("neof="+n \
+" r_olr_u850=" +sprintf("%4.3f",r_olr_u850(n)) \
+" r_olr_u200=" +sprintf("%4.3f",r_olr_u200(n)) \
+" r_u850_u200="+sprintf("%4.3f",r_u850_u200(n)) )
end do
print("==============")
;************************************************
; Compute cross correlation of the multivariate EOF; EOF 1 vs EOF 2 ;计算多元EOF的互相关；EOF 1与EOF 2
;************************************************
mxlag = 25
rlag_01 = esccr(eof_ts_cdata(0,:),eof_ts_cdata(1,:), mxlag) ; (N,mxlag+1) 计算x最右边维和y最右边维的交叉相关系数
rlag_10 = esccr(eof_ts_cdata(1,:),eof_ts_cdata(0,:), mxlag) ; (N,mxlag+1)
ccr_12 = new ( (/2*mxlag+1/), float)
ccr_12(mxlag:) = rlag_10(0:mxlag)
ccr_12(0:mxlag) = rlag_01(::-1) ; reverse order
;;print(ccr_12)
;************************************************
; Normalize the multivariate EOF 1&2 component time series 归一化多元EOF 1和2分量时间序列
; Compute (PC1^2+PC2^2): values > 1 indicate "strong" periods;计算（PC1^2+PC2^2）：值>1表示“强”周期
;************************************************
eof_ts_cdata(0,:) = eof_ts_cdata(0,:)/stddev(eof_ts_cdata(0,:));stddev:计算数组x中所有数值的标准差
eof_ts_cdata(1,:) = eof_ts_cdata(1,:)/stddev(eof_ts_cdata(1,:))
mjo_ts_index = eof_ts_cdata(0,:)^2 + eof_ts_cdata(1,:)^2
mjo_ts_index_smt = runave(mjo_ts_index, 91, 0) ; 91-day running mean 对数组x的第0维等权重滑动平均
nGood = num(.not.ismissing(mjo_ts_index)) ; # non-missing
nStrong = num(mjo_ts_index .ge. 1.0)
print("nGood="+nGood+" nStrong="+nStrong+" nOther="+(nGood-nStrong))
;************************************************
; Write PC results to netCDF for use in another example.将PC结果写入netCDF，以便在另一个示例中使用
;************************************************
mjo_ts_index!0 = "time"
mjo_ts_index&time = time
mjo_ts_index@long_name = "MJO PC INDEX"
mjo_ts_index@info = "(PC1^2 + PC2^2)"
PC1 = eof_ts_cdata(0,:)
PC1!0 = "time"
PC1&time = time
PC1@long_name = "PC1"
PC1@info = "PC1/stddev(PC1)"
PC2 = eof_ts_cdata(1,:)
PC2!0 = "time"
PC2&time = time
PC2@long_name = "PC2"
PC2@info = "PC2/stddev(PC2)"
diro = "/mnt/c/project/MJO/mjo14/"
filo = "MJO_PC_INDEX.nc"
system("/bin/rm -f "+diro+filo) ; remove any pre-existing file 删除任何预先存在的文件
ncdf = addfile(diro+filo,"c") ; open output netCDF file
; make time an UNLIMITED dimension 让时间成为无限的维度
filedimdef(ncdf,"time",-1,True) ; recommended for most applications
; output variables directly 直接输出变量
ncdf->MJO_INDEX = mjo_ts_index
ncdf->PC1 = PC1
ncdf->PC2 = PC2
;------------------------------------------------------------
; PLOTS ;绘图
;------------------------------------------------------------
yyyymmdd = cd_calendar(time, -2)
yrfrac = yyyymmdd_to_yyyyfrac(yyyymmdd, 0.0)
delete([/ yrfrac@long_name, lon@long_name /])
day = ispan(-mxlag, mxlag, 1) ;创造一维的等差数列，公差为1，/-25 -24...24 25/
;day@long_name = "lag (day)"
pltPath = pltDir+pltName
wks = gsn_open_wks(pltType,pltPath) ;建立一个图形文件
gsn_define_colormap(wks,"default") ;更换默认色板
plot = new(3,graphic)
;************************************************
; Multivariate EOF plots 多变量EOF图
;************************************************
rts = True
rts@gsnDraw = False ; don't draw yet
rts@gsnFrame = False ; don't advance frame yet
rts@gsnScale = True ; force text scaling
rts@vpHeightF = 0.40 ; Changes the aspect ratio
rts@vpWidthF = 0.85
rts@vpXF = 0.10 ; change start locations
rts@vpYF = 0.75 ; the plot
rts@xyLineThicknesses = (/2, 2, 2/)
rts@xyLineColors = (/"black","red","blue"/)
rts@gsnYRefLine = 0. ; reference line
rts@trXMaxF = max(lon)
rts@trXMinF = min(lon)
rts@pmLegendDisplayMode = "Always" ; turn on legend
rts@pmLegendSide = "Top" ; Change location of
rts@pmLegendParallelPosF = 1.16 ; move units right
rts@pmLegendOrthogonalPosF = -0.50 ; move units down
rts@pmLegendWidthF = 0.15 ; Change width and
rts@pmLegendHeightF = 0.15 ; height of legend.
rts@lgLabelFontHeightF = 0.0175
rtsP = True ; modify the panel plot
; rtsP@gsnMaximize = True ; large format
rtsP@gsnPanelMainString = "Multivariate EOF: 15S-15N: "+yrStrt+"-"+yrLast
do n=0,neof-1
rts@xyExplicitLegendLabels = (/"OLR: "+sprintf("%4.1f", pcv_eof_u200(n)) +"%" \
,"U850: "+sprintf("%4.1f", pcv_eof_u850(n))+"%" \
,"U200: "+sprintf("%4.1f", pcv_eof_olr(n))+"%" /)
rts@gsnLeftString = "EOF "+(n+1)
rts@gsnRightString = sprintf("%3.1f",ceof@pcvar(n)) +"%"
plot(n) = gsn_csm_xy (wks,lon,ceof(:,n,:),rts)
end do
gsn_panel(wks,plot(0:1),(/2,1/),rtsP) ; now draw as one plot
;-----------------------------------------
; The following doesn't work with some older versions of NCL 以下内容不适用于一些旧版本的NCL
; With old versions, the user must delete each individually. 对于旧版本，用户必须单独删除每个版本
;-----------------------------------------
delete([/ rts@xyExplicitLegendLabels, rts@pmLegendDisplayMode \
, rts@xyLineThicknesses , rts@gsnLeftString \
, rts@gsnRightString , rts@xyLineColors \
, rts@trXMaxF , rts@trXMinF /] )
lag = ispan(-mxlag,mxlag,1)
lag@long_name = "lag (days)"
plot(0) = gsn_csm_xy (wks, lag ,ccr_12,rts)
rtsP@gsnPanelMainString = "Cross Correlation: Multivariate EOF: 15S-15N: " \
+ yrStrt+"-"+yrLast
rtsP@gsnPaperOrientation = "portrait" ; force portrait
gsn_panel(wks,plot(0),(/1,1/),rtsP) ; now draw as one plot
;************************************************
; MJO "strong" index
;************************************************
rts@gsnYRefLine = 1.0
rts@gsnYRefLineColor = "black"
rts@xyMonoDashPattern = True
rts@xyLineColors = (/"black", "blue"/)
rts@xyLineThicknesses = (/1, 2/)
rts@pmLegendDisplayMode = "Always" ; turn on legend
rts@pmLegendWidthF = 0.12 ; Change width and
rts@pmLegendHeightF = 0.10 ; height of legend.
rts@pmLegendParallelPosF = 0.86 ; move units right
rts@pmLegendOrthogonalPosF = -0.40 ; move units down
rts@xyExplicitLegendLabels = (/"daily", "91-day runavg" /)
mjo_ind_plt = new ( (/2,ntim/), typeof(mjo_ts_index))
mjo_ind_plt(0,:) = mjo_ts_index
mjo_ind_plt(1,:) = (/ mjo_ts_index_smt /)
plot(0) = gsn_csm_xy(wks, yrfrac,mjo_ind_plt,rts)
rtsP@gsnPanelMainString = "MJO Index: (PC1^2+ PC2^2) : 15S-15N: "+yrStrt+"-"+yrLast
gsn_panel(wks,plot(0),(/1,1/),rtsP) ; now draw as one plot
end五、mjoclivar16
提取位相，所需数据为olr异常场，插值后的u850异常场和v850异常场，最后得到的即文章开始的两张位相图
;***********************************************************
;
; mjoclivar_16.ncl
;
;***********************************************************
; Generate life cycle composites based upon daily phase space
; If the MJO_INDEX is < 1.0 it is not included
;
; Source: Eun-Pa Lim: Bureau of Meteorology, Australia
; July, 2016
;***********************************************************
;
; These files are loaded by default in NCL V6.2.0 and newer
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
latS = -20
latN = 20
ymdStrt = 20000101 ; start yyyymmdd
ymdLast = 20191231 ; last
yrStrt = ymdStrt/10000
yrLast = ymdLast/10000
pltSubTitle = "Anomalous: OLR, U850, V850"
pltDir= "/mnt/c/project/MJO/mjo16/" ; plot directory
pltType= "png" ; send graphics to PNG file
pltName= "mjoclivar" ; yrStrt+"_"+yrLast
diro= "/mnt/c/project/MJO/data/" ; output directory
diri= "/mnt/c/project/MJO/data/" ; input directory
filo="OLR.day.anomalies.2000-2019.nc"
filu="U850.day.anomalies.2000-2019.144.nc"
filv="V850.day.anomalies.2000-2019.144.nc"
;************************************************
; create BandPass Filter
;************************************************
ihp = 2 ; bpf=>band pass filter
nWgt = 201
sigma = 1.0 ; Lanczos sigma
fca = 1./100.
fcb = 1./20.
wgt = filwgts_lanczos (nWgt, ihp, fca, fcb, sigma );计算滤波器的权重系数
;***********************************************************
; Find the indicies (subscripts) corresponding to the start/end times 查找与开始/结束时间相对应的标记（下标）
;***********************************************************
f = addfile (diri+filu , "r")
TIME = f->time ; days since ...
YMD = cd_calendar(TIME, -2) ; entire (time,6)
iStrt = ind(YMD.eq.ymdStrt) ; index start
iLast = ind(YMD.eq.ymdLast) ; index last
delete(TIME)
delete(YMD )
time = f->time(iStrt:iLast) ; days since ...
u = f->U850_anom(iStrt:iLast,{latS:latN},:)
;***********************************************************
; Read anomalies frpm other fields 读取其他字段的异常
;***********************************************************
f = addfile (diri+filv , "r")
v = f->V850_anom(iStrt:iLast,{latS:latN},:)
f = addfile (diri+filo , "r")
x = f->OLR_anom(iStrt:iLast,{latS:latN},:)
dimx = dimsizes( x )
ntim = dimx(0)
nlat = dimx(1)
mlon = dimx(2)
;************************************************
; Apply the band pass filter to the original anomalies
;************************************************
a = wgt_runave_leftdim (x, wgt, 0);对数组x第0维滑动平均？（leftdim是啥
u = wgt_runave_leftdim (u, wgt, 0)
v = wgt_runave_leftdim (v, wgt, 0)
;***********************************************************
; Open PC components file created in 'mjo_14.ncl'
;***********************************************************
dirMJO = "/mnt/c/project/MJO/mjo14/" ; input directory
fMJO = "MJO_PC_INDEX.nc" ; created in mjo_14.ncl
f = addfile (dirMJO+fMJO, "r")
;***********************************************************
; Find the indices corresponding to the start/end times 查找与开始/结束时间对应的索引
;***********************************************************
TIME = f->time ; days since ...
YMD = cd_calendar(TIME, -2) ; entire (time,6)
iStrt = ind(YMD.eq.ymdStrt) ; index start
iLast = ind(YMD.eq.ymdLast) ; index last
delete(TIME)
delete(YMD )
delete(time)
;***********************************************************
; Read the data for the desired period 读取所需时间段的数据
;***********************************************************
pc1 = f->PC1(iStrt:iLast)
pc2 = f->PC2(iStrt:iLast)
mjo_indx= f->MJO_INDEX(iStrt:iLast)
time = pc1&time
ymdhms = cd_calendar(time, 0)
imon = floattoint( ymdhms(:,1) ) ; convenience
iday = floattoint( ymdhms(:,2) ) ; subscripts must be integer
;***********************************************************
; Place each array into an appropriate array 将每个阵列放入适当的阵列中
;***********************************************************
nPhase = 8
angBnd = new( (/2,nPhase/), "float")
angBnd(0,:) = fspan( 0,315,nPhase)
angBnd(1,:) = fspan( 45,360,nPhase)
r2d = 180./(4.*atan(1.0))
ang = atan2(pc2,pc1)*r2d ; phase space
nn = ind(ang.lt.0)
ang(nn) = ang(nn) + 360 ; make 0 to 360
;----------------------------------------------------------
; 0 <= ang < 45 --> MJO Phase 5 (i.e. +ve PC1 & +ve PC2)
;
; print(pc1(:19)+" "+pc2(:19)+" "+ang(:19))
;----------------------------------------------------------
nDays = new (nPhase, "integer")
pLabel = "P"+ispan(1,nPhase,1)+": "
;------------------------------------------------------------
; PLOTS
;------------------------------------------------------------
; pltPath = pltDir+pltName+"_rv_16"
pltPath = pltDir+pltName
wks = gsn_open_wks(pltType,pltPath)
plot = new(nPhase,graphic) ; create graphic array
res = True
res@gsnDraw = False ; don't draw yet
res@gsnFrame = False ; don't advance frame yet
res@mpFillOn = False ; turn off map fill
res@mpMinLatF = latS ; zoom in on map
res@mpMaxLatF = latN
res@mpCenterLonF = 210.
res@cnFillOn = True ; turn on color fill
res@cnFillPalette = "ViBlGrWhYeOrRe" ; set color map
res@cnLinesOn = False ; True is default
res@cnLineLabelsOn = False ; True is default
res@lbLabelBarOn = False ; turn off individual lb's
res@gsnScalarContour = True ; contour 3rd array
res@gsnMajorLatSpacing = 15
res@gsnMajorLonSpacing = 60
res@tmXBLabelFontHeightF = 0.01
res@tmYLLabelFontHeightF = 0.01
; common contours
;mnmxint = nice_mnmxintvl( min(x) , max(x), 16, False)
res@cnLevelSelectionMode = "ManualLevels"
res@cnMinLevelValF = -40 ; -100; mnmxint(0)
res@cnMaxLevelValF = 40 ; 80; mnmxint(1)
res@cnLevelSpacingF = 5 ; 20; mnmxint(2)
;print(res)
res@vcMinDistanceF = 0.01 ; thin the vector density
res@vcRefMagnitudeF = 2.0 ; define vector ref mag
res@vcRefLengthF = 0.025 ; define length of vec ref
res@vcRefAnnoOrthogonalPosF = -1.0 ; move ref vector
res@vcRefAnnoArrowLineColor = "black" ; change ref vector color
res@vcRefAnnoArrowUseVecColor = False ; don't use vec color for ref
; panel plot only resources
resP = True ; modify the panel plot
resP@gsnMaximize = True ; large format
resP@gsnPanelLabelBar = True ; add common colorbar
resP@lbLabelFontHeightF = 0.01
resP@gsnPanelBottom = 0.05 ; add some space at bottom
resP@pmLabelBarWidthF = 0.8 ; label bar width
resP@pmLabelBarHeightF = 0.05
resP@gsnPanelFigureStringsFontHeightF = 0.0125 ; bit larger than default
;resP@pmLabelBarOrthogonalPosF = 0.015 ; move labelbar up a bit
txres = True
txres@txFontHeightF = 0.01
txid = gsn_create_text(wks, pltSubTitle, txres)
amres = True
;amres@amParallelPosF = 0.575
amres@amOrthogonalPosF = 0.75
amres@amJust = "CenterCenter"
;amres@amResizeNotify = True
;*******************************************
; Loop over each phase
;*******************************************
res@gsnLeftString = ""
res@gsnRightString = ""
do nSeason=1,2
if (nSeason.eq.1) then
resP@gsnPanelMainString = yrStrt+"-"+yrLast+": May to Oct"
else
resP@gsnPanelMainString = yrStrt+"-"+yrLast+": Nov to Apr"
end if
do n=0,nPhase-1
na = n+4 ; temporary adjustment for 0 <= ang < 45 represents MJO phase 5 not MJO phase 1
if(na.gt.7) then
na = na - 8
end if
; print(na)
if (nSeason.eq.1) then
nt = ind(mjo_indx.gt.1.0 .and. \
(imon.ge.5 .and. imon.le.10).and. \
ang.ge.angBnd(0,n) .and. ang.lt.angBnd(1,n))
else
nt = ind(mjo_indx.gt.1.0 .and. \
(imon.ge.11 .or. imon.le. 4).and. \
ang.ge.angBnd(0,n) .and. ang.lt.angBnd(1,n))
end if
if (.not.all(ismissing(nt))) then
xAvg = dim_avg_Wrap( x(lat|:,lon|:,time|nt) )
uAvg = dim_avg_Wrap( u(lat|:,lon|:,time|nt) )
vAvg = dim_avg_Wrap( v(lat|:,lon|:,time|nt) )
nDays(na) = dimsizes(nt)
res@tmXBLabelsOn = False ; do not draw lon labels
res@tmXBOn = False ; lon tickmarks
if (n.eq.(nPhase-1)) then ;
res@tmXBLabelsOn = True ; draw lon labels
res@tmXBOn = True ; tickmarks
end if
plot(na) = gsn_csm_vector_scalar_map(wks,uAvg,vAvg,xAvg,res)
end if
delete(nt) ; will change next iteration
end do
resP@gsnPanelFigureStrings= pLabel+nDays
gsn_panel(wks,plot,(/nPhase,1/),resP) ; now draw as one plot
end do
end
1131
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
