;======================================================================
; ESMF_regrid_29.ncl
;
; Concepts illustrated:
;   - Regrid data on a WRF grid to a rectilinear grid
;   - Regrid data on a rectilinear grid back to the original WRF grid
;   - Generate weights files for regridding use later.
;======================================================================
; This example regrids WRF data to a rectilinear grid, and then takes
; the regridded data and regrids it back to the original WRF grid.
; Three plots are generated and paneled on one page.
;
; This script also generates weight files for use in ESMF_wgts_29
; and other similar scripts.
;======================================================================
; These files are loaded by default in NCL V6.2.0 and newer
; 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/wrf/WRFUserARW.ncl"
;
; This file will be loaded by default in NCL V6.5.0 and newer
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"

begin
;======================================================================
; Generate WRF to Rectilinear regrid weights 
;======================================================================

    InterpMethod= "patch"                      ; define interpolation method

;---Input file
    srcDirName  = "./"
    srcFileName = "wrfout_d01_2003-07-13_12:00:00"
    srcFilePath =  srcDirName + srcFileName

;---Retrieve either one level, or all levels. Use '-1' for all.
    sfile = addfile(srcFilePath,"r")

;---Get the variable to be regridded               
    ua3d = wrf_user_getvar(sfile,"ua",1)       ; On mass grid 
    printVarSummary(ua3d)                      ; (bottom_top,south_north,west_east)

    ua   = ua3d(0,:,:)                         ; (south_north,west_east)
                                               ; only need 'horizontal' spatial points 
    printVarSummary(ua)                        ; (bottom_top,south_north,west_east)
    printMinMax(ua,0)

;---Get the source lat/lon grid (mass grid)
    lat2d = sfile->XLAT(0,:,:)               ; (south_north,west_east)
    lon2d = sfile->XLONG(0,:,:)

    dims  = dimsizes(lat2d)
    nlat  = dims(0)
    nlon  = dims(1)

;---Create the destination rectilinear lat[*]/lon[*] arrays
    lat = fspan(min(lat2d), max(lat2d) ,nlat)
    lon = fspan(min(lon2d), max(lon2d) ,nlon)

;---Set up options for regridding

;---Title and filename options

    Opt                = True
    Opt@SrcTitle       = "WRF grid"   ; optional
    Opt@WgtFileName    = "WRF_to_Rect.WgtFile_"+InterpMethod+".nc"

;---Source grid options
    Opt@SrcFileName    = "WRF.SCRIP_grid_description.nc"         ; Name of source and
    Opt@SrcRegional    = True
    Opt@SrcGridLat     = lat2d
    Opt@SrcGridLon     = lon2d

;---Destination grid options
    Opt@DstFileName    = "Rectilinear.SCRIP_grid_description.nc" ; destination files
    Opt@DstGridType    = "rectilinear"
    Opt@DstGridLat     = lat
    Opt@DstGridLon     = lon
    Opt@DstRegional    = True

;---Specify other options  
    Opt@ForceOverwrite = True
    Opt@InterpMethod   = InterpMethod

;---Perform the regrid: WRF ==> rectilinear (_reclin)
    ua_regrid = ESMF_regrid(ua, Opt)     ; Do the regridding for ua   

;---Reset 0 values to missing values.
    ua_regrid@_FillValue     = default_fillvalue(typeof(ua_regrid))
    ua_regrid                = where(ua_regrid.eq.0.0,ua_regrid@_FillValue,\
                                     ua_regrid)

;---Print rectilinear variable information
    printVarSummary(ua_regrid)
    print("ua_regrid: min="+min(ua_regrid)+"   max="+max(ua_regrid))
    nmsg = num(ismissing(ua_regrid))
    print("nmsgRectilinearGrid="+nmsg)

;======================================================================
; Generate Rectilinear to WRF regrid weights 
;======================================================================

;---For clarity, delete the above options and start again.
    delete(Opt)

;---Options for regridding rectilinear to WRF (curvilinear) grid
    Opt                 = True

