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

;*********************************************************************
; Source Grid (ASR file) 
;*********************************************************************
fnlDir   = "/home/vanucia/"              
fnlDir   = "./"
fnlFile  = "shf_asr-2000-classic.nc"      
f        = addfile(fnlDir + fnlFile, "r") 
varSrc0 = f->LH 
varnames = (/"LH"/)

;**********************************************************************
;The flux data does not have XLAT and XLONG in it. I have included the 
;this file of the static data. Be sure it is in your directory with this
; script. 
;**********************************************************************
     g   = addfile("asr15km.static.2d.2007070100.nc","r")
;latSrc   = g->XLAT(::6,::6)    ; Use these if you have smoothed ASR to 
;lonSrc   = g->XLONG(::6,::6)   ; every 6 grid points and are now regridding
latSrc   = g->XLAT
lonSrc   = g->XLONG

printVarSummary(latSrc)
printVarSummary(lonSrc)
printVarSummary(varSrc0)


;***********************************************************************
; Destination Grid (Global Grid you are regridding to). Be sure to change
; the names of the latitude and longitude that are in your 1x1 deg 
; global grid file. 
;***********************************************************************
wrfDir   = "/home/vanucia/"; directory of global 
wrfDir   = "./"; directory of global 
wrfFile  = "shf_era-1x1.nc" ; filename
h        = addfile(wrfDir + wrfFile, "r") 
latDst   = h->lat
lonDst   = h->lon

;************************************************************************
; Options to pass to ESMF_regrid <-Should leave these as is. 
; Your sources is ASR (Opt@SrcRegional = True) and your Destination is
; global Opt@DstRegional    = False
;************************************************************************
Opt                = True
Opt@SrcGridType = "curvilinear" 
Opt@DstGridType = "rectilinear" 
Opt@SrcGridLat     = latSrc
Opt@SrcGridLon     = lonSrc
Opt@DstGridLat     = latDst
Opt@DstGridLon     = lonDst  
Opt@SrcRegional    = True  ; because the source grid (ASR)  is a regional grid
Opt@DstRegional    = False ; if regridding to global; set to True if smoothing ASR
Opt@InterpMethod   = "bilinear" 
Opt@NoPETLog       = True

;Opt@RemoveSrcFile  = True
;Opt@RemoveDstFile  = True
;Opt@RemoveWgtFile  = True

;Opt@SkipSrcGrid  = True     ; Once you have the weights file,
;Opt@SkipDstGrid  = True     ; you can set these to True.
;Opt@SkipWgtGen   = True

Opt@Debug          = True 
Opt@CopyVarCoords  = True   ; True is the default in later versions of NCL
                            ; if True, the lat-lon data will be attached as
Opt@ForceOverwrite = True


varDst0 = ESMF_regrid(varSrc0,Opt)
printVarSummary(varDst0)
printMinMax(varSrc0,0)       
printMinMax(varDst0,0)       

; Write to netcdf

WRITE_NETCDF = False
if(WRITE_NETCDF) then
  fout = "test.nc"    
  system("rm -f " + fout)
  ncdf = addfile(fout, "c")
  ncdf->$varnames(0)$ = varDst0
  print("")
  print("Regridded data saved to " + fout)
end if

PLOT_DATA = True

if(.not.PLOT_DATA) then
  exit
end if

;----------------------------------------------------------------------
; Plot the data for extra assurance that everything looks good.
;----------------------------------------------------------------------
varSrc0@lat2d = latSrc     ; Needed for plotting.
varSrc0@lon2d = lonSrc

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

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

;
; Source grid covers smaller area than destination grid,
; so use source grid limits for the plotting. You can
; comment out these four lines to get the full global
; lat/lon area.
;
res@mpMinLatF             = min(latSrc)
res@mpMaxLatF             = max(latSrc)
res@mpMinLonF             = min(lonSrc)
res@mpMaxLonF             = max(lonSrc)

mnmxint = nice_mnmxintvl( min(varDst0), max(varDst0), 18, False)
res@cnLevelSelectionMode = "ManualLevels"
res@cnMinLevelValF       = mnmxint(0)
res@cnMaxLevelValF       = mnmxint(1)
res@cnLevelSpacingF      = mnmxint(2)

nt = 0   ; time step to plot

;---Resources for plotting regridded data
dims = tostring(dimsizes(varDst0(nt,:,:)))

res@gsnAddCyclic  = True
res@tiMainString  = "Regridded to rectilinear grid using " + Opt@InterpMethod + \
                     "method (" +  str_join(dims," x ") + ")"

plot_regrid = gsn_csm_contour_map(wks,varDst0(nt,:,:),res)

;---Resources for plotting original data

dims = tostring(dimsizes(varSrc0(nt,:,:)))
res@gsnAddCyclic = False     ; original data is regional
res@tiMainString = "Original curvilinear grid (" + \
                    str_join(dims," x ") + ")"

plot_orig = gsn_csm_contour_map(wks,varSrc0(nt,:,:),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