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

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

begin
    ;; plot sst difference (2000-2020) - (1980-2000)
    nmem = 5
    ocngrid = "/nird/projects/NS9039K/shared/pgchiu/diag_norcpm/grid_tnx1v4_20170622.nc"

    ;; get dsst of all members
    do i = 1,nmem
        ii = sprinti("%2.2d",i)
        infile1 = "../../N2_FF_w2_1980-2000/01_sst_Atl3-cor/sst_"+ii+".nc"
        infile2 = "../../N2_FF_w2_2000-2020/01_sst_Atl3-cor/sst_"+ii+".nc"
        f = addfile(infile1,"r")
        sst1 := f->sst
        sst1 := dim_avg_n_Wrap(sst1,0)
        f = addfile(infile2,"r")
        sst2 := f->sst
        sst2 := dim_avg_n_Wrap(sst2,0)
        if(.not.isvar("dsst"))then
            dims = dimsizes(sst1)
            dsst = conform_dims(array_append_record(nmem,dims,0),sst1,(/1,2/))
            dsst(0,:,:) = sst1
            dsst!0 = "member"
            dsst&member = ispan(1,nmem,1)
        end if
        dsst(i-1,:,:) = sst2-sst1
    end do

    ;; assign coordinate to sst
    f = addfile(ocngrid,"r")
    dsst@lat2d = f->plat
    dsst@lon2d = f->plon

    ;; common res
    res = True
        res@cnFillOn    = True
res@cnFillMode  = "CellFill"
res@cnLinesOn   = False
res@lbOrientation     = "Horizontal"
res@gsnAddCyclic = True
;res@gsnRightString = "p siglvl=0.05"
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.
res@mpFillDrawOrder = "PostDraw"
res@mpCenterLonF = 210.


    ;; plot dsst of each member
    do i = 1, nmem
        ii = sprinti("%2.2d",i)
        wks = gsn_open_wks("ps","dsst_N2_FF_w2_"+ii)
        res@gsnLeftString = "N2_FF_w2 mem "+ii+"dsst 2000-2020 - 1980-2000"
        plot = gsn_csm_contour_map(wks,dsst(i-1,:,:),res)
    end do
    ;; plot ensemble mean dsst 
    wks = gsn_open_wks("ps","dsst_N2_FF_w2_ens")
    dsst_ens = dim_avg_n_Wrap(dsst,0)
    res@gsnLeftString = "N2_FF_w2 ensmean dsst 2000-2020 - 1980-2000"
    plot = gsn_csm_contour_map(wks,dsst_ens,res)
end