; Script for regridding from WRF curvelinear grid to rectlinear grid of gpm data

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

;---Read the source and destination lat/lon grids
 
    srcFileName = "wrf.nc" 
    dstFileName = "precip.nc" 
    sfile = addfile(srcFileName,"r") 
    dfile = addfile(dstFileName,"r") 
    src_lat = sfile->XLAT(8,:,:) 
    src_lon = sfile->XLONG(8,:,:) 

    wrf_times = wrf_user_list_times(sfile)         ; "2008-09-29_18:30:00", etc

    dst_lat = dfile->lat 
    dst_lon = dfile->lon 

;---Read the data you want to regrid 
;    r1 = sfile->RAINNC
;    r2 = sfile->RAINC
;    rain = r1 + r2
;    var@long_name = "Total Rain"            
;    var@units = r1@units
;    copy_VarCoords(r1, var)

;---For test purposes, use a variable that has more "coverage" (than, say, precipitation)
    var = wrf_user_getvar(sfile,"T2",-1)
 
;---Options for regridding 
    Opt = True 

;---source grid information 
    Opt@SrcRegional = True 
    Opt@SrcGridLat = src_lat 
    Opt@SrcGridLon = src_lon 

;---destination grid information 
    Opt@DstRegional = True 
    Opt@DstGridLat = dst_lat 
    Opt@DstGridLon = dst_lon 
 
    Opt@ForceOverwrite = True 
    Opt@CopyVarCoords = True ; Whether to copy coordinate information 
                                 ; to regridded variable 
    Opt@InterpMethod   = "bilinear"   ; use "conserve" for precipitation
    
; ---- Regrid 
   var_regrid = ESMF_regrid(var,Opt) ; Regrid "var" to new grid 

   printVarSummary(var_regrid)
   printMinMax(var,0)
   printMinMax(var_regrid,0)

   var@lat2d = src_lat
   var@lon2d = src_lon

   wks = gsn_open_wks("png","wrf_to_rect_t2m")

   res                       = True

   res@gsnMaximize           = True

   res@gsnDraw               = False
   res@gsnFrame              = False

   res@cnFillOn              = True
   res@cnLinesOn             = False
   res@cnLineLabelsOn        = False

   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)
   res@mpDataBaseVersion     = "MediumRes"
   res@pmTickMarkDisplayMode = "Always"

;---Fix the contour levels across all time steps
   mnmxint = nice_mnmxintvl( min(var), max(var), 18, False)
   res@cnLevelSelectionMode = "ManualLevels"
   res@cnMinLevelValF       = mnmxint(0)
   res@cnMaxLevelValF       = mnmxint(1)
   res@cnLevelSpacingF      = mnmxint(2)

;--- Regional data, don't add longitude cyclic point
   res@gsnAddCyclic  = False

;---Title stuff
   res@tiMainFont     = "helvetica"
   res@pmTitleZone    = 4
   res@gsnRightString = ""
   res@gsnLeftString  = ""

;---Compare the plots in a panel
   pres                    = True
   pres@gsnMaximize        = True
   pres@gsnPanelLabelBar   = True
   pres@pmLabelBarWidthF   = 0.8
   pres@lbLabelFontHeightF = 0.015

   ntim = dimsizes(var(:,0,0))
   src_dims = str_join(""+dimsizes(var(0,:,:)),"x")
   dst_dims = str_join(""+dimsizes(var_regrid(0,:,:)),"x")

;---Loop across each time step and create a panel plot
   do nt=0,ntim-1,2    ; do every other time step for smaller GIF file.
     res@tiMainString  = "original data : " + src_dims
     plot_orig = gsn_csm_contour_map(wks,var(nt,:,:),res)  

     res@tiMainString  = "regridded data : " + dst_dims
     plot_regrid = gsn_csm_contour_map(wks,var_regrid(nt,:,:),res)

     pres@gsnPanelMainString = var@description + " (" + var@units + ") - " + wrf_times(nt)
     gsn_panel(wks,(/plot_orig,plot_regrid/),(/1,2/),pres)
  end do

;---Convert to animated GIF.
  system("convert -delay 25 wrf_to_rect*.png wrf_to_rect_t2m.gif")
end