;----------------------------------------------------------------------
; ESMF_rect_to_unstruct.ncl
;
; This is an NCL/ESMF template file for regridding from a 
; rectilinear grid to a unstructured grid. It uses
; ESMF_regrid  to do the regridding.
;
; The assumption is that the rectilinear grid will be represented
; as coordinate arrays of the "var" variable. If this isn't the
; case, you will need to read the lat/lon rectilinear arrays 
; separately and set them via the SrcLat/SrcLon attributes (see
; later code).
;
;
; 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-beta 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
;--Data file containing source grid
    src_file = "/work/MOD3EVAL/css/HTAP_Emissions/annual/TRANSPORT/nox/edgar_HTAP_NOx_emi_TRANSPORT_2010.0.1x0.1.nc" 
    sfile    = addfile(src_file,"r")

;---Get variable to regrid
    varname = "emi_nox"                 
    var     = sfile->$varname$         ; Assumption is that "var"
                                       ; contains coordinate arrays.
    printVarSummary(var)               ; Make sure it has coord arrays.
    printMinMax(var,0)             ; looks okay.

    src_lat  = sfile->lat                     
    src_lon  = sfile->lon                    

;---Data file containing destination grid
;    dst_file = "/home/pou/ncl/x4.163842.grid.nc"                      ;;---Change (likely)
    dst_file = "/work/MOD3DEV/grc/MPAS/MPAS_INIT/INPUTS/90-25/static.90-25.usgs.nc"                      
    dfile    = addfile(dst_file,"r")
    
    dst_lat  = dfile->latCell                     ;;---Change (likely)
    dst_lon  = dfile->lonCell                     ;;---Change (likely)

    r2d     = 180.0d/(atan(1)*4.0d)  ; Radian to Degree
    dst_lon = dst_lon*r2d
    dst_lat = dst_lat*r2d

;---Set up regridding options
    Opt                   = True

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

    Opt@SrcFileName       = "/work/MOD3EVAL/css/HTAP_Emissions/mpas/source_grid_file_htap01.nc"
    Opt@DstFileName       = "/work/MOD3EVAL/css/HTAP_Emissions/mpas/destination_grid_file_90-25.nc"
    Opt@WgtFileName       = "rect_to_unstruct_90-25.nc"

    Opt@Debug             = True

;
; These next two lines only needed if "var" doesn't
; contain coordinate arrays.
;
   Opt@SrcGridLat        = sfile->lat        ;;--Change (maybe)
   Opt@SrcGridLon        = sfile->lon        ;;--Change (maybe)


    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.

    Opt@DstGridType       = "unstructured"
    Opt@DstGridLat        = dst_lat           ; destination grid
    Opt@DstGridLon        = dst_lon
 
    Opt@DstRegional       = False              ;;--Change (maybe)
    Opt@DstMask2D = where(.not.ismissing(dst_lat).and.\
                          .not.ismissing(dst_lon),1,0) ; Necessary if lat/lon
                                                       ; has missing values.

    Opt@DstESMF           = True

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

    Opt@SkipSrcGrid       = False 
    Opt@SkipDstGrid       = False 
    Opt@SkipWgtGen        = False 

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

    printVarSummary(var_regrid)           ; Check that everything
    printMinMax(var_regrid,0)             ; looks okay.

;
;   write out the regridded data
;

    outfile = "/work/MOD3EVAL/css/HTAP_Emissions/mpas/nox_transport_regrid_output_90-25.nc"
    fout     = addfile(outfile,"c")
    varname = "emi_nox"                 
    fout->$varname$   = var_regrid 

;----------------------------------------------------------------------
; Plotting section
;
; This section creates filled contour plots of both the original
; data and the regridded data, and panels them.
;----------------------------------------------------------------------
    wks = gsn_open_wks("ps","rect_to_unstruct")

    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(src_lat)
    res@mpMaxLatF             = max(src_lat)
    res@mpMinLonF             = min(src_lon)
    res@mpMaxLonF             = max(src_lon)

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

;---Resources for plotting regridded data
    res@gsnAddCyclic  = False            ;;---Change (maybe)
    res@tiMainString  = "Unstructured grid (" + Opt@InterpMethod + ")"
    
    res@sfXArray = var_regrid@lon1d
    res@sfYArray = var_regrid@lat1d

    plot_regrid = gsn_csm_contour_map(wks,var_regrid,res)

;---Resources for plotting original data

    delete(res@sfXArray)
    delete(res@sfYArray)

    res@gsnAddCyclic = True            ;;---Change (maybe)
    res@tiMainString = "Original rectilinear 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