;
; The first three loads are not needed in NCL V6.4.0
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
;;  src_dir  = "/home/umalmonj/"
;;  cnt_dir  = "/global/scratch/umalmonj/WRF/juris/CTRL/2D/"
    src_dir  = "./"
    cnt_dir  = "./"

;---Data file containing source grid
    src_file      = "Jan_Apr_monthly_total_2d.nc"
    src_constants = "wrfout_conus_constants.nc"
    sfile         = addfile(src_dir+src_file,"r")
    constants     = addfile(cnt_dir+src_constants,"r")

;---Get variable to regrid
    prec_acc_nc = sfile->PREC_ACC_NC
    var = dim_sum_n_Wrap(prec_acc_nc,0) ; this sums the variable for the 4 months
    src_lat  = constants->XLAT(0,:,:)
    src_lon  = constants->XLONG(0,:,:)
    src_lon = where(src_lon.lt.0,src_lon+360,src_lon)

;---Data file containing destination grid
;;  dst_dir  = "/global/scratch/umalmonj/WRF/juris/CaPA/"
    dst_dir  = "./"
    dst_file = "AccumulatedtotalPrecipitation.nc"
    dfile    = addfile(dst_dir+dst_file,"r")
    dst_lat  = dfile->gridlat_0
    dst_long = dfile->gridlon_0
    dst_lon  = where(dst_long.lt.0,dst_long+360,dst_long)

;---"bilinear" is the default. "patch" and "conserve" are other options.
    interp_method    = "conserve"       ; might want to use conservative for precip
    wgt_filename     = "WRF_to_CMC_RDPA_APCP_" + interp_method + ".nc"

;---Set up regridding options
    Opt                   = True

    Opt@InterpMethod      = interp_method
    Opt@WgtFileName       = wgt_filename

;---Source grid
    Opt@SrcGridLat        = src_lat
    Opt@SrcGridLon        = src_lon
    Opt@SrcRegional       = True
    Opt@SrcGridMask       = where(.not.ismissing(var),1,0) ; Necessary if source grid has msg vals

;---Destination grid
    Opt@DstGridLat        = dst_lat
    Opt@DstGridLon        = dst_lon
    Opt@DstRegional       = True
    Opt@DstGridMask = where(.not.ismissing(dst_lat).and.\
                          .not.ismissing(dst_lon),1,0) ; Necessary if destination lat/lon has msg vals
;
; there is an attribute in the dest. file called "corners", but there are four
; diff. numbers for both the lon and lat variables
; just let NCL figure out the corners...
;    Opt@DstLLCorner = (/12.21208,217.1075/)
;    Opt@DstURCorner = (/79.82558,147.6248/)
;
    Opt@ForceOverwrite    = True
    Opt@PrintTimings      = True
    Opt@Debug             = True

;;; USE THIS CALL IF YOU ALREADY HAVE THE WEIGHTS FILE!!
;;  var_regrid = ESMF_regrid_with_weights(var,wgt_filename,True) ; Do the regridding using a weights file

    var_regrid = ESMF_regrid(var,Opt)     ; Create the weights file and 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("x11","WRF_regrid_lc_"+interp_method)

    res                       = True

    res@gsnMaximize           = True

    res@gsnDraw               = False
    res@gsnFrame              = False

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

    res@lbLabelBarOn          = False    ; Turn on later in panel

    res@cnLevelSelectionMode = "ExplicitLevels"
    res@cnLevels             = (/0,50,100,200,300,500,750,1000,1100,1200,\
                                 1300,1400,1500,1600,1800,2000,2200,2400,\
                                 2500,2600,2700,2800/)

    res@gsnAddCyclic         = False
    res@mpProjection         = "LambertConformal"  ; choose projection
    res@mpLambertParallel1F  = 33.0                ; two parallels
    res@mpLambertParallel2F  = 45.0
    res@mpLambertMeridianF   = -95.0               ; central meridian
    res@mpLimitMode          = "LatLon"
    res@mpMinLatF            = 15.0                ; map area
    res@mpMaxLatF            = 80.0                ; latitudes
    res@mpMinLonF            = -150.0              ; and
    res@mpMaxLonF            = -40.0               ; longitudes
    res@gsnLeftString        = "Accumulated precip"    ; var@long_name is a bit too long

;---Resources for plotting regridded data
    res@tiMainString = "Curvilinear grid (" + Opt@InterpMethod + ")"
    plot_regrid       = gsn_csm_contour_map(wks,var_regrid,res)

;---Resources for plotting original data
    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/),(/1,2/),pres)

;---Add some grid lines showing the destination grid and panel again
    nstep                 = 10      ; If you draw every line, it will be too dense!
    gres                  = True
    gres@gsLineColor      = "HotPink"
    gres@gsLineThicknessF = 2.0
    gres@gsnCoordsLat     = dst_lat(::nstep,::nstep)
    gres@gsnCoordsLon     = dst_lon(::nstep,::nstep)
    gres@gsnCoordsAsLines = True
    gres@gsnCoordsAttach  = True 
    gsn_coordinates(wks,plot_orig,var_regrid(::nstep,::nstep),gres)
    gsn_coordinates(wks,plot_regrid,var_regrid(::nstep,::nstep),gres)

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