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

;;DIAG_NORCPM; FIGFILENAME: fig
;;DIAG_NORCPM; TITLE: fig title
;;DIAG_NORCPM; PLOTRES:  ;; plot resource
;;DIAG_NORCPM; FILE1AVG: 
;;DIAG_NORCPM; FILE1AVGVAR: INVAR1
;;DIAG_NORCPM; FILE1STD: file1std.not.set
;;DIAG_NORCPM; FILE1STDVAR: INVAR1
;;DIAG_NORCPM; INVAR1: 
;;DIAG_NORCPM; NUMBEROF1: 0
;;DIAG_NORCPM; FILE2AVG: 
;;DIAG_NORCPM; FILE2AVGVAR: INVAR2
;;DIAG_NORCPM; FILE2STD: file2std.not.set
;;DIAG_NORCPM; FILE2STDVAR: INVAR2
;;DIAG_NORCPM; INVAR2: 
;;DIAG_NORCPM; NUMBEROF2: 0
;;DIAG_NORCPM; PSIGLVL: 0.05
;;DIAG_NORCPM; ISOCNGRID: False


begin
    title = "FLUT NOS-SON dif"
    figfn = "OLRTOA_dif"
    ;; read avg, should be (y,x)
    f1 = addfile("mean_noresm2_lm-snow-rad02_OLRTOA.nc","r")
    avg1 = f1->FLUT
    f2 = addfile("mean_CTRL_OLRTOA.nc","r")
    avg2 = f2->FLUT

    davg = avg1
    davg = avg1 - avg2

    dims1 = dimsizes(avg1)
    dims2 = dimsizes(avg2)
    if (any(dims1.ne.dims2))then
        print(dims1)
        print(dims2)
        print("dims does not match,exit... ")
        exit
    end if

    ;; read stddev if presented, should be (y,x)
    shading = False
    if (isfilepresent("anovar_noresm2_lm-snow-rad02_OLRTOA.nc").and.isfilepresent("anovar_CTRL_OLRTOA.nc"))then
        f1 = addfile("anovar_noresm2_lm-snow-rad02_OLRTOA.nc","r")
        std1 = f1->FLUT
        f2 = addfile("anovar_CTRL_OLRTOA.nc","r")
        std2 = f2->FLUT
        ;; number of elements used to cal stddev and avg
        n1 = (2014-1980+1)*12
        n2 = (2014-1980+1)*12
        var1 = std1
        var1 = std1*std1
        var2 = std2
        var2 = std2*std2
        shading = True
    end if


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

    if(shading)then
        psig = 0.05                     ;; test significance level
        iflag = True
        tval_opt = False                   ;; False for p-value only, True for additional t-value
        prob = ttest(avg1,var1,n1,avg2,var2,n2,iflag,tval_opt)
        copy_VarCoords(davg,prob)
    end if

    ;; 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@mpCenterLonF = 210.  ;; do not cut ocean basin
        


    if (shading)then
    res2 = True ;; for shade on filled contour
        res2@cnLinesOn = False
        res2@cnLineLabelsOn = False
        ;;res2@cnFillMode = "CellFill"
        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(prob)/)
        
    resShade = True
        resShade@gsnShadeFillType = "pattern"
        resShade@gsnShadeHigh   = 17
        resShade@gsnShadeLow = 17
    end if

    if (False)then
        gridf = addfile("/nird/projects/NS9039K/shared/pgchiu/diag_norcpm/grid_norcpm1_ocn.nc","r")
        plat = gridf->plat
        plon = gridf->plon
        davg@lon2d := plon
        davg@lat2d := plat
        if(shading)then
            prob@lon2d = plon
            prob@lat2d = plat
        end if
    end if

    latbe = (/-40.,60./)
    lonbe = (/120.,300./)
    res@mpMinLatF = latbe(0)
    res@mpMaxLatF = latbe(1)
    res@mpMinLonF = lonbe(0)
    res@mpMaxLonF = lonbe(1)
    res@cnFillPalette = "hotcold_18lev"
res@cnLevelSelectionMode = "ExplicitLevels"
res@cnLevels = (/-15,-10,-5,-3,-2,2,3,5,10,15/)

    

    wks = gsn_open_wks("ps",figfn)

    plot_fill = gsn_csm_contour_map(wks,davg,res)

    if(shading)then
        plot_shade = gsn_csm_contour(wks,prob,res2)
        plot_shade = gsn_contour_shade(plot_shade,psig,max(prob),resShade)
        overlay(plot_fill,plot_shade)
    end if
    draw(plot_fill)
    frame(wks)

end