萌新求助

做了一个关于臭氧的多元线性回归，最后输出的变量b@Yest是系数，求助怎么把这个改成合成的臭氧序列，多谢
附上程序
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
; ==============================================================
; Open the file: Read only the user specified period
; ==============================================================
lat0s = -10
lat0e = 10
lon0 = 0
lon1 = 360
lev0s = 10               ;hPa                ;hPa
lev0e = 300
LEV   = 20 ;

  f0      = addfile("MLS_O3_monthly_2004-2015.nc","r")
  f1      = addfile ("u.prod.assim.instM_3d_asm_Cp.1979-2015.SUB.nc", "r")
con0  = f0->O3(12:131,{lev0s:lev0e},{lat0s:lat0e},:)
con1  = f1->u(312:431,{lev0s:lev0e},{lat0s:lat0e},:)

mon0 = con0&time
lev0=con0&levelist

;+++++++++++++++++++
;calulate MonAnomTLLL
;+++++++++++++++++++
o3Clm  = clmMonTLLL(con0)
o3Anom =calcMonAnomTLLL(con0,o3Clm)

;*********************
;Interpolates U pressure levels
;*********************
u = int2p_n (con1&levelist,con1,lev0,1,1)
u!0 = "month"
u!1 = "level"
u!2 = "lat1"
u!3 = "lon1"
u&month = con1&time
u&level = lev0
u&lat1  = con1&latitude
u&lon1  = con1&longitude
nlev0 = dimsizes(lev0)
nmon  = dimsizes(u&month)

mon1 = con1&time
lat1 = con1&latitude
lev1 = con1&levelist
nlev1 = dimsizes(lev1)       ; longitudes to span -180 to 177.5: facilitate coordinate subscripting

;+++++++++++++++++++
;calulate MonAnomTLLL
;+++++++++++++++++++
uClm  = clmMonTLLL(u)
u_Anom =calcMonAnomTLLL(u,uClm)

; ==============================================================
;EOF_QBO
; ==============================================================
  neof   = 3        ; number of EOFs
  optEOF = True      
  optEOF@jopt = 0   ; This is the default; most commonly used; no need to specify.
;;optEOF@jopt = 1   ; **only** if the correlation EOF is desired
  optETS = False
; =================================================================
; create weights:  sqrt(cos(lat))   [or sqrt(gw) ]
; =================================================================
  rad    = 4.*atan(1.)/180.
  clat   = con1&latitude           
  clat   = sqrt( cos(rad*clat) )                 ; gw for gaussian grid

; =================================================================
; Reorder (lat,lon,time) the *weighted* input data
; Access the area of interest via coordinate subscripting
; =================================================================
  xw     = u_Anom(level|:,lat1|:,lon1|:,month|:)
  x      = u_Anom(month|:,level|:,lat1|:,lon1|:)

  eof      = eofunc_Wrap(xw, neof, optEOF)     
;  eof_ts   = eofunc_ts_Wrap (xw, eof, optETS)
  eof_ts   = new((/neof,dimsizes(con1&time),nlev0/),float)

do ip=0,nlev0-1
  eof_ts(:,:,ip)   = eofunc_ts_Wrap (xw(ip,:,:,:), eof(:,ip,:,:), optETS)
end do
eof_ts!2 = "levs"
eof_ts&levs = lev0

  printVarSummary( eof )                         ; examine EOF variables
  printVarSummary( eof_ts )
; =================================================================
; ENSO
; =================================================================

  MEI = asciiread("M ENSO Indx1979_2015.txt",-1,"float")         ;!time
pMEI = MEI(312:431)
printVarSummary(pMEI)
; =================================================================
; SOLAR CYCLE
; =================================================================
solar = asciiread("solar_cycleSN_m_tot_V2.0_1979_2015.txt",-1,"float")
psolar=solar(312:431)      ;!time
printVarSummary(psolar)

; =================================================================
; trend
; =================================================================
o3=dim_avg_n_Wrap(dim_avg_n_Wrap(con0(:,:,:,:),2),2)
trend= o3
do ip=0,nlev0-1
  rc=regline(mon0,o3(:,ip))
  trend(:,ip)=rc*(mon0-rc@xave) + rc@yave          ; use solid line
end do


; =================================================================
; regres
; =================================================================
X=new((/120,5,nlev0/),float)

  do ip=0,nlev0-1
   X(:,0,ip)=trend(:,ip)
   X(:,1,ip)=psolar(:)
   X(:,2,ip)=eof_ts(0,:,ip)
   X(:,3,ip)=eof_ts(1,:,ip)
   X(:,4,ip)=pMEI(:)
  end do
po3Anom=dim_avg_n_Wrap(dim_avg_n_Wrap(o3Anom(:,:,:,:),2),2)

Y=po3Anom
opt = True
opt@print_anova = True

printVarSummary(Y)

sig=new((/120,5,nlev0/),float)
b=new((/6,nlev0/),float)

