NCL一次读取多个文件出多图

是冉冉升起的冉

begin
;指定年份1月份50 hPa 位势高度场
level=50;指定高度
; ;======================= PLOT ===============================

     res                                     = True
     res@gsnFrame                            = False
     res@gsnDraw                             = False
    ; res@vpWidthF                            = .6
     res@mpMinLatF=45
     ; ;------------------ Control the contour
     res@cnLinesOn                           =  False          ; Contour lines on/off: True/False     g
     res@cnFillOn                            =  True              ; Shading on/off: True/False     
    ; res@cnLevelSpacingF      = 10 ; contour spacing
    ; res@cnMinLevelValF       = -100                 ; min level
    ; res@cnMaxLevelValF       = 100                ; max level

     ;res@mpGridLatSpacingF    = 60.                ; change latitude  line spacing
     ;res@mpGridLonSpacingF    = 90.                ; change longitude line spacing

     ;res@gsnPolarLabelFontHeightF = 0.01
     res@lbLabelBarOn             = False
     res@gsnLeftStringFontHeightF = 0.029
     res@gsnStringFont = 22
     res@gsnLeftStringFontHeightF=0.05
     res@vpHeightF=1
     res@vpWidthF=1
    res@gsnPolarLabelFontHeightF = .05
    res@gsnLeftString = " "
     res@gsnPolar             = "NH"               ; specify the hemisphere
     res@lbBoxEndCapStyle      = "TriangleBothEnds"

yr=(/1991,1994,1995,1999,2002,2003,2006,2009,2013,2014,2015/);指定年份
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
do i=0,10,1
f=addfile("/data/hgt/hgt."+yr(i)+".nc","r")

TIME   = f->time
YYYY   = cd_calendar(TIME,-1)/100                 ; entire file
MON    = cd_calendar(TIME,-1)-YYYY*100
;f1    =addfile("/hgt/hgt.2012.nc","r")
     hc=f->hgt(0:35,{level},{0:90},:)  ;1月份
     ;hc=dim_rmvmean_n_Wrap(h,0)
  printVarSummary(hc)
   ;;;;;;创建工作台
    wks = gsn_open_wks("png", "/hgt"+yr(i)+"")
     ; gsn_define_colormap(wks,"BlueWhiteOrangeRed")
     gsn_define_colormap(wks,"BlueDarkRed18")               ; choose colormap

     plot                                    = new(16,graphic)   
     r=new((/16,16,37,144/),float)
     hh=new((/16,37,144/),float)
     ;h1=new((/1,37,144/),float)
     h=new((/16,16,37,144/),float)

     do j=0,30,2
             m=j/2
             hh(m,:,:)=hc(j,:,:)
             r(m,:,:,:)=hh
             ;h(m,:,:,:)=dim_standardize_n_Wrap(r(m,:,:,:),opt,1)
             printVarSummary(r)

     end do
;res@gsnCenterString=""+name+" "+level+" hPa"
     do m=0,15,1

       day=m*2+1
       res@gsnLeftString=""+day+" JAN "+yr(i)+""
      plot(m)= gsn_csm_contour_map_polar(wks,r(m,m,:,:),res) ; create plot object
     end do
     resp = True
     resp@gsnPanelRowSpec          = True
     resp@gsnPanelLabelBar         = True
     resp@pmLabelBarOrthogonalPosF = -0.007
     resp@pmLabelBarWidthF         = 0.7   ;长度
     resp@pmLabelBarHeightF        = 0.05  ;宽度
     resp@lbLabelFontHeightF       = 0.012
     gsn_panel(wks,plot,(/4,4,4,4/), resp)
end do
end