;; Empty NCL template
;; run code in recipes, use it carefully

;;DIAG_NORCPM; RUNTHESECODES: print("    no ncl code run here")

begin
    obsfn = "amip_tos_NINO3_ts.nc"
    obsvn = "tos"
    modelfn = "sst_superensmean_NINO3_ts.nc"
    modelvn = "sst"
    figfn = "superensmean_sst_NINO3_ts"
    title = "super ensmean NINO3"
    rmann = True
    detr  = False
    ;; the methods to present members
    ;;  values are "line", "line25"
    members = "shade25"

    f = addfile(obsfn,"r")
    obsts = f->$obsvn$
    f = addfile(modelfn,"r")
    modelts = f->$modelvn$

    obsta = cd_calendar(obsts&time,0)
    nt = dimsizes(obsts)
    mdims = dimsizes(modelts)


    if(rmann)then  ;; remove annual cycle
        if (dimsizes(mdims).gt.1)then
            do i = 0, mdims(0)-1
                modelts(i,:) = rmAnnCycle1D(modelts(i,:))
            end do
        else
            modelts = rmAnnCycle1D(modelts)
        end if
        obsts = rmAnnCycle1D(obsts)
    end if
    if(detr)then  ;; remove trend
        modelts = dtrend(modelts,False)
        obsts = dtrend(obsts,False)
    end if

    if (dimsizes(mdims).gt.1)then
        ;; mean correlation is meanless
        ;;r = 0.
        ;;do i = 0,mdims(0)-1
        ;;    r = r+ escorc(tofloat(obsts),tofloat(modelts(i,:)))/mdims(0) 
        ;;end do
        ;;rstring = "mean R = "+sprintf("%3.3f",r) 
        r = escorc(tofloat(obsts),tofloat(dim_avg_n(modelts,0))) ;; need be float, do not know why
        rstring = "ensmean R = "+sprintf("%3.3f",r) 
    else
        r = escorc(tofloat(obsts),tofloat(modelts)) ;; need be float, do not know why
        rstring = "R = "+sprintf("%3.3f",r) 
    end if

    xbvalues = ispan(0,nt-1,60)
    xbminorvalues = ispan(0,nt,12)

    wks = gsn_open_wks("ps",figfn)
    plotres = True
        plotres@vpHeightF = 0.2
        plotres@gsnLeftString = title
        plotres@gsnRightString = rstring
        plotres@xyLineColors = (/"black","darkgreen"/)
        plotres@xyLineThicknesses = (/3,2/)
        plotres@gsnYRefLine = 0.
        plotres@trYMaxF = 6.
        plotres@trYMinF = -4.
        plotres@tmXBMode = "Explicit"
        ;plotres@tmLabelAutoStride = True
        plotres@tmXBValues = xbvalues
        plotres@tmXBLabels = obsta(xbvalues,0)
        plotres@tmXBLabelJust = "BottomLeft"
        plotres@xyDashPatterns = (/0,0/)
        plotres@tmXTOn = False
        plotres@tmXBMinorValues = xbminorvalues

        plotres@pmLegendDisplayMode = "Always"
        plotres@pmLegendSide = "Top"
        plotres@pmLegendParallelPosF   = .7                   ; move units right
        plotres@pmLegendOrthogonalPosF = -0.40               ; move units down
        plotres@pmLegendWidthF         = 0.05                ; Change width and
        plotres@pmLegendHeightF        = 0.05                ; height of legend.
        plotres@lgPerimOn              = False               ; turn off box around
        plotres@lgLabelFontHeightF     = .015                ; label font height
        plotres@xyExplicitLegendLabels = (/"Model","Obs"/)
    
        plotres@gsnFrame = False
        plotres@gsnDraw  = False

        lineres = True
            lineres@gsLineColor = "grey44"
            lineres@tfPolyDrawOrder = "PreDraw"

    if (dimsizes(mdims).gt.1)then
        plot = gsn_csm_xy(wks,ispan(0,nt-1,1),(/dim_avg_n(modelts,0),obsts/),plotres)
        if (members .eq. "lines")then
            do i = 0, mdims(0)-1
                gsn_polyline(wks,plot,ispan(0,nt-1,1),modelts(i,:),lineres)
            end do
        end if
        if (members .eq. "line25" .or. members .eq. "shade25")then
            upper = new(mdims(1),typeof(modelts))
            lower = new(mdims(1),typeof(modelts))
            xi = fspan(1,100,mdims(0))
            xo = (/25.,75./)
            do t = 0, mdims(1)-1
                fi = modelts(:,t)
                qsort(fi)
                res = linint1(xi,fi,False,xo,0)
                lower(t) = res(0)
                upper(t) = res(1)
            end do
            if(members.eq. "line25" )then
                gsn_polyline(wks,plot,ispan(0,nt-1,1),upper,lineres)
                gsn_polyline(wks,plot,ispan(0,nt-1,1),lower,lineres)
            end if
            if(members.eq. "shade25" )then
                do t = 0, mdims(1)-1
                    gsn_polyline(wks,plot,(/t,t/),(/upper(t),lower(t)/),lineres)
                end do
            end if
        end if
    else
        plot = gsn_csm_xy(wks,ispan(0,nt-1,1),(/modelts,obsts/),plotres)
    end if
    draw(plot)
    frame(wks)
end