load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
    
    srcFileName = "triple2grid2d_source.cdf"      ; source (orig) grid
    srcVar      = "Water_Path"

    srcfile     = addfile(srcFileName,"r")
    srclat      = srcfile->Latitude
    srclon      = srcfile->Longitude

    x_src = srcfile->$srcVar$    ; Rank >= 2

    nmsg1 = num(x_src.eq.-888)  ; ***special for this variable***
    nmsg2 = num(x_src.eq.-99)
    print("nmsg1="+nmsg1+"   nmsg2="+nmsg2)

    x_src@_FillValue = 1e20 
    x_src = where(x_src.lt.0, x_src@_FillValue, x_src)
    printVarSummary(x_src)
    printMinMax(x_src, 0)

    opt = True
    opt@PrintStat = True
    statb = stat_dispersion(x_src, opt )    ; look at data

;----------------------------------------------------------------------
; Fill-in the missing values
;----------------------------------------------------------------------
    X_SRC = x_src    ; make a copy for plotting purposes
    poisson_grid_fill(X_SRC, False, 1, 1500, 1e-2, 0.6, 0)

    nrep = 10          ; arbitrary
    do nr=1,nrep
        X_SRC_1 = (/ smth9(X_SRC, 0.50,  0.25, False) /) ; heavy local smoothing
    end do

;----------------------------------------------------------------------
; PLOT                    
;----------------------------------------------------------------------
    plot = new (3, "graphic")

    wks = gsn_open_wks("ps","JAKOB_TEST")
   ;gsn_define_colormap(wks,"wh-bl-gr-ye-re")  ; Change color map
    gsn_define_colormap(wks,"cmp_flux")  ; Change color map

;---Resources to share between both plots
    res                     = True              ; Plot modes desired.
    res@gsnDraw             = False
    res@gsnFrame            = False 
    res@gsnSpreadColors     = True              ; use full range of colormap
    res@gsnAddCyclic        = False

    res@cnFillOn            = True              ; color plot desired
    res@cnLinesOn           = False             ; turn off contour lines
    res@cnLineLabelsOn      = False             ; turn off contour labels
    res@cnFillMode          = "RasterFill"      ; turn raster on      
    res@cnLevelSelectionMode= "ExplicitLevels"
    res@cnLevels            = (/40,80,120,160,200,240,280,320,400,600,800,1000,1200,1400/)

    res@trGridType          = "TriangularMesh"  ; allow missing coordinates
    res@lbLabelBarOn        = False

    res@mpFillOn            = False
    res@mpMinLatF           = min(srclat)
    res@mpMaxLatF           = max(srclat)
    res@mpMinLonF           = min(srclon)
    res@mpMaxLonF           = max(srclon)

    x_src@lat2d = srclat
    x_src@lon2d = srclon
    plot(0)     = gsn_csm_contour_map(wks,x_src  ,res)  

    X_SRC@lat2d = srclat
    X_SRC@lon2d = srclon
    res@gsnLeftString = "after poisson_grid_fill"
    plot(1)     = gsn_csm_contour_map(wks,X_SRC,res)

    X_SRC_1@lat2d = srclat
    X_SRC_1@lon2d = srclon
    res@gsnLeftString = "Repeated smth9: nrep="+nrep
    plot(2)     = gsn_csm_contour_map(wks,X_SRC_1,res)

;---Compare the plots in a panel
    pres                   = True
    pres@gsnMaximize       = True
    pres@gsnPanelLabelBar  = True
    pres@pmLabelBarWidthF  = 0.7
    pres@lbLabelAutoStride = False
    gsn_panel(wks, plot,(/3,1/),pres)