do ip=0,nlev0-1
b(:,ip)=reg_multlin_stats(Y(:,ip),X(:,:,ip),opt) ; partial regression coef
sig(:,0,ip)=b(2,ip)*eof_ts(0,:,ip)+b(0,ip)
sig(:,1,ip)=b(3,ip)*eof_ts(1,:,ip)+b(0,ip)
sig(:,2,ip)=b(4,ip)*pMEI(:)*(1e8)+b(0,ip)
end do
sig!0="month"
sig!2 = "level"
sig&month = con1&time
sig&level = lev0

plot_QBO1=sig(:,0,{LEV})
plot_QBO2=sig(:,1,{LEV})
plot_MEI =sig(:,2,{LEV})

;************************************************
data     = new((/5,nmon/),float)
data(0,:) = plot_QBO1
data(1,:) = plot_QBO2
data(2,:) = plot_MEI
data(3,:) =  po3Anom(:,{20})
data(4,:) = b@Yest


;plot_trend=dim_avg_n_Wrap(trend,0)
; data_o3      = new((/3,nlev0/),float)
; data_o3(0,:) = plot_trend
; data_o3(1,:) =  0
;************************************************
wks   = gsn_open_wks ("png","1")                  ; open ps file
; wks   = gsn_open_wks  ("png","2")                  
;************************************************
res                   = True  
res@vpHeightF   = 0.4               ; change aspect ratio of plot
res@vpWidthF           = 0.7                     ; plot mods desired
  res@tiYAxisString    = " "    ; y-axis label     

res@tmXBValues=(/13,   37,  61,  85,  109,121/)
res@tmXBLabels=(/"2005","2007","2008","2010","2012","2013"/)
; res@tiMainString      = "o3 trend"+lat0s+"_"+lat0e+"_MERRA"             ; add title
;res@trYReverse        = True                       ; reverse Y-axis
; add a legend
;res@pmLegendDisplayMode    = "Always"              ; turn on legend(璇存槑绾?
;res@pmLegendSide           = "Top"                 ; Change location of
;res@pmLegendParallelPosF   = 0.90                   ; move units right
;res@pmLegendOrthogonalPosF = -0.8                  ; more neg = down
res@pmLegendWidthF         = 0.12                  ; Change width and
res@pmLegendHeightF        = 0.25                  ; height of legend.
res@lgLabelFontHeightF     = .02                   ; change font height
;res@lgPerimOn              = False                 ; no box around
; labels for the legend
; res@xyExplicitLegendLabels = (/lat0s+"-"+lat0e," "/)
;************************************************
; Change to log scaling for Y Axis
;************************************************
;res@xyYStyle          = " Linear"
;res@tmYLMode          = "Explicit"                ; explicit labels
;res@tmYLValues        = (/300,200,100,50,30,10,5/)
;res@tmYLLabels        = ""+res@tmYLValues         ; make strings
;res@xyLineColors      =(/"blue","darkgreen","black"/)
;res@xyLineThicknesses =(/3.0,3.0,3.0,3.0,3.0/)
;res@xyComputeYMin = True
res@gsnFrame          = False                      ; don't advance frame yet
res@gsnDraw          = False

plot=new(5,graphic)
do i=0,4
plot(0)  = gsn_csm_xy (wks,mon0,data(0,:),res)           ; data(0,:) = plot_QBO1
plot(1)  = gsn_csm_xy (wks,mon0,data(1,:),res)            ; data(1,:) = plot_QBO2
plot(2)  = gsn_csm_xy (wks,mon0,data(2,:),res)            ;data(2,:) = plot_MEI
plot(3)  = gsn_csm_xy (wks,mon0,data(3,:),res)             ; data(3,:) =  po3Anom(:,{20})
plot(4)  = gsn_csm_xy (wks,mon0,data(4,:),res)             ; data(4,:) = b@Yest
end do
;---Arrays to hold text annotation ids
  txid_tr = new(5,graphic)
  amid_tr = new(5,graphic)

  txres                       = True
  txres@txPerimOn             = True
  txres@txFontHeightF         = 0.04

;---Top right string
  amres_tr                  = True
  amres_tr@amParallelPosF   =  0.5    ; This is the right edge of the plot.
  amres_tr@amOrthogonalPosF = -0.5    ; This is the top edge of the plot.
  amres_tr@amJust           = "TopRight"

Ti=(/"QBO1","QBO2","MEI","po3Anom","b@Yest"/)
do i=0,4
  tr_label = Ti(i)

    txres@txBackgroundFillColor = "Goldenrod"
    txid_tr(i) = gsn_create_text(wks, tr_label, txres)

    txres@txBackgroundFillColor = "LawnGreen"      ; just for something different

;---Attach text strings to plot
    amid_tr(i) = gsn_add_annotation(plot(i), txid_tr(i), amres_tr)
end do

pres=True
;pres@gsnPanelCenter=False
;pres@gsnPanelRowSpec=True

gsn_panel(wks,plot,(/5,1/),pres)
pres@lbLabelBarOn = True
pres@gsnPanelLabelBar = True

end