;; 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; INFILEMAP: 
;;DIAG_NORCPM; INVARMAP: 
;;DIAG_NORCPM; INFILETS: 
;;DIAG_NORCPM; INVARTS: 
;;DIAG_NORCPM; RMANNUAL: False
;;DIAG_NORCPM; DETREND: False
;;DIAG_NORCPM; SEASONDJF: False
;;DIAG_NORCPM; ADDBOX: False
;;DIAG_NORCPM; BOXRES: ;; box plot resource
;;DIAG_NORCPM; ISOCNGRID: False
;;DIAG_NORCPM; VARPREPROC: ;; no var preproc

begin
    title = "NUG_Ano ens. cor. U850-Atl3 SST P1 1980-2000"
    figfn = "cor_ens_U850-Atl3SST_P1"
    addbox = True  ;; plot latlon box on map
    ;; need additional latb, late, lonb, lone

    ;; var(t,y,x)
    f1 = addfile("NUG_Ano_U850_ensmean_1980-2020.nc","r")
    var1 = f1->U850
    f2 = addfile("../01_sst_regional-cor/sst_ens_Atl3_ts.nc","r")
    var2 = f2->sst

    ;; ensemble mean
    var2 := dim_avg_n_Wrap(var2,0)
    ;; filter year
    yearb = 1980
    yeare = 2000
    lt1 = f1->time
    t1 = toint(lt1)
    copy_VarAtts(lt1,t1)
    years1 = cd_calendar(t1,-1)/100
    lyears1 = years1.ge.yearb .and. years1.le.yeare
    var1 := var1(ind(lyears1),:,:)
    lt2 = f2->time
    t2 = toint(lt2)
    copy_VarAtts(lt2,t2)
    years2 = cd_calendar(t2,-1)/100
    lyears2 = years2.ge.yearb .and. years2.le.yeare
    var2 := var2(ind(lyears2))
    

    ;; strange, sometime ncl read NaN as integer
    if (typeof(var1@_FillValue).eq."integer")then
        var1@_FillValue = -1e20
    end if
    if (typeof(var2@_FillValue).eq."integer")then
        var2@_FillValue = -1e20
    end if
    ;; escorc cannot handle nan as fillvalue
    if (isnan_ieee(var1@_FillValue))then
        var1@_FillValue = -1e20
        replace_ieeenan (var1, var1@_FillValue, 0)
    end if
    if (isnan_ieee(var2@_FillValue))then
        var2@_FillValue = -1e20
        replace_ieeenan (var2, var2@_FillValue, 0)
    end if

    dims1 = dimsizes(var1)
    dims2 = dimsizes(var2)

    if (False)then
        ;; assume DJF if used as seasonal mean
        ;; 1st year is JF and last year is D only
        ;; Use year in coordinate to align data
        ;; And discard last year
        t1a = cd_calendar(var1&time,0)
        t2a = cd_calendar(var2&time,0)

        yb = max((/min(t1a(:,0)),min(t2a(:,0))/))
        ye = min((/max(t1a(:,0)),max(t2a(:,0))/)) -1
        it1 = ind(t1a(:,0).ge.yb .and. t1a(:,0).le.ye)
        it2 = ind(t2a(:,0).ge.yb .and. t2a(:,0).le.ye)
        if (any(t1a(it1,0).ne.t2a(it2,0)))then
            print("ERROR: Year not align")
            exit 
        end if
        if(dimsizes(dims1).eq.1)then
            var1 := var1(it1)
        else
            var1 := var1(it1,:,:)
        end if
        if(dimsizes(dims2).eq.1)then
            var2 := var2(it2)
        else
            var2 := var2(it2)
        end if

    end if

    if (True)then ;; remove annual cycle
        var1 = rmMonAnnCycTLL(var1)
        if (dimsizes(dims2).eq.1)then
            tvar2 = new((/dims2,1,1/),typeof(var2))
            ;; avoid warning
            tvar2!1 = "dim1"
            tvar2!2 = "dim2"
            tvar2(:,0,0) = var2
            tvar2 = rmMonAnnCycTLL(tvar2)
            var2(:) = tvar2(:,0,0)
        else
            var2 = rmMonAnnCycTLL(var2)
        end if

        ;;title = title+" (ACC)"
    end if
    if (False)then ;; remove trend
        var1 = dtrend_n(var1,False,0)
        var2 = dtrend_n(var2,False,0)
    end if

    if(dimsizes(dims2).eq.1)then  ;; fewer dimensions must be seconded
        r = escorc_n(var1,var2,0,0) ;; 
    else
        ;; check dims, should be identical
        if (any(dims1.ne.dims2))then
            print("dims not match")
            print("  dims MAP: "+str_join(dims1,","))
            print("  dims TS:  "+str_join(dims2,","))
            print("Maybe the grid file is wrong?")
            print("exit...")
            exit
        end if
        ;; array with missing values must be first
        if (any(ismissing(var1)))then
            r = escorc_n(var1,var2,0,0) ;; 
        else
            r = escorc_n(var2,var1,0,0) ;; 
        end if
    end if

    ;; statistic tests
    ;; from https://www.ncl.ucar.edu/Document/Functions/Built-in/escorc_n.shtml

    n = dims1(0)
    df = n-2                    ;; degree for freedom

    ;; Fischer z-transformation
    z = 0.5*log((1+r)/(1-r))    ;; z-statistic
    se = 1.0/sqrt(n-3)          ;; standard error of z-statistic
    zlow = z-1.96*se            ;; low and high z vzlues
    zhi  = z+1.96*se            ;; 95%: 1.96, 99%: 2.58
    ;; inverse z-transform; return to r space (-1 to +1)
    rlow = (exp(2*zlow)-1)/(exp(2*zlow)+1)
    rhi  = (exp(2*zhi )-1)/(exp(2*zhi )+1)

    ;; student test
    t    = r*sqrt((n-2)/(1-r^2))
    p    = student_t(t, df)
    psig = 0.05                       ;; test significance level
                                      ;; 

    ;; 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@gsnRightString = "p siglvl=0.05"
        ;res@cnLevels = (/-.9,-.8,-.7,-.5,-.2,.2,.5,.7,.8,.9/)
        ;res@cnLevelSelectionMode = "ExplicitLevels"


    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
        res2@gsnAddCyclic = True
        ;; 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/datalake/NS9039K/users/pgchiu/diag_norcpm/grid_tnx1v4_20170622.nc","r")
        plat = gridf->plat
        plon = gridf->plon
        r@lon2d = plon
        r@lat2d = plat
        p@lon2d = plon
        p@lat2d = plat
    else
        copy_VarCoords(var1(0,:,:),r)
        copy_VarCoords(var1(0,:,:),p)
    end if

    res@mpCenterLonF = 210.
    res@gsnAddCyclic = True
    res@cnFillColors = (/(/130,33,239/), (/0,0,150/), (/0,0,204/), (/63,104,224/), (/30,142,255/), (/0,191,255/), (/160,209,255/), (/209,244,255/), (/255,255,255/), (/255,255,255/), (/255,255,255/), (/255,255,255/), (/255,255,198/), (/255,224,51/), (/255,170,0/), (/255,109,0/), (/255,0,0/), (/198,0,0/), (/160,35,35/), (/255,104,178/)/)/255.
    res@cnLevelSelectionMode = "ExplicitLevels"
    res@cnLevels = ispan(-9,9,1)/10.
            
    

    wks = gsn_open_wks("ps",figfn)

    plot_fill = gsn_csm_contour_map(wks,r,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)

    if (addbox)then 
        latbe = (/-3.,3./)
        lonbe = (/-20.,0./)
        xb = min(lonbe)
        xe = max(lonbe)
        yb = min(latbe)
        ye = max(latbe)

        boxres = True
        ;; box plot resource
        
        ;; counter clockwise
        x = (/xb,xb,xe,xe,xb/)
        y = (/yb,ye,ye,yb,yb/)
        gsn_polyline(wks,plot_fill,x,y,boxres)
    end if
    frame(wks)

end