;;DIAG_NORCPM; PLOTRES: ;; PLOTRES section, empty
;;DIAG_NORCPM; PYEARB: 0    ;; filter year begin
;;DIAG_NORCPM; PYEARE: 9999 ;; filter year end
;;DIAG_NORCPM; YRSANOMALY: False  ;; subtract all year clm

;; plot seasonal longitude-pressure vector figure

;; ref https://www.ncl.ucar.edu/Applications/Scripts/vector_5.ncl
;; input files, should be (time,plev,1,lon), time need be multiyear monthly data

ifileu = "U_ensmean.nc"
ifilew = "OMEGA_ensmean.nc"
vnu = "U"
vnw = "OMEGA"
figfn = "walker_SON_P1"
title = "Walker circ. SON P1 1980-2000"
season = "SON"
yb = 1980
ye = 2000
yrsano = False  ;; subtract all year clm

;; for vector length
if(yrsano)then
    vcrefm = 5.
    vcrefl = 0.045
else
    vcrefm = 50.
    vcrefl = 0.045
end if

;; read data
f = addfile(ifileu,"r")
u = f->$vnu$(:,:,0,:)
time = f->time
f = addfile(ifilew,"r")
w = f->$vnw$(:,:,0,:)

if(yrsano)then
    useaall = month_to_season(u,season)
    wseaall = month_to_season(w,season)
    useaclm = dim_avg_n_Wrap(useaall,0)
    wseaclm = dim_avg_n_Wrap(wseaall,0)
end if

time = time -10 ;; avoid Jan data at 01Feb
t = cd_calendar(time,-1) ;; YYYYMM
needt = ind(t.ge.(yb*100) .and. t.le.(ye*100))
u := u(needt,:,:)
w := w(needt,:,:)

usea = month_to_season(u,season)
wsea = month_to_season(w,season)
useamean = dim_avg_n_Wrap(usea,0)
wseamean = dim_avg_n_Wrap(wsea,0)

if(yrsano)then
    useamean = useamean - useaclm
    wseamean = wseamean - wseaclm
end if


;;scale = abs(avg(useamean)/avg(wseamean))
scale = 400. ;; tuned
wseamean100 = wseamean*scale

res = True
    res@tiMainString = title
    res@vcRefMagnitudeF = vcrefm             ; define vector ref mag
    res@vcRefLengthF    = vcrefl             ; define length of vec ref
    res@vcGlyphStyle    = "CurlyVector"      ; turn on curly vectors
    res@vcMinDistanceF  = 0.02               ; thin out vectors
    res@vcMapDirection  = False

    res@vpHeightF = 0.2
    res@vpWidthF  = 0.7
    res@vcRefAnnoOn = False
    res@cnFillPalette = "BlAqGrWh2YeOrReVi22"
    res@cnLevelSelectionMode = "ManualLevels"
    if(yrsano)then
        res@cnMaxLevelValF  =  0.01
        res@cnMinLevelValF  = -0.01
        res@cnLevelSpacingF =  0.001
    else
        res@cnMaxLevelValF  =  0.1
        res@cnMinLevelValF  = -0.1
        res@cnLevelSpacingF =  0.01
    end if
    res@cnLinesOn = False
    res@cnFillOn = True

;; PLOTRES section, empty

wks = gsn_open_wks("ps",figfn)
plot = gsn_csm_pres_hgt_vector(wks,wseamean,useamean,-wseamean100,res)
                                      ;; omega, negative means up