;----------------------------------------------------------------------
; ESMF_curv_to_curv.ncl
;
; This is an NCL/ESMF template file for regridding from one
; curvilinear grid to another curvilinear grid. It uses 
; ESMF_regrid to do the regridding.
;
; Both grids are assumed to be contained in separate NetCDF files.
; 
; Search for lines with ";;---Change (likely)" or ";;---Change (maybe)".
; These are the lines you will likely or maybe have to change.
;
; Of course, you'll probably want to change other aspects of this
; code, like the options for plotting (titles, colors, etc).
;
; For more information on ESMF_regrid, see: 
;
;    http://www.ncl.ucar.edu/Document/Functions/ESMF/ESMF_regrid.shtml
;----------------------------------------------------------------------
; This example uses the ESMF application "ESMF_RegridWeightGen" to 
; generate the weights.
;
; For more information about ESMF:
;
;        http://www.earthsystemmodeling.org/
;
; This script uses built-in functions that are only available in 
; NCL V6.1.0 and later.
;----------------------------------------------------------------------
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"

begin
    latS = 40
    latN = 46
    lonW = -80
    lonE = -71
;---Data file containing source grid
    diri = "./"
    src_file = systemfunc("cd "+diri+" ; ls MCD19A2.A2010029*.hdf")
    fili = systemfunc("cd "+diri+" ; ls MCD19A2.A2010029*.hdf")
    sfile    = addfile(diri+fili+".he2", "r")

;---Get variable to regrid
    var_name = "Optical_Depth_047_grid1km"           ;;---Change (likely)
    var      = short2flt(sfile->$var_name$)            ;;---Change (likely)
    var      = where(var.ge. -0.1 .and. var.le.5.0, var@FillValue, var) 
    src_lat  = sfile->GridLat_grid1km                   ;;---Change (maybe)
    src_lon  = sfile->GridLon_grid1km                   ;;---Change (maybe)

;---Data file containing destination grid
;    dst_file = systemfunc("cd "+diri+" ; ls MOD04_L2*.hdf")                    ;;---Change (likely)
;    dfili = systemfunc("cd "+diri+" ; ls MOD04_L2*.hdf")
;    dfile    = addfile(diri+dfili+".he2", "r")

;---Set up regridding options
    Opt                   = True

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

    Opt@SrcGridLat        = src_lat           ; source grid
    Opt@SrcGridLon        = src_lon
    Opt@SrcRegional       = True              ;;--Change (maybe)
    Opt@SrcInputFileName  = src_file          ; optional, but good idea
    Opt@SrcMask2D         = where(.not.ismissing(var),1,0) ; Necessary if has
                                                           ; missing values.
;****************************************************************************
;changing grid description
;***************************************************************************
    Opt@DstGridType = "0.5deg" ; destination grid description
    Opt@DstTitle    = "Global 0.5 degree resolution"
    Opt@DstLLCorner = (/ latS, lonW /)
    Opt@DstURCorner = (/ latN, lonE /)
    Opt@DstRegional       = True              ;;--Change (maybe)
;                                                      ; has missing values.
    Opt@DstMask2D = where(.not.ismissing(var) ,1,0)

    Opt@ForceOverwrite    = True
    Opt@PrintTimings      = True
    Opt@Debug             = True

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

    printVarSummary(var_regrid)

;----------------------------------------------------------------------
; Plotting section
;
; This section creates filled contour plots of both the original
; data and the regridded data, and panels them.
;----------------------------------------------------------------------
    var@lat2d = src_lat     ; Needed for plotting.
    var@lon2d = src_lon     ; var_regrid already has these

    wks = gsn_open_wks("ps","curv_to_curv")

    res                       = True

    res@gsnMaximize           = True

    res@gsnDraw               = False
    res@gsnFrame              = False

    res@cnFillOn              = True
    res@cnLinesOn             = False
    res@cnLineLabelsOn        = False
    res@cnFillMode            = "RasterFill"

    res@lbLabelBarOn          = False    ; Turn on later in panel

    res@mpMinLatF             = min(dst_lat)
    res@mpMaxLatF             = max(dst_lat)
    res@mpMinLonF             = min(dst_lon)
    res@mpMaxLonF             = max(dst_lon)

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

;---Resources for plotting regridded data
    res@gsnAddCyclic = False            ;;---Change (maybe)
    res@tiMainString = "Curvilinear grid (" + Opt@InterpMethod + ")"

    plot_regrid = gsn_csm_contour_map(wks,var_regrid,res)

;---Resources for plotting original data
    res@gsnAddCyclic = False            ;;---Change (maybe)
    res@tiMainString = "Original curvilinear grid"

    plot_orig = gsn_csm_contour_map(wks,var,res)  

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

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