; This is the *only* file (library) that has to be explicitly loaded
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"

;********************************************************
;********************************************************

    src_dir  = "./"
    src_file = "maxele.63.nc"
    src_path = src_dir + src_file
    sfile = addfile(src_path,"r")

    lat1D = sfile->y                  ; y(node); node=1813443
    lon1D = sfile->x      
    vertices = sfile->element-1

    var   = sfile->zeta_max           ; var(node), type=>double"
    printVarSummary(var)
    printVarSummary(lat1D)
    printMinMax(lat1D,0)
    printMinMax(lon1D,0)

    dst_dir  = "./"
    dst_file = "geo_em.d01.nc"
    dst_path = dst_dir + dst_file
    dfile    = addfile(dst_path,"r")

    XLAT_M               = dfile->XLAT_M(0,:,:)    ; Will be useful later              
    XLONG_M              = dfile->XLONG_M(0,:,:) 
    printVarSummary(XLAT_M)
    printVarSummary(XLONG_M)
    printMinMax(XLAT_M,0)
    printMinMax(XLONG_M,0)

;---Set up regridding options

    Opt                = True

;---"bilinear" is the default. "patch" and "conserve" are other options.
    Opt@InterpMethod     = "bilinear"        ;;---Change (maybe)
    Opt@WgtFileName      = "unstruct_to_rect.nc"

    Opt@SrcGridLat       = lat1D
    Opt@SrcGridLon       = lon1D
    Opt@SrcRegional      = True
    Opt@SrcTriangularMesh = sfile->element-1
    Opt@SrcInputFileName = src_file          ; optional, but good idea

    Opt@SrcMask2D        = where(.not.ismissing(var),1,0) ; Necessary if has
                                                          ; missing values.

    Opt@DstGridLon       = XLONG_M 
    Opt@DstGridLat       = XLAT_M 
    Opt@DstRegional      = True            ; WRF taget grid is regional
;   Opt@DstMask2D = where(.not.ismissing(dst_lat).and.\
;                         .not.ismissing(dst_lon),1,0) ; Necessary if lat/lon
                                                      ; has missing values.
                                                      ; has missing values.
    Opt@ForceOverwrite   = True
    Opt@Debug            = True
    Opt@PrintTimings     = True

    var_regrid = ESMF_regrid(var,Opt)     ; Do the regridding

    var_regrid@long_name  = "maximum water surface elevation above geoid"  ; assign attributes 
    var_regrid@units      = "m"
    var_regrid@_FillValue = var_regrid@missing_value
    var_regrid@coordinates= "XLAT_M XLONG_M"

    delete( [/ var_regrid@lat2d, var_regrid@lon2d \
             , var_regrid@mesh,  var_regrid@location /] ); not appropriate or needed

    copy_VarCoords(XLAT_M, var_regrid)
    printVarSummary(var)
    printVarSummary(var_regrid)
    printMinMax(var,0)
    printMinMax(var_regrid,0)

;----------------------------------------------------------------------
; Plotting section
;
; This section creates filled contour plots of both the original
; data and the regridded data, and panels them.
;----------------------------------------------------------------------
    var_regrid@lat2d = XLAT_M                    ; @lat2d  for plotting only
    var_regrid@lon2d = XLONG_M                   ; @lon2d   "     "      "
   
    wks = gsn_open_wks("png","unstruct_to_rect")

;---Resources to share between both plots
    res                      = True              ; Plot mods desired
    res@gsnDraw              = False
    res@gsnFrame             = False

    res@gsnMaximize          = True              ; Maximize plot

    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@mpMinLatF            = min(XLAT_M)
    res@mpMaxLatF            = max(XLAT_M)
    res@mpMinLonF            = min(XLONG_M)
    res@mpMaxLonF            = max(XLONG_M)

   res@cnLevelSelectionMode = "ManualLevels"
   res@cnMinLevelValF       = -6
   res@cnMaxLevelValF       = 6
   res@cnLevelSpacingF      = 0.25

    res@lbLabelBarOn         = False      ; Will turn on in panel later

    res = wrf_map_resources(dfile,res)
    res@gsnAddCyclic          = False            ; regional data: not cyclic 
    res@tfDoNDCOverlay        = True             ; set True for native projections

    res@mpFillOn                    = False
    res@mpGeophysicalLineColor      = "black"    ; wrf_map_resources uses "gray"
    res@mpUSStateLineColor          = "black"
    res@mpGeophysicalLineThicknessF = 2.0        ; wrf_map_resources uses 0.5
    res@mpUSStateLineThicknessF     = 2.0

;---Resources for plotting regridded data
    res@gsnAddCyclic = False      
    res@tiMainString  = "Data on WRF grid (" + Opt@InterpMethod + ")"

    plot_regrid = gsn_csm_contour_map(wks,var_regrid,res)

;---Resources for plotting original data
    res@sfXArray     = lon1D
    res@sfYArray     = lat1D
    res@tiMainString = "Original unstructured grid (" + \
                        dimsizes(lon1D) + " cells)"

    plot_orig = gsn_csm_contour_map(wks,var,res)  

;---Draw both plots in a panel
    pres                   = True
    pres@gsnMaximize       = True
    pres@gsnPanelLabelBar  = True

    gsn_panel(wks,(/plot_orig,plot_regrid/),(/2,1/),pres)