;----------------------------------------------------------------
;            Curvilinear-to-Curvilinear (C2C) Regrid Template
; src=>Source [Original] Grid;   dst=>Destination [Target] Grid
;----------------------------------------------------------------
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"
    
;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; Part 1: Grid definitions and Weight Generation
;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; (a) Input files containing grid information
; (b) Specify a variable on the source file. 
;----------------------------------------------------------------------

    srcFileName = "triple2grid2d_source.cdf"      ; source (orig) grid
    dstFileName = "triple2grid2d_destination.grb" ; destination (target)
    srcVar      = "Water_Path"

    method      = "patch"         ; "bilinear"(default), "patch", and "conserve"

;----------------------------------------------------------------
; Grid information only 
;----------------------------------------------------------------

;---Retrieve grid data from Source Grid (Curvilinear)
    srcfile     = addfile(srcFileName,"r")
    srclat      = srcfile->Latitude
    srclon      = srcfile->Longitude

;---Retrieve grid data from Destination Grid (Curvilinear)
    dstfile     = addfile(dstFileName,"r")
    dstlat      = dstfile->g10_lat_0
    dstlon      = dstfile->g10_lon_1

;---SCRIP files that (generally) are generated; Used for input to regrid
    srcGridName = "src_c2c_"+method+"_SCRIP.nc"          ; src grid descripion 
    dstGridName = "dst_c2c_"+method+"_SCRIP.nc"          ; dst grid description
    wgtFile     = "wgt_c2c_"+method+"_SCRIP.nc"

;---Set to True if you want to skip any of these steps
    SKIP_src_SCRIP_GEN = False
    SKIP_dst_SCRIP_GEN = False
    SKIP_wgt_SCRIP_GEN = False

;----------------------------------------------------------------------
; DEBUG: the +360 should not be necessary but it does not hurt
;----------------------------------------------------------------------
    printVarSummary(srclat)
    printMinMax(srclat,0)

    printVarSummary(srclon)
    printMinMax(srclon,0)   ; min=-139.814   max=-57.4228

    srclon = srclon + 360   ; match destination grid (0 - 360)

    printVarSummary(dstlat)
    printMinMax(dstlat,0)

    printVarSummary(dstlon)
    printMinMax(dstlon,0)   ; longitude: min=254.887   max=309.577

;----------------------------------------------------------------------
; Convert source grid to SCRIP file.
;----------------------------------------------------------------------
   ;x2d = srcfile->$srcVar$    ; [*][*]   possibly used for mask

    if(.not.SKIP_src_SCRIP_GEN) then
;---Convert to an SCRIP Convention file.
      Opt                = True
      Opt@ForceOverWrite = True
      Opt@PrintTimings   = True
      if (isvar("x2d")) then
          Opt@Mask2D     = where(.not.ismissing(x2d),1,0)
      end if

      curvilinear_to_SCRIP(srcGridName,srclat,srclon,Opt)
      delete(Opt)        ;---Clean up
    end if

;----------------------------------------------------------------------
; Convert destination curvilinear grid to SCRIP file
;----------------------------------------------------------------------

    if(.not.SKIP_dst_SCRIP_GEN) then
      Opt                = True
      Opt@ForceOverWrite = True
      Opt@PrintTimings   = True
      Opt@Title          = "GRIB curvilinear grid."

      curvilinear_to_SCRIP(dstGridName,dstlat,dstlon,Opt)
      delete(Opt)        ;---Clean up
    end if

;----------------------------------------------------------------------
; Generate interpolation weights for  source grid to destination grid  
;----------------------------------------------------------------------
    if(.not.SKIP_wgt_SCRIP_GEN) then
      Opt                      = True
      Opt@SrcESMF              = False
      Opt@DstESMF              = False
      Opt@ForceOverWrite       = True
      Opt@PrintTimings         = True
      Opt@method               = method 

      ESMF_regrid_gen_weights(srcGridName, dstGridName, wgtFile, Opt)
      delete(Opt)        ;---Clean up
    end if

;----------------------------------------------------------------------
; Part 2:  Evaluation 
; Regrid variable from source grid to destination grid.
;----------------------------------------------------------------------

    Opt              = True
;   Opt@Debug        = True
    Opt@PrintTimings = True

    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 = -9999.
    x_src = where(x_src.lt.0, x_src@_FillValue, x_src)
    printVarSummary(x_src)
    printMinMax(x_src, 0)

    x_regrid = ESMF_regrid_with_weights(x_src ,wgtFile,Opt)
    copy_VarAtts_except(x_src  ,x_regrid,"_FillValue")  ; ????
    printVarSummary(x_regrid)
    printMinMax(x_regrid, 0)

;----------------------------------------------------------------------
; PLOT                    
;----------------------------------------------------------------------
    wks = gsn_open_wks("ps","ESMF_c2c_"+method)
    gsn_define_colormap(wks,"rainbow")             ; Change color map

;---Resources to share between both plots
    res                     = True              ; Plot modes desired.
    res@gsnDraw             = False
    res@gsnFrame            = False 
  ;;res@gsnMaximize         = True              ; Maximize plot
    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            = ispan(270,300,2)

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

    res@mpFillOn            = False
    res@mpMinLatF           = min((/min(srclat),min(dstlat)/))-1
    res@mpMaxLatF           = max((/max(srclat),max(dstlat)/))+1
    res@mpMinLonF           = min((/min(srclon),min(dstlon)/))-1
    res@mpMaxLonF           = max((/max(srclon),max(dstlon)/))+1

;---Regridded data
    res@tiMainString = "REGRID "+method+" (" + \
                       str_join(tostring(dimsizes(x_regrid))," x ") + \
                       ") (SCRIP)"

    x_regrid@lat2d   = dstlat
    x_regrid@lon2d   = dstlon
    plot_regrid = gsn_csm_contour_map(wks,x_regrid,res)


;---Original data
    x_src@lat2d      = srclat
    x_src@lon2d      = srclon
    res@tiMainString = "SOURCE (" + \
                      str_join(tostring(dimsizes(x_src))," x ")+")"

    plot_orig = gsn_csm_contour_map(wks,x_src  ,res)  

;---Compare the plots in a panel
    pres                   = True
    pres@gsnMaximize       = True
    pres@gsnPanelLabelBar  = True
    pres@lbLabelAutoStride = True

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