- 积分
- 318
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2021-3-21
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
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 = -40.
latN = 40.
lonL = 0.
lonR = 180.
neof=2
ymdStrt = 19810101 ; start yyyymmdd
ymdLast = 20101231 ; last
yrStrt = ymdStrt/10000
yrLast = ymdLast/10000
pltDir = "./" ; plot directory
pltType = "ps"
pltName = "monsoon_mveof"
diri = "./" ; input directory
filolr = "olr.day.anomalies.1981-2010.nc"
;filu850 = "uwnd.day.850.anomalies.1981-2010.nc"
;filv850 = "vwnd.day.850.anomalies.1981-2010.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 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
;***********************************************************
;f =addfile("olr.day.anomalies.1981-2010.nc","r")
fu=addfile("uwnd.day.850.anomalies.1981-2010.nc","r")
fv=addfile("vwnd.day.850.anomalies.1981-2010.nc","r")
TIME = f->time
olr=f->olr({iStrt:iLast},{latS:latN},{lonL:lonR})
olr=dim_standardize_n_Wrap(olr, 0, 0)
delete(TIME)
TIME = fu->time
uwnd=short2flt(fu->uwnd({iStrt:iLast},{latS:latN},{lonL:LonR}))
u=dim_standardize_n_Wrap(uwnd, 0, 0)
delete(TIME)
TIME = fv->time
vwnd=short2flt(fv->vwnd({iStrt:iLast},{latS:latN},{lonL:LonR}))
v=dim_standardize_n_Wrap(vwnd, 0, 0)
delete(TIME)
;*******************Write into one dataset*****************************
; Combine the normalized data into one variable ;
;**********************************************************************
dimw = dimsizes(olr)
mtim = dimw(0)
mlat = dimw(1)
mlon = dimw(2)
cdata = new ( (/3*mlat,3*mlon,mtim/), typeof(olr), getFillValue(olr))
do m1=0,mlat-1
do m2=0,mlon-1
cdata(m1, m2, :)=(/olr(:,m1,m2)/)
copy_VarMeta(olr,cdata(m1, m2, :) )
cdata(m1+mlat, m2+mlon, :)=(/u (:,m1,m2)/)
copy_VarMeta(olr, cdata(m1+mlat, m2+mlon, :))
cdata(m1+2*mlat,m2+2*mlon,:)=(/v (:,m1,m2)/)
copy_VarMeta(olr, cdata(m1+2*mlat,m2+2*mlon,:))
end do
end do
;***************************************t*****************************
; Compute **combined** EOF; Sign of EOF is arbitrary ;
;**********************************************************************
eof_cdata = eofunc(cdata , neof, False)
eof_ts_cdata = eofunc_ts(cdata,eof_cdata,False)
printVarSummary(eof_cdata)
printVarSummary(eof_ts_cdata)
nvar=3 ;olr,u,v
ceof = new( (/nvar,neof,mlat,mlon/), typeof(cdata), getFillValue(cdata))
do n=0,neof-1
ceof(0,n,:,:) = eof_cdata(n,0:mlat-1,0:mlon-1) ; olr
ceof(1,n,:,:) = eof_cdata(n,mlat:mlat*2-1,mlon:2*mlon-1) ; u
ceof(2,n,:,:) = eof_cdata(n,2*mlat:,2*mlon:) ; v
end do
ceof!0 = "var"
ceof!1 = "eof"
ceof!2 = "lat"
ceof!3 = "lon"
ceof&lat = olr&lat
ceof&lon = olr&lon
printVarSummary(ceof)
printMinMax(ceof(0,1,:,:),False)
;***************************************t*****************************
; plot now ;
;**********************************************************************
wks=gsn_open_wks("ps","MVEOF")
plot = new(neof,graphic)
plot2= new(neof,graphic)
gsn_define_colormap(wks,"MPL_RdYlGn")
res = True
res@gsnDraw = False ; don't draw yet
res@gsnFrame = False ; don't advance frame yet
; res@gsnPolar = "SH"
res@gsnAddCyclic = False
res@gsnMaximize = True
res@mpFillOn = True ; turn off map fill
res@mpOutlineOn = False
res@mpMinLatF = -40
res@mpMaxLatF = 40
res@mpMaxLonF = 180
res@mpMinLonF = 0
res@cnFillOn = True ; turn on color fill
;res@cnFillPalette = "BlueWhiteOrangeRed"
res@cnLinesOn = False ; True is default
res@cnLineLabelsOn = False ; True is default
res@lbLabelBarOn = False ; turn off individual lb's
;res@cnLevelSelectionMode = "ManualLevels"
;res@cnMaxLevelValF = -0.02
;res@cnMinLevelValF = 0.02
;res@cnLevelSpacingF = 0.005
res@cnLevelSelectionMode = "ExplicitLevels"
res@cnLevels = (/-0.02,-0.015,-0.01,-0.005,0,0.005,0.01,0.015,0.02/)
res@cnFillColors = (/5,25,38,56,0,0,80,94,108,121/)
; panel plot only resources
vcres = True
vcres@gsnDraw = False
vcres@gsnFrame = False
vcres@vcRefLengthF =0.045
vcres@vcMinDistanceF =0.017
;vcres@vcGlyphStyle = "CurlyVector"
resP = True ; modify the panel plot
resP@gsnMaximize = True ; large format
resP@gsnPanelLabelBar = True ; add common colorbar
resP@gsnPaperOrientation = "portrait" ; force portrait
resP@txString = ""
;*******************************************
; first plot
;*******************************************
do n=0,neof-1
res@tiMainString = "MVEOF "+(n+1)
res@gsnLeftString = "olr and 850hPa winds"
res@gsnRightString = ""
plot(n) =gsn_csm_contour_map(wks,ceof(0,n,:,:),res)
plot2(n)=gsn_csm_vector (wks,ceof(1,n,:,:), ceof(2,n,:,:), vcres)
overlay(plot(n),plot2(n))
draw(plot(n))
frame(wks)
end do
gsn_panel(wks,plot,(/neof,1/),resP) ; now draw as one plot
end
|
-
|