;======================================================================
; ESMF_regrid_6.ncl
;
; Concepts illustrated:
;   - Interpolating from one grid to another using ESMF_regrid
;   - Interpolating data from a CMIP5 grid to a 1X1 degree rectilinear grid
;======================================================================
; This example is identical to ESMF_all_6.ncl, except it does the
; regridding in one call to "ESMF_regrid".  See ESMF_wgts_6.ncl
; for a faster example of regridding using an existing weights file.
;======================================================================
; This example uses the ESMF application "ESMF_RegridWeightGen" to 
; generate the weights.
;
; For more information about ESMF:
;
;        http://www.earthsystemmodeling.org/
;======================================================================
; This script regrids a CMIP5 grid to a 1.0 degree world grid and
; plots sea water potential temperature on the new grid.
;
; It uses SCRIP for both the CMIP5 and 1.0 degree world grid.
;======================================================================
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

     fi   = addfile("/data1/linx/Hadley/HadISST_ice.nc","r")
     hsic = fi->sic
     hlon = fi->longitude
     hlat = fi->latitude


;---Interpolation methods
    ;methods      = (/"bilinear","patch","conserve"/)
     method       = "bilinear"

;---Input file 
    dir_name = systemfunc("ls -F /data1/linx/CMIP5_SIC/historical | grep '/' ")
    nd       = dimsizes(dir_name)
    print(dir_name)
 
do i = 4,5
          path        = "/data1/linx/CMIP5_SIC/historical/" + dir_name(i) 
          srcFileName = systemfunc("ls " + path + "sic*.nc")
          print(i + ".     " + srcFileName)
          wgtFile     = "ESMF_regrid_wgt_" + method + ".nc"

;---Get data and lat/lon grid from CMIP5 Grid
    sfile        = addfile(srcFileName,"r")
    thetao       = sfile->sic
    thetao@lat2d = sfile->lat
    thetao@lon2d = sfile->lon

    ncor  = dimsizes(dimsizes(thetao@lat2d))

    Opt                  = True
    Opt@SrcFileName      = "source_grid_file.nc"       ; source file name
    Opt@DstFileName      = "destination_grid_file.nc"   ; destination file name
    Opt@ForceOverwrite   = True

    Opt@SrcGridMask      = where(.not.ismissing(thetao(0,:,:)),1,0)
    Opt@DstGridMask      = where(.not.ismissing(hsic(0,:,:)),1,0)

    ;if(ncor.eq.1)then
    ;   Opt@SrcGridType   = "rectilinear"
    ;else
    ;:end if

   ; Opt@DstGridType      = "1x1"              ; Destination grid
   ; Opt@DstTitle         = "World Grid 1-degree Resolution"
   ; Opt@DstLLCorner      = (/-89.5d, -179.5d/)    ;     (/-89.75d,   0.00d /)
   ; Opt@DstURCorner      = (/ 89.5d,  179.5d/)    ;     (/ 89.75d, 359.75d /) 

    Opt@DstGridLat       = hlat     
    Opt@DstGridLon       = hlon

    Opt@PrintTimings   = True
    Opt@Debug          = True

    print(Opt)
   
;----------------------------------------------------------------------
; Setup for graphics

;----------------------------------------------------------------------
    ch   = stringtochar(dir_name(i))
    nc   = dimsizes(ch)   
    dname= chartostring(ch(0:nc-3))
    wks  = gsn_open_wks("png", dname + "_ESMF_regrid")
    gsn_define_colormap(wks,"rainbow")             ; Change color map

;---Resources to share between both plots
    res                     = True              ; Plot modes desired.

    res@gsnDraw             = False             ; Will panel later
    res@gsnFrame            = False             ; Will panel later

    res@gsnMaximize         = True              ; Maximize plot

    res@cnFillOn            = True              ; color plot desired
    res@cnLinesOn           = False             ; turn off contour lines
    res@cnLineLabelsOn      = False             ; turn off contour labels
    res@cnFillMode          = "RasterFill"      ; turn raster on      
  
    ;res@cnLevelSelectionMode = "ExplicitLevels"
    ;res@cnLevels             = ispan(270,300,2)
  
    res@mpFillOn            = False

    res@trGridType         = "TriangularMesh"  ; allow missing coordinates
    res@gsnAddCyclic       = False

    res@pmLabelBarWidthF   = 0.7
    res@lbLabelBrOn       = False   ; Will do this in panel

    res@gsnAddCyclic       = False

;---Resources for paneling
    pres                  = True
    pres@gsnMaximize      = True
    pres@gsnPanelLabelBar = True
    pres@lbLabelFontHeightF = 0.01

;----------------------------------------------------------------------
; Loop across each method and generate interpolation weights for 
; CMIP5 Grid to World Grid    
;----------------------------------------------------------------------
    ; plot_regrid = new(dimsizes(methods),graphic)

    
      print("Generating interpolation weights from CMIP5 to")
      print("World 1 degree grid using the " + method + " method.")

      Opt@WgtFileName  = wgtFile
      Opt@InterpMethod = method

;----------------------------------------------------------------------
; Interpolate data from CMIP5 to World 1-degree grid.
;----------------------------------------------------------------------
    
      thetao_regrid = ESMF_regrid(thetao,Opt)
     printVarSummary(thetao_regrid)

;----------------------------------------------------------------------
; Plotting section
;----------------------------------------------------------------------

;---Resources for plotting original data
      res@tiMainString = "Data on original CMIP5 grid (" + \
                         str_join(tostring(dimsizes(thetao(0,:,:)))," x ") + ")"
      plot_orig        = gsn_csm_contour_map(wks,thetao(0,:,:),res)

;---Resources for plotting regridded data
      res@tiMainString = "CMIP5 to 1x1-degree grid (" + \
                         method + ") (" + \
                         str_join(tostring(dimsizes(thetao_regrid(0,:,:)))," x ") + ")"
      plot_regrid      = gsn_csm_contour_map(wks,thetao_regrid(0,:,:),res)

;---Panel two plots
      gsn_panel(wks,(/plot_orig,plot_regrid/),(/2,1/),pres)

;---Clean up before next time in loop.
      delete(thetao)
      delete(Opt)      
      delete(ch)
      delete(nc)
  
;    output netcdf file

     
       system("rm -f " + dname + "_sic_regrid.nc")   ; remove any pre-existing file
       ncdf = addfile(dname + "_sic_regrid.nc" ,"c")  ; open output netCDF file

    ;===================================================================
    ; create global attributes of the file (optional)
    ;===================================================================
       fAtt               = True            ; assign file attributes
       fAtt@title         = "CMIP5 SIC data to HadISST Grid: " + dname
       ;fAtt@source_file   =  "original-file.nc"
       ;fAtt@Conventions   = "None"
       fAtt@creation_date = systemfunc ("date")
       fileattdef( ncdf, fAtt )            ; copy file attributes

    ;===================================================================
    ; make time an UNLIMITED dimension; recommended  for most applications
    ;===================================================================
       filedimdef(ncdf,"time",-1,True) 

    ;===================================================================
    ; output variables directly; NCL will call appropriate functions
    ; to write the meta data associated with each variable
    ;===================================================================
       ncdf->sic = thetao_regrid
 
      delete(thetao_regrid)
      delete(dname)
            

end do


end