; 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")
;print(sfile)
    lat1D = sfile->y                  ; y(node); node=1813443
    lon1D = sfile->x      
    var   = sfile->zeta_max           ; var(node), type=>double"
    printVarSummary(var)

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

;---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      = False             ;;--Change (maybe)
    Opt@SrcInputFileName = src_file          ; optional, but good idea

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

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

    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 elevationabove 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_regrid)

;---Code for writing netcdf format file with lat and long cordinate

    setfileoption("nc", "Format",  "NetCDF4")
    ncDir  = "./"
    ncName = "ZETA_MAX_D01.nc"
    ncPath = ncDir+ncName
    system("/bin/rm -f " + ncPath) ; remove if exists
    fnew = addfile(ncPath, "c")

;===================================================================
; create global attributes of the file (optional)
;===================================================================
    fAtt               = True            ; assign file attributes
    fAtt@title         = "Unstructured ADCIRC Grid to WRF Grid Creation"
    fAtt@source_file   =  sfile
    fAtt@Conventions   = "None"
    fAtt@NCL_tag       = "ESMF: "+Opt@InterpMethod 
    fAtt@creation_date = systemfunc ("date")
    fileattdef( fnew, fAtt )            ; copy file attributes

    fnew->ZETA_MAX = var_regrid                ; name on file different from name in code
    fnew->XLAT_M   = XLAT_M           
    fnew->XLONG_M  = XLONG_M          
 
;----------------------------------------------------------------------
; 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)

   ;mnmxint = nice_mnmxintvl( 0, max(var), 18, False)
   ;res@cnLevelSelectionMode = "ManualLevels"
   ;res@cnMinLevelValF       = mnmxint(0)
   ;res@cnMaxLevelValF       = mnmxint(1)
   ;res@cnLevelSpacingF      = mnmxint(2)

    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)