NCL初学者画海平面气压异常场

一、目的

计算并绘制2021年4月的海平面气压距平场（气候平均分别取为 1961-1990年（30a平均）、1991-2020年（30a平均）和1961- 2020年（60a平均）的结果）；

二、代码实现

1.数据提取并按要求切片

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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
;************************************************
;               求距平和气候态的代码
;************************************************

;************************************************
; define parameters
;************************************************
yrStrt                     = 1961
yrLast                     = 2020
nyr                        = yrLast-yrStrt+1
nmon                       = 12
ntime                      = nmon*nyr
lat1                       = 0
lat2                       = 90
;************************************************
;             1961-2020
;************************************************
fdir1                      = "/home/22597/hwmf/"
fname1                     = "slp.mon.mean.nc"
infile1                    = addfile(""+fdir1+""+fname1+"", "r")
TIME                       = infile1->time
slp = infile1->slp
YYYY                       = cd_calendar(TIME,-1)/100
iYYYY                      = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
;print(iYYYY)
XX_obs                     = infile1->slp(iYYYY,::-1,:)
XX_ob                      = XX_obs(3:ntime-1:12,{lat1:lat2},:)
f_avg                      = dim_avg_n_Wrap(XX_ob,0)
;************************************************
;             1961-1990
;************************************************
yrStrt0                     = 1961
yrLast0                     = 1990
nyr0                        = yrLast0-yrStrt0+1
nmon0                       = 12
ntime0                      = nmon0*nyr0
nyr0                        = yrLast0-yrStrt0+1
TIME0                       = infile1->time
YYYY0                       = cd_calendar(TIME0,-1)/100
iYYYY0                      = ind(YYYY0.ge.yrStrt0 .and. YYYY0.le.yrLast0)
XX_obs0                     = infile1->slp(iYYYY0,::-1,:)
XX_ob0                      = XX_obs0(3:ntime0-1:12,{lat1:lat2},:)
f_avg0                      = dim_avg_n_Wrap(XX_ob0,0)
;************************************************
;             1991-2020
;************************************************
yrStrt1                     = 1991
yrLast1                     = 2020
nyr1                        = yrLast1-yrStrt1+1
nmon1                       = 12
ntime1                      = nmon1*nyr1
nyr1                        = yrLast1-yrStrt1+1
TIME1                       = infile1->time
YYYY1                       = cd_calendar(TIME1,-1)/100
iYYYY1                      = ind(YYYY1.ge.yrStrt1 .and. YYYY1.le.yrLast1)
XX_obs1                     = infile1->slp(iYYYY1,::-1,:)
XX_ob1                      = XX_obs1(3:ntime1-1:12,{lat1:lat2},:)
f_avg1                      = dim_avg_n_Wrap(XX_ob1,0)
;****************************************************************
;                        求距平2021年4月
;****************************************************************
yrStrt7                     = 2020
yrLast7                     = 2021
nyr7                        = yrLast7-yrStrt7+1
nmon7                       = 12
ntime7                      = nmon7*nyr7
nyr7                        = yrLast7-yrStrt7+1
TIME7                       = infile1->time
YYYY7                       = cd_calendar(TIME7,-1)/100
iYYYY7                      = ind(YYYY7.ge.yrStrt7 .and. YYYY7.le.yrLast7)
XX_obs7                     = infile1->slp(iYYYY7,::-1,:)
XX_ob7                      = XX_obs7(3:ntime7-1:12,{lat1:lat2},:)

至此完成了数据的提取。

2.求距平

slp214 = XX_ob7(1, :, :)
anoslp6190 = slp214-f_avg0
anoslp6120 = slp214-f_avg
anoslp9120 = slp214-f_avg1

anoslp6190!0 = "lat" ;第一维的名称
anoslp6190&lat = f_avg0&lat ;第一维的值，直接从文件中的变量u中提取
anoslp6190!1 = "lon" ;第二维的名称
anoslp6190&lon = f_avg0&lon


anoslp6120!0 = "lat" ;第一维的名称
anoslp6120&lat = f_avg&lat ;第一维的值，直接从文件中的变量u中提取
anoslp6120!1 = "lon" ;第二维的名称
anoslp6120&lon = f_avg&lon

anoslp9120!0 = "lat" ;第一维的名称
anoslp9120&lat = f_avg1&lat ;第一维的值，直接从文件中的变量u中提取
anoslp9120!1 = "lon" ;第二维的名称
anoslp9120&lon = f_avg1&lon

3.由于数据被重新切片过，要赋予其新的经纬度属性。

slp214 = XX_ob7(1, :, :)
anoslp6190 = slp214-f_avg0
anoslp6120 = slp214-f_avg
anoslp9120 = slp214-f_avg1

anoslp6190!0 = "lat" ;第一维的名称
anoslp6190&lat = f_avg0&lat ;第一维的值，直接从文件中的变量u中提取
anoslp6190!1 = "lon" ;第二维的名称
anoslp6190&lon = f_avg0&lon


anoslp6120!0 = "lat" ;第一维的名称
anoslp6120&lat = f_avg&lat ;第一维的值，直接从文件中的变量u中提取
anoslp6120!1 = "lon" ;第二维的名称
anoslp6120&lon = f_avg&lon

anoslp9120!0 = "lat" ;第一维的名称
anoslp9120&lat = f_avg1&lat ;第一维的值，直接从文件中的变量u中提取
anoslp9120!1 = "lon" ;第二维的名称
anoslp9120&lon = f_avg1&lon

4.绘图程序

wks = gsn_open_wks("png","short-term climate forecast anomaly")     
gsn_define_colormap(wks,"BlueYellowRed")
res = True 
res@gsnAddCyclic = False
res@cnLevelSelectionMode = "ManualLevels"
res@cnMinLevelValF = -20.0 ; Min contour
res@cnMaxLevelValF = 20.0 ; Max contour
res@cnLevelSpacingF = 1 ; Spacing
res@cnFillOn = True ; Turn on contour fill
res@cnLinesOn = False ; Turn off contour lines
res@lbBoxLinesOn = False ; Turn off lbar box lines
res@mpFillOn                   = False
res@mpCenterLonF=180
res@mpMinLatF=0
res@mpMaxLatF=90
res@mpMinLonF=-180
res@mpMaxLonF=180
;res@tiMainString = "1991-2020 Northern hemisphere climate"  
plot = new(3,"graphic")
plot(0) = gsn_csm_contour_map(wks,anoslp6120,res)
plot(1) = gsn_csm_contour_map(wks,anoslp6190 ,res)
plot(2) = gsn_csm_contour_map(wks,anoslp9120,res)


; ************************************************
; create panel
; ************************************************
resP                  = True                   ; modify the panel plot
resP@gsnFrame         = False                  ; don't advance panel plot
; resP@gsnPanelLabelBar = True                   ; add common colorbar
resP@gsnPanelMainString ="short-term climate forecast"      ; new way of setting main title (as of NCL V6.4.0)
; resP@txString           = "A common title"     ; old way of setting main title
resP@gsnPanelBottom   = 0.05                   ; add space at bottom
resP@gsnPanelFigureStrings= (/"1961-2020 anomaly","1961-1990 anomaly","1991-2020 anomaly"/) ; add strings to panel
resP@amJust   = "TopLeft"
gsn_panel(wks,plot,(/3,1/),resP)  
frame(wks)

三、图片展示