;----------------------------------------------------------------------
; ESMF_emissions_FV_96x144.ncl
;
; Concepts illustrated:
;   - Interpolating from one rectilinear grid to another using ESMF_regrid
;   - Create netCDF containing regridded variable(s) 
;----------------------------------------------------------------------
;---Create netCDF containing the regridded variable(s)
    netCDF = True

;---Specify remap method  
    remap_method = "conserve"     ; "bilinear" , "patch", "conserve"

;---Source grid information 
    SrcDir   = "./"
    SrcName  = "emissions-cmip6_so4_a1_contvolcano_vertical_2000-2015.nc"
    SrcPath  = SrcDir + SrcName
    SrcVar   = "emiss_volcanoes"    ; any variable 

    sfile    = addfile(SrcPath,"r")
    lat_in   = sfile->lat
    lon_in   = sfile->lon
    var_in   = sfile->$SrcVar$

;---Data file containing destination grid
;---or 
;---Create destination latitude and longitude 

    NLAT25x19    =  96                              ; RES   = "2.5x1.9"
    MLON25x19    = 144
    LAT25x19     = latGau    (NLAT25x19, "lat", "latitude" , "degrees_north")
    LON25x19     = lonGlobeFo(MLON25x19, "lon", "longitude", "degrees_east" )

;---Destination grid information: Will contain regridded data

    DstDir   = "./"
    DstGrid  = "FV_96x144"
    DstRoot  = "Emissions_CMIP6_to"
    DstPath  = DstDir + DstRoot+"_"+DstGrid  

;---Specify name of weight file to be generated; name of destination grid; dst directory

    WgtDir   = "./"
    WgtName  = WgtDir+DstRoot+"_"+DstGrid+".Weight."+remap_method+".nc"

;---Specify name of file containing the regridded variable(s) 

    if (netCDF) then
        RgdDir = "./" 
        RgdName= RgdDir+DstRoot+"_"+DstGrid+".Regrid."+remap_method+".nc"
        RgdPath= RgdDir + RgdName
    end if

;---Sample plot options
    pltDir       = "./"
    pltType      = "png"
    pltPath      = pltDir+DstRoot+"_"+DstGrid

;---Specify any variable to regrid

;+++
; End user input
;+++

;---Set up regridding options
    Opt                = True
    Opt@InterpMethod   = remap_method

    Opt@WgtFileName    = WgtDir+WgtName
    Opt@SrcGridLat     = lat_in                      ; [*]
    Opt@SrcGridLon     = lon_in                      ; [*]
    Opt@SrcInputFileName = SrcPath   
    Opt@RemoveSrcFile  = True      ; remove SCRIP source grid description file

    Opt@DstGridLat     = LAT25x19
    Opt@DstGridLon     = LON25x19
    Opt@RemoveDstFile  = True      ; remove SCRIP destination grid description file

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

;---Call the regridding function

    var_regrid = ESMF_regrid(var_in,Opt)
    printVarSummary(var_regrid)
    printMinMax(var_regrid,0)
    print("----------------")

;----------------------------------------------------------------------
; Write the regridded variable to a netCDF file 
;----------------------------------------------------------------------
    if (netCDF) then
        var_regrid@ESMF_method = remap_method

        system("rm -f " + RgdPath)
        ncdf = addfile(RgdPath,"c")

        fAtt               = True            ; assign file attributes
        fAtt@title         ="ESMF Regrid: "+remap_method
        fAtt@source_file   = SrcName               
        fAtt@Conventions   = "None"
        fAtt@creation_date = systemfunc ("date")
        fAtt@NCL           = "NCAR Command Language: "+get_ncl_version()
        fileattdef( ncdf, fAtt )             ; copy file attributes

        filedimdef(ncdf,"time",-1,True)      ; make time "UNLIMITED" 
        ncdf->$SrcVar$     = var_regrid      ; Write regridded variable
    end if

;----------------------------------------------------------------------
; Plotting section
; This section creates filled contour plots of both the original
; data and the regridded data, and panels them.
;----------------------------------------------------------------------

    dims_in     = dimsizes(var_in)
    rank_in     = dimsizes(dims_in)
    ntim        = dims_in(0)

    dims_regrid = dimsizes(var_regrid)
    rank_regrid = dimsizes(dims_regrid)

    if (rank_in.gt.3) then
        klev = dims_in(1)     
        nt   = 100                                ; arbitrary
        kl   = 10                                 ; arbitrary
        LEV  = var_in&altitude(kl)
    end if

    wks = gsn_open_wks(pltType,pltPath)

;---Resources to share between both plots
    res                      = True              ; Plot modes desired.
    res@gsnDraw              = False
    res@gsnFrame             = 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@cnFillPalette        = "WhiteBlueGreenYellowRed"               

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

    res@gsnLeftString        = SrcVar            ; long_name is too long!
    res@gsnRightString       = var_in@units      

;---Resources for plotting regridded data

    res@tiMainString  = "Destination_grid: (" + \
                        str_join(tostring(dims_regrid(rank_regrid-2:))," x ") + ") using '" + \
                        Opt@InterpMethod + "' method"

    if (rank_in.gt.3) then
        res@gsnCenterString      = "level="+LEV
        plot_regrid = gsn_csm_contour_map(wks,var_regrid(nt,kl,:,:),res)
    else 
        plot_regrid = gsn_csm_contour_map(wks,var_regrid(nt,:,:),res)
    end if

;---Get resources used for contouring the above grid; apply to next grid

    res@cnLevelSelectionMode = "ManualLevels"
    getvalues plot_regrid@contour
      "cnMinLevelValF"       : res@cnMinLevelValF   
      "cnMaxLevelValF"       : res@cnMaxLevelValF 
      "cnMaxLevelValF"       : res@cnMaxLevelValF 
      "cnLevelSpacingF"      : res@cnLevelSpacingF
    end getvalues

    if (rank_in.ge.3) then
        res@tiMainString = "Source grid (" + str_join(tostring(dims_in(rank_in-2:))," x ") +")"

        res@gsnCenterString      = "level="+LEV
        plot_in = gsn_csm_contour_map(wks,var_in(nt,kl,:,:),res)
    else 
        plot_in = gsn_csm_contour_map(wks,var_in(nt,:,:),res)
    end if

;---Draw both plots in a panel
    pres                   = True
    pres@gsnMaximize       = True
    pres@gsnPanelLabelBar  = True
    gsn_panel(wks,(/plot_in,plot_regrid/),(/2,1/),pres)