;   The grid descriptions have been generated in 'Part A'
;   but we still need to provide their names.
    Opt@SkipSrcGrid     = True
    Opt@SkipDstGrid     = True
    Opt@DstFileName     = "WRF.SCRIP_grid_description.nc"         ; Name of source and
    Opt@SrcFileName     = "Rectilinear.SCRIP_grid_description.nc" ; destination files
    Opt@ForceOverwrite  = True

    Opt@SrcTitle        = srcFileName                  ; source grid
    Opt@SrcMask2D       = where(ismissing(ua_regrid),0,1)
    Opt@SrcRegional     = True

    Opt@DstGridType     = "curvilinear"
    Opt@DstTitle        = "Rectilinear_to_WRF"
    Opt@DstRegional     = True

    Opt@InterpMethod    = InterpMethod
    Opt@WgtFileName     = "Rect_to_WRF.WgtFile_"+InterpMethod+".nc"

    ua_reclin           = ESMF_regrid(ua_regrid,Opt)

;---Print regridded WRF variable information
    printVarSummary(ua_reclin)
    printMinMax(ua_reclin,0)
    nmsg = num(ismissing(ua_reclin))
    print("# missing in regridded WRF grid="+nmsg)

;======================================================================
; Create plots of original data and regridded data
;======================================================================
    ua@lat2d = lat2d    ; Needed for plotting
    ua@lon2d = lon2d

    wks = gsn_open_wks("png","ESMF_regrid")    ; send graphics to PNG file

    plot = new(3,graphic)                      ; create a plot array
    
    res                      = True
    res@gsnDraw              = False           ; don't draw
    res@gsnFrame             = False           ; don't advance frame
    res@gsnAddCyclic         = False           ; regional data
  
    res@cnInfoLabelOn        = False           ; turn off cn info label
    res@cnFillOn             = True            ; turn on color
   ;res@cnFillMode           = "RasterFill"
    res@cnLinesOn            = False
    res@cnLineLabelsOn       = False
    res@lbLabelBarOn         = False           ; turn off individual cb's
  
    res@cnLevelSelectionMode =  "ManualLevels"   
    res@cnMinLevelValF       = -28
    res@cnMaxLevelValF       =  12
    res@cnLevelSpacingF      =  2
    
    res@mpMinLatF            = min(lat2d)      ; range to zoom in on
    res@mpMaxLatF            = max(lat2d)
    res@mpMinLonF            = min(lon2d)
    res@mpMaxLonF            = max(lon2d)
    

    res@gsnLeftString        = ""
    res@gsnRightString       = ""

    res@pmTitleZone          = 4              ; moves title down
    res@tiMainFont           = "Helvetica"    ; set plot title font

    res@pmTickMarkDisplayMode = "Always"      ; nicer map tickmarks
      
    nt = 3                   ; arbitrary time
    kl = 20                  ; arbitrary level
   
    res@tiMainString  = "Original WRF data"
    plot(0) = gsn_csm_contour_map(wks,ua,res) 

    res@tiMainString  = "Regridded to rectilinear"
    plot(1) = gsn_csm_contour_map(wks,ua_reclin,res) 

    res@tiMainString  = "Rectilinear regridded back to WRF"
    plot(2) = gsn_csm_contour_map(wks,ua_reclin,res) 

;======================================================================
; Create panel plot
;======================================================================
    resP                     = True                ; modify the panel plot
    resP@gsnMaximize         = True                ; make ps, eps, pdf large
    resP@gsnPanelMainString  = "Interpolation method used: " + InterpMethod
    resP@gsnPanelMainFont    = "Helvetica-bold"    ; set panel title font
    resP@gsnPanelLabelBar    = True                ; add common colorbar
    resP@pmLabelBarWidthF    = 0.8                 ; increase width of colorbar
    resP@gsnPanelRowSpec     = True                ; 1 plot 1st row, 2 plots 2nd row
    gsn_panel(wks,plot,(/1,2/),resP)               ; now draw as one plot
end