; Job submission is via ~/CHINA/VPRM/sh_scripts/submit_regrid_job.sh

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"


;----------------------------------------------------------------------
; This function reads the necessary variables for regridding off
; the weights file, and returns them in a List object, so we can reuse 
; them later.  This is so we don't have to keep opening and closing 
; the weights file multiple times, which can be an expensive operation 
; if it's a large file.
;----------------------------------------------------------------------
function get_regrid_vars(wgts_file_name)
begin
  fid = addfile(wgts_file_name,"r")
  row = fid->row
  col = fid->col
  S   = fid->S
  row = row-1    ; dimensions are Fortran based
  col = col-1
  dst_grid_dims = fid->dst_grid_dims(::-1)  ; dimensions are Fortran ordered
  return([/S,row,col,dst_grid_dims/])
end

;----------------------------------------------------------------------
; This function calls the sparse_matrix function to directly
; apply weights for regridding the given variable "x". The necessary
; weights information is stored in a list called "thelist".
;
; "thelist" must contain these variables, in this order:
;   thelist[0] :  S             - the weights to apply
;   thelist[1] :  row           - the row indexes (0-based)
;   thelist[2] :  col           - the columns indexes (0-based)
;   thelist[3] :  dst_grid_dims - the output grid dimensions
;----------------------------------------------------------------------
function regrid_var(x,thelist,dst_grid_dims)
begin
  xregrid = tofloat(reshape(sparse_matrix_mult(thelist[1],thelist[2],thelist[0],\
                    ndtooned(x),product(thelist[3])),dst_grid_dims))

  return(xregrid)
end

begin

; ALL DIRECTORIES DEFINED BELOW **************************************************
    src_dir     = "./" 	;where your hourly avg files to regrid are located
    diro 	= "./" 	; where you want your netcdf regridded files to be put
    wgt_dir 	= "./" 	; where your weights file is located (previously generated 
			; by the weight generation script, bilinear interpolation.
; ALL DIRECTORIES DEFINED ABOVE *************************************************

;---Data file containing source grid; open as "nc" file.........................
    yr_inp = 2006 ;field set in shell script by sed for monthly regrids
    mo_inp = 1 ; field set in shell script by sed for monthly regrids
    search_term = yr_inp+"-"+sprinti("%0.2i", mo_inp)
    src_files 	= systemfunc("cd "+src_dir+";"+"ls aux*"+search_term+"*nc")
    sfiles    	= addfiles(src_dir+src_files,"r")

;---Specify variables to regrid..................................................
    var_names = (/"T2","SWDOWN"/)

;---Set up variable to hold regridding options...................................
    wgt_file_name       = wgt_dir+"T2SWDOWN_to_MODIS_bilinear.faster.nc"
    Opt                	= True

;---Set up templates.............................................................

var      = sfiles[0]->$var_names(0)$(:,:,0)
;;src_lat  = sfiles[0]->$"XLAT"$
;;src_lon  = sfiles[0]->$"XLONG"$
printVarSummary(var)

	;---Set bounds
;;	minlat = floor(min(src_lat))
;;	minlon = floor(min(src_lon))
;;	maxlat = ceil(max(src_lat))
;;	maxlon = ceil(max(src_lon))

	;---Get grid dimensions, template for building array
;
; Regrid one of the variables in order to get the lat/lon. There are
; better and more efficient ways to do this; this is the lazy way.
;
        rlist      = get_regrid_vars(wgt_file_name)
	var_regrid = regrid_var(var,rlist,(/6275,7692/))
        ESMF_copy_varmeta(var,var_regrid,wgt_file_name,False)
        lat  = var_regrid&lat
        lon  = var_regrid&lon
        delete(var_regrid)       ; no longer needed
	ZDim = ispan(0,23,1)     ; for assigning Z time coordinate
        nz   = dimsizes(ZDim)
        nlat = dimsizes(lat)
        nlon = dimsizes(lon)


; Loop through each file............................................................
do sfile = 0, dimsizes(src_files)-1
;
; Read variable and reorder. This is necessary, because the
; regridding expects lat/lon to be the rightmost two dimensions.

  T2_llt = sfiles[sfile]->$var_names(0)$
  T2     = T2_llt(ncl3|:,south_north|:,west_east|:)
  delete(T2_llt)

  SWDOWN_llt = sfiles[sfile]->$var_names(1)$
  SWDOWN     = SWDOWN_llt(ncl4|:,south_north|:,west_east|:)
  delete(SWDOWN_llt)

; Do the regridding.............................................................
  print("Regridding T2...")
  T2_array = (/regrid_var(T2,rlist,(/nz,nlat,nlon/))/)

;; Print out information about original and regridded T2
  printVarSummary(T2_array)
  printMinMax(T2,0)
  printMinMax(T2_array,0)
  print("Regridding SWDOWN...")
  SWDOWN_array = (/regrid_var(SWDOWN,rlist,(/nz,nlat,nlon/))/)

;; Print out information about original and regridded SWDOWN
  printVarSummary(SWDOWN_array)
  printMinMax(SWDOWN,0)
  printMinMax(SWDOWN_array,0)

        ;----------------------------------------------------------------------
        ; WRITE DAILY REGRIDDED FILES (HRLY AVERAGES) OUT TO NETCDF
        ; This section writes the regridded files to netcdf.
        ;----------------------------------------------------------------------

        filo     = diro+"REGRID_"+src_files(sfile)
        ; SET NETCDF FILE PROPERTIES...........................................
        system("/bin/rm -f " +filo)
        setfileoption("nc","Format","NetCDF4")
	setfileoption("nc","CompressionLevel",4)
        ; CREATE NETCDF FILE WITH ABOVE PROPERTIES ............................       
        ncdf        = addfile(filo ,"c")
        
        ; make time and UNLIMITED dimension. 
        ;; This doesn't really make sense here, because you don't have a variable
        ;; with a "time" dimension. You have one with a "ZDim" dimension, which
        ;; I guess is supposed to be time?
        ; (This is recommended for most applications)
;;        filedimdef(ncdf,"time",-1,True)    

        ; fill in missing attributes ..........................................

;;---Attach coordinate arrays and attributes only once (first time in loop)
        if(sfile.eq.0) then
  	    T2_array!0    = "ZDim"    ; should this be "time"?
	    T2_array!1    = "lat"
	    T2_array!2    = "lon"
	    T2_array&ZDim = ZDim      ; T2_array&time = ZDim???
	    T2_array&lat  = lat
	    T2_array&lon  = lon
	    T2_array@projection  = "Geographic"
	    T2_array@description = T2@description
	    ;T2_array@_FillValue = -9999
	    T2_array@units = T2@units

            SWDOWN_array!0    = "ZDim"    ; should this be "time"?
            SWDOWN_array!1    = "lat"
            SWDOWN_array!2    = "lon"
            SWDOWN_array&ZDim = ZDim      ; SWDOWN_array&time = Zdim?
            SWDOWN_array&lat = lat
            SWDOWN_array&lon = lon
	    SWDOWN_array@projection = "Geographic"
            SWDOWN_array@description = SWDOWN@description
           ;SWDOWN_array@_FillValue = -9999
            SWDOWN_array@units = SWDOWN@units
        end if
        ncdf->T2	       	= T2_array       ; this will write all coordinate arrays and attributes
        ncdf->SWDOWN	     	= SWDOWN_array

end do ; ends sfile

end ; ends script