;;diag_norcpm;;amip_tos_NINO3_ts.nc: 
;;diag_norcpm;;amip_tos_Atl3_ts.nc: 
;;diag_norcpm;;tos: 
;;diag_norcpm;;sst_ens_NINO3_ts.nc:
;;diag_norcpm;;sst_ens_Atl3_ts.nc:
;;diag_norcpm;;sst:
;;diag_norcpm;;True: False
;;diag_norcpm;;DETREND: False
;;diag_norcpm;;0.5:  0.5
;;diag_norcpm;;-0.5: -0.5
;;diag_norcpm;;NUG_FF_5m NINO3 Atl3 lag cor. 1980-2000: NUG_FF_5m NINO3 Atl3 lag cor. 1980-2000
;;diag_norcpm;;lagcor_NINO3_Atl3:  fig
begin
    obs1fn = "amip_tos_NINO3_ts.nc"
    obs2fn = "amip_tos_Atl3_ts.nc"
    obsvn = "tos"
    model1fn = "sst_ens_NINO3_ts.nc"
    model2fn = "sst_ens_Atl3_ts.nc"
    modelvn = "sst"
    tag1 = "NINO3"
    tag2 = "Atl3"
    figfn = "lagcor_NINO3_Atl3"
    title = "NUG_FF_5m NINO3 Atl3 lag cor. 1980-2000"
    rmann = True
    do_detrend = DETREND
    monthrange = 18  ;; 

    f = addfile(obs1fn,"r")
    obs1ts = f->$obsvn$
    f = addfile(obs2fn,"r")
    obs2ts = f->$obsvn$
    f = addfile(model1fn,"r")
    model1ts = f->$modelvn$
    f = addfile(model2fn,"r")
    model2ts = f->$modelvn$

    ;; assume 2 obs ts are in same time 
    obsta = cd_calendar(obs1ts&time,0)
    ;; assume 2 model ts are in same time 
    modelta = cd_calendar(model1ts&time,0)
    
    mdims = dimsizes(model1ts)

    if(rmann)then  ;; remove annual cycle
        if (dimsizes(mdims).gt.1)then
            do i = 0, mdims(0)-1
                model1ts(i,:) = rmAnnCycle1D(model1ts(i,:))
                model2ts(i,:) = rmAnnCycle1D(model2ts(i,:))
            end do
        else
            model1ts = rmAnnCycle1D(model1ts)
            model2ts = rmAnnCycle1D(model2ts)
        end if
        obs1ts = rmAnnCycle1D(obs1ts)
        obs2ts = rmAnnCycle1D(obs2ts)
    end if
    if(do_detrend)then  ;; remove annual cycle
        if (dimsizes(mdims).gt.1)then
            do i = 0, mdims(0)-1
                model1ts(i,:) = dtrend_n(model1ts(i,:),False,0)
                model2ts(i,:) = dtrend_n(model2ts(i,:),False,0)
            end do
        else
            model1ts = dtrend_n(model1ts,False,0)
            model2ts = dtrend_n(model2ts,False,0)
        end if
        obs1ts = dtrend_n(obs1ts,False,0)
        obs2ts = dtrend_n(obs2ts,False,0)
    end if

    if (dimsizes(mdims).gt.1)then
        modelens1ts = dim_avg_n_Wrap(model1ts,0)
        modelens2ts = dim_avg_n_Wrap(model2ts,0)
    else
        modelens1ts = model1ts
        modelens2ts = model2ts
    end if

    r_obs   = new(2*monthrange+1,float)
    r_obs(monthrange:)     = esccr(tofloat(obs2ts),  tofloat(obs1ts),  monthrange)   ;; left part, 2 lead 1
    r_obs(:monthrange)     = esccr(tofloat(obs1ts),  tofloat(obs2ts),  monthrange)   ;; right part,1 lead 2
    r_obs(:monthrange)     = r_obs(:monthrange:-1)

    r_modelens   = new(2*monthrange+1,float)
    r_modelens(monthrange:)     = esccr(tofloat(modelens2ts),  tofloat(modelens1ts),  monthrange)   ;; left part, 2 lead 1
    r_modelens(:monthrange)     = esccr(tofloat(modelens1ts),  tofloat(modelens2ts),  monthrange)   ;; right part,1 lead 2
    r_modelens(:monthrange)     = r_modelens(:monthrange:-1)

    if (dimsizes(mdims).gt.1)then
        r_model = new((/mdims(0),2*monthrange+1/),float)
        do i = 0,mdims(0)-1
            r_model(i,monthrange:)   = esccr(tofloat(model2ts(i,:)),  tofloat(model1ts(i,:)),  monthrange)
            r_model(i,:monthrange)   = esccr(tofloat(model1ts(i,:)),  tofloat(model2ts(i,:)),  monthrange)
            r_model(i,:monthrange)   = r_model(i,:monthrange:-1)
        end do
    else
        r_model = new(2*monthrange+1,float)
        r_model(monthrange:)   = esccr(tofloat(model2ts),  tofloat(model1ts),  monthrange)
        r_model(:monthrange)   = esccr(tofloat(model1ts),  tofloat(model2ts),  monthrange)
        r_model(:monthrange)   = r_model(:monthrange:-1)
    end if

    xbvalues = ispan(-monthrange,monthrange,6)
    xbminorvalues = ispan(-monthrange,monthrange,1)

    wks = gsn_open_wks("ps",figfn)
    plotres = True
        plotres@vpHeightF = 0.2
        plotres@gsnLeftString = title
        plotres@xyLineColors = (/"darkgreen","black"/)
        plotres@xyLineThicknesses = (/3,3/)
        plotres@gsnYRefLine = 0.
        plotres@gsnXRefLine = 0.
        plotres@trYMaxF = 0.5
        plotres@trYMinF = -0.5
        ;;plotres@tmXBMode = "Explicit"
        ;plotres@tmLabelAutoStride = True
        ;;plotres@tmXBValues = xbvalues
        ;;plotres@tmXBLabels = obsta(xbvalues,0)
        plotres@tmXBLabelJust = "CenterCenter"
        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 = (/"Obs","Model(EnsMean)"/)

    lineres = True
        lineres@gsLineColor = "black"

    textres = True
        textres@txFontHeightF = 0.02

    ;; PLOTRES

    plotres@gsnFrame = False
    if (dimsizes(mdims).gt.1)then
        ;;plot = gsn_csm_xy(wks,xbminorvalues,(/r_obs,dim_avg_n(r_model,0)/),  plotres)
        plot = gsn_csm_xy(wks,xbminorvalues,(/r_obs,r_modelens/),  plotres)
        do i = 0,mdims(0)-1
            gsn_polyline(wks,plot,xbminorvalues,r_model(i,:),lineres)
        end do
    else
        plot = gsn_csm_xy(wks,xbminorvalues,(/r_obs,r_model/),  plotres)
    end if

    gsn_text(wks,plot,tag2+ " lead "+tag1, 10,-0.42,textres)
    gsn_text(wks,plot,tag1+ " lead "+tag2,-10,-0.42,textres)

    frame(wks)
end