;; plot time correlation map between 2 file variables
;; variable should be (time,y,x)

;;DIAG_NORCPM; FIGFILENAME: fig
;;DIAG_NORCPM; FIGTITLE: TITLE
;;DIAG_NORCPM; PLOTRES:  ;; plot resource
;;DIAG_NORCPM; INPUTFILE: 
;;DIAG_NORCPM; INPUTVAR: 
;;DIAG_NORCPM; RMANNUAL: False
;;DIAG_NORCPM; YEARB: 0
;;DIAG_NORCPM; YEARE: 0
;;DIAG_NORCPM; DETREND: False
;;DIAG_NORCPM; SEASONDJF: False
;;DIAG_NORCPM; ISOCNGRID: False
;;DIAG_NORCPM; FILLNA: True
;;DIAG_NORCPM; VARPREPROC: ;; no var preproc

begin
    title = "AMIP tos variance P2 (annual removed)"
    figfn = "var_amip_tos_p2"
    yearb = 2000
    yeare = 2020

    ;; var(t,y,x)
    f = addfile("amip_tos.nc","r")
    var = f->tos

    ;; no var preproc

    if (True)then
        guess     = 1                ; use zonal means
        is_cyclic = True             ; cyclic [global]
        nscan     = 1500             ; usually much less than this
        eps       = 1.e-2            ; variable dependent
        relc      = 0.6              ; relaxation coefficient
        opt       = 0                ; not used
        poisson_grid_fill( var, is_cyclic, guess, nscan, eps, relc, opt)
    end if

    if (yearb.ne.0 .and. yeare.ne.0)then ;; filter years
        time = var&time
        lt1 = var&time
        t1 = toint(lt1)
        copy_VarAtts(lt1,t1)
        years1 = cd_calendar(t1,-1)/100
        lyears1 = years1.ge.yearb .and. years1.le.yeare
        var := var(ind(lyears1),:,:)
    end if

    if (True)then ;; remove annual cycle
        var = rmMonAnnCycTLL(var)
    end if

    if (False)then ;; remove trend
        var = dtrend_n(var,False,0)
    end if

    var_var = dim_variance_n_Wrap(var,0)

    ;; do plot
    res = True
        res@gsnDraw     = False
        res@gsnFrame    = False
        res@cnFillOn    = True
        res@cnFillMode  = "CellFill"
        res@cnLinesOn   = False
        res@lbOrientation     = "Horizontal"
        res@gsnAddCyclic = True
        res@gsnLeftString = title
        res@cnLevelSelectionMode = "ExplicitLevels"
        res@cnLevels = ispan(0,19,1)/10.
        res@mpFillDrawOrder = "PostDraw"
        res@cnFillPalette = "WhiteBlueGreenYellowRed"
        res@mpCenterLonF = 210.

    ;; not use here, yet?
        res2 = True ;; for shade on filled contour
            res2@cnLinesOn = False
            res2@cnLineLabelsOn = False
            res2@gsnDraw = False
            res2@gsnFrame = False
            res2@gsnLeftString = ""
            res2@gsnRightString = ""
            res2@cnInfoLabelOn     = False
            ;; shading only fill between contours
            res2@cnLevelSelectionMode = "ExplicitLevels"
            ;;res2@cnLevels = (/psig,max(p)/)
            
        resShade = True
            resShade@gsnShadeFillType = "pattern"
            resShade@gsnShadeHigh   = 17
            resShade@gsnShadeLow = 17

    ;; assign coordinates
    if (False)then
        gridf = addfile("/nird/projects/NS9039K/shared/pgchiu/diag_norcpm/grid_norcpm1_ocn.nc","r")
        plat = gridf->plat
        plon = gridf->plon
        var_var@lon2d = plon
        var_var@lat2d = plat
    end if

    ;; plot resource

    wks = gsn_open_wks("ps",figfn)

    plot_fill = gsn_csm_contour_map(wks,var_var,res)

    ;;plot_shade = gsn_csm_contour(wks,p,res2)
    ;;plot_shade = gsn_contour_shade(plot_shade,psig,max(p),resShade)
    ;;overlay(plot_fill,plot_shade)
    draw(plot_fill)

    frame(wks)

end