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

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

procedure plotice(n)
begin
    season = getenv("SEASON")
    trendtitle = season+" trend dif. "+getenv("DIFSTR")

    if(n.eq.1)then
        vn = "aice"
        title = "Sea ice fraction "+trendtitle
        unit = "%/10yr"
        cnlevels = fspan(-8.,8.,9)
    end if
    if(n.eq.2)then
        vn = "hi"
        title = "Sea ice thickness "+trendtitle
        unit = "m/10yr"
        cnlevels = fspan(-.5,.5,9)
    end if

    case_tropac = getenv("CASE1")
    case_global = getenv("CASE2")
    figfn = "dif_"+getenv("DIFSTRFN")+"_"+vn+"_"+season+"_trend_ensmean"
    fn1 = case_tropac+"_ens_"+vn+"_"+season+"_trend.nc"
    fn2 = case_global+"_ens_"+vn+"_"+season+"_trend.nc"

    f = addfile(fn1,"r")
    var1 = f->$vn$
    var1 = var1*10 ;; annual trend to decadal

    f = addfile(fn2,"r")
    var2 = f->$vn$
    var2 = var2*10 ;; annual trend to decadal

    dvar = var1
    dvar = var1-var2

    if (True)then
        gridf = addfile("~/pgc_NS9039K/gridfiles/grid_gx1v6.nc","r")
        plat = gridf->plat
        plon = gridf->plon
        dvar@lat2d = plat
        dvar@lon2d = plon
    end if

    res = True
        res@gsnPolar = "NH"
        res@cnFillOn = True
        res@mpMinLatF = 40.
        res@gsnLeftString = title
        res@gsnRightString = unit
        ;res@cnFillPalette = "BlueDarkRed18"
        cmap = read_colormap_file("BlueDarkRed18")
        cmap = cmap(::-1,:)
        res@cnFillPalette = cmap


        res@cnLevelSelectionMode = "ExplicitLevels"
        res@cnLevels = cnlevels


    wks = gsn_open_wks("ps",figfn)
    plot = gsn_csm_contour_map_polar(wks,dvar(0,:,:),res)
end

begin
    plotice(1)
    plotice(2)
end