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

; WRITE_NETCDF = False  
   
; TEXT FILE CONVERTING ROTATED AND SPHERICAL LATS AND LONS
filename = "hirham_rot_geo_koor.dat"
coords = asciiread(filename,(/100,110,4/),"float")
lat2d = coords(:,:,0) ; reading rotated latitude from ascii file
lon2d = coords(:,:,1) ; reading rotated longitude from ascii file
printMinMax(lat2d, True)
printMinMax(lon2d, True)

;---NCEP-NCAR reanalysis
    srcFileName = "slp.1948.nc"
    sfile = addfile(srcFileName,"r")

;---Read variable to be regridded
    SLP = sfile->slp 
    printVarSummary(SLP)   ; note this grid is rectilinear
 
;---Set the regridding options
    Opt                = True
    Opt@SrcTitle       = "NCEP-NCAR reanalysis"   ; optional

    Opt@WgtFileName    = "NCEP_to_Rect.nc"

    Opt@ForceOverwrite = True

    Opt@SrcFileName     = "Rectilinear.nc"      ; Name of source and
    Opt@DstFileName     = "WRF_SCRIP.nc"    ; destination files

    Opt@DstGridLat           = lat2d      ; curvilinear grid
    Opt@DstGridLon           = lon2d

    Opt@InterpMethod         = "bilinear"
    Opt@SrcRegional          = False
    Opt@DstRegional          = True
    Opt@Debug                = True

    SLP_regrid = ESMF_regrid(SLP,Opt)     ; Do the regridding for SLP

;
; The source and destination grid description files and
; weight file will be the same for the next call to 
; ESMF_grid, so no need to regenerate them.
;
    Opt@SkipSrcGrid   = True
    Opt@SkipDstGrid   = True
    Opt@SkipWgtGen    = True
 
;NOTE: I'M NOT SURE THIS PART IS NEEDED, SO I COMMENTED IT OUT.
;---Reset 0 values to missing values.
;;    SLP_regrid@_FillValue = default_fillvalue(typeof(SLP_regrid))
;;
;;    SLP_regrid            = where(SLP_regrid.eq.0.0,SLP_regrid@_FillValue,\
;;                                  SLP_regrid)
    printVarSummary(SLP_regrid)
    printMinMax(SLP_regrid,0)
    printMinMax(SLP,0)

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

    MAP_ZOOM = True     ; Whether to zoom in on area of rotated grid

    wks_name = "interpolate_slp"
    if(MAP_ZOOM) then
      wks_name = wks_name + "_zoom"
    end if
    wks_slp = gsn_open_wks("png",wks_name)

    res                       = True

    res@gsnDraw               = False
    res@gsnFrame              = False

    if(MAP_ZOOM) then
      res@mpMinLatF           = min(lat2d)
      res@mpMaxLatF           = max(lat2d)
      res@mpMinLonF           = min(lon2d)
      res@mpMaxLonF           = max(lon2d)
      res@pmTickMarkDisplayMode = "Always"   ; better map tickmarks
      res@pmTitleZone           = 4          ; move title closer to plot

      mnmxint_slp = nice_mnmxintvl( min(SLP_regrid), max(SLP_regrid), 18, False)
    else
      mnmxint_slp = nice_mnmxintvl( min(SLP), max(SLP), 18, False)
    end if

    res@cnFillOn              = True
    res@cnLinesOn             = False
    res@cnLineLabelsOn        = False
    res@lbLabelBarOn          = False
    res@cnLevelSelectionMode  = "ManualLevels"

    nrec      = 0
    dims_orig = tostring(dimsizes(SLP(nrec,:,:)))
    dims_rgrd = tostring(dimsizes(lat2d))


;---SLP
    res@cnMinLevelValF  = mnmxint_slp(0)
    res@cnMaxLevelValF  = mnmxint_slp(1)
    res@cnLevelSpacingF = mnmxint_slp(2)/2.   ; Create more levels
    res@tiMainFontHeightF = 0.02
    res@gsnStringFontHeightF = 0.015

;---SLP regridded
    res@tiMainString  = "curvilinear grid (" + str_join(dims_rgrd," x ") + \
                       ", " + Opt@InterpMethod + ")"
                        
    res@gsnAddCyclic          = False
    slp_regrid = gsn_csm_contour_map(wks_slp,SLP_regrid(nrec,:,:),res)

;---SLP original
    res@tiMainString = "SLP on original rectilinear grid (" +  \
                        str_join(dims_orig," x ") + ")"

    res@gsnAddCyclic          = True
    slp_orig = gsn_csm_contour_map(wks_slp,SLP(nrec,:,:),res)  

;---Compare the plots in a panel
    pres                   = True
    pres@gsnMaximize       = True
    pres@gsnPanelLabelBar  = True

    if(MAP_ZOOM) then
      gsn_panel(wks_slp,(/slp_orig,slp_regrid/),(/1,2/),pres)
    else
      gsn_panel(wks_slp,(/slp_orig,slp_regrid/),(/2,1/),pres)
    end if

;---Trim
    png_name = wks_name + ".png"
    system("convert -trim " + png_name + " " + png_name)

end