;----------------------------------------------------------------------
; ESMF_regrid_template.ncl
;
; This is a template file for use with ESMF regridding. This is the
; template that uses the single "ESMF_regrid" function to do the 
; regridding. There's another template, "ESMF_template.ncl", that
; does the regridding in several individual steps. This can be 
; useful if you need to skip any of the steps.
;
; The ESMF_regrid function does the following:
;
;   1. Generates a description file (SCRIP or ESMF) for the source grid.
;
;   2. Generates a description file (SCRIP or ESMF) for the destination
;      grid.
;
;   3. Generates the weights file, using the source and destination
;      files.
;
;   4. Applies the weights to the data you want to regrid.
;
;   5. Attaches metadata to regridded variable, if appropriate.
;
;----------------------------------------------------------------------
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

   INTERP_METHOD = "bilinear" 
;   INTERP_METHOD = "neareststod" 

;---Input file that contains lat/lon grid of source grid, and some data to regrid
;    srcFileName = "../CESM_CISM/LGM_CLIMATE/HUY3_SCEAN1D/HUY3-SCEAN1D_0.nc"    
    srcFileName = "../../CESM_CISM/LGM_CLIMATE/HUY3_SCEAN1D/PMIP4_SETUP_CAM/Orog_USGS-huy3-S1D_21_10min.nc"    

;---Input file that contains lat/lon grid of destination grid, if applicable.
;--- from CESM input data
    dstFileName = "/Volumes/labdata/sarah/10min-30sec_CAM/gmted2010_modis-rawdata_htopo.nc"
    
;---Open the source data file and get some data to regrid
    setfileoption("nc","Format","NetCDF4Classic")
    sfile = addfile(srcFileName,"r")
	
    x     = sfile->htopo(:,:)            ; You may need to subscript this.
;    p     = sfile->landmask(0,:,:)            ; You may need to subscript this.

;---Open the destination data file
    setfileoption("nc","Format","NetCDF4Classic")
    dfile = addfile(dstFileName,"r")

;---Set up options for regridding
    Opt                = True

    Opt@SrcGridLat     = sfile->lat    ; Be sure to use appropriate names 
    Opt@SrcGridLon     = sfile->lon    ;here
    Opt@SrcRegional    = True         ; Necessary if source grid not global
    Opt@SrcTitle       = "NCEP Grid"   ; Optional 
;    Opt@SrcGridMask      = where(.not.ismissing(x),1,0) ; Necessary if has
                                                      ; missing values.
;    Opt@DstGridMask      = where((d_mask .eq. 1).or.(d_mask .eq. 4),1,0)

;
; If your destination grid is contained in a file, then use this code.
; Otherwise, if you just want to regrid to a generic grid, like a
; "1x1 degree grid, then skip to the next set of resources.
;
    dstlat             = dfile->lat    ; Be sure to use appropriate names 
    dstlon             = dfile->lon    ; here.
    Opt@DstGridLon     = dstlon
    Opt@DstGridLat     = dstlat
;    Opt@DstMask2D = where(.not.ismissing(dstlat).and.\
;                          .not.ismissing(dstlon),1,0) ; Necessary if lat/lon
                                                      ; has missing values.
;
; Use this if you just want to regrid to some generic grid, like a "1x
; "1x1 degree grid.
;
;;  Opt@DstGridType    = "1x"

    Opt@DstRegional    = True          ; Necessary if destination grid not global

    Opt@ForceOverwrite = True          ; Optional, but recommended
    Opt@PrintTimings   = True          ; Optional

    Opt@InterpMethod   = INTERP_METHOD    ; "patch", "conserve"

     x_regrid = ESMF_regrid(x,Opt)
;    p_regrid = ESMF_regrid(p,Opt)


;      setfileoption("nc","Format","LargeFile")				; need to produce nc files > 2 GB 
      rgrdFileName = "../CESM_CISM/LGM_CLIMATE/HUY3_SCEAN1D/PMIP4_SETUP_CAM/30sec/htopo_USGS-huy3-S1D_21_30sec.nc"      
;      rgrdFileName = "/Volumes/labdata/sarah/30sec_HUY3/huy3-S1D_21_alt.nc"      
      system("rm -f " + rgrdFileName)
      setfileoption("nc","Format","NetCDF4Classic")				; need to produce nc files > 2 GB 
      rgrd_nc = addfile(rgrdFileName,"c")

;---Create variable to hold global file attributes
      global = True
      copy_VarAtts(dfile, global)

      if (isatt(dfile,"title")) then
      global@TITLE = "REMAPPED: " + dfile@title
      end if
      global@remap         = "NCL: ESMF_regrid_with_weights (NCL version '" + \
                             get_ncl_version() + "')"
      global@remap_method  = INTERP_METHOD 
      global@creation_date = systemfunc("date")
    
      fileattdef(rgrd_nc, global)        ; copy global file attributes
      filedimdef(htopo,"lat",-1,True)   ; force an unlimited dimension
;
; Write variables to file. Coordinate arrays will be written 
; automatically
;

      rgrd_nc->htopo         = x_regrid
;     rgrd_nc->landmask    = p_regrid
       

end