;	This script is written by ncl-talk@ucar.edu
;	It was edited to use "ESMF_regrid" to interpolate IMERG data (nc format) to station points.
;	Email from ncl-talk@ucar.edu:
;	As I said in an earlier email, you can make the writing of the ASCII files faster if you don't use a double do loop.  For example, this code:
;
;    do i=0, dimsizes(lat_NW)-1
;      do j=0, dimsizes(time)-1
;        write_table("GPM_R24_Interpolated_slow.txt","a",[/time(j),lat_NW(i),lon_NW(i),GPM_fod(j,i)/], \
;                     " %8i %7.4f %7.4f %5.1f")
;      end do
;    end do
;
;can be replaced with:
;
;    ntim = dimsizes(time)
;    nlat = dimsizes(lat_NW)
;    dims = (/nlat,ntim/)
;    time1d = ndtooned(conform_dims(dims,time,1))
;    lat1d  = ndtooned(conform_dims(dims,lat_NW,0))
;    lon1d  = ndtooned(conform_dims(dims,lon_NW,0))
;    write_table("GPM_R24_Interpolated_fast.txt","a",[/time1d,lat1d,lon1d,transpose(GPM_fod)/], \
;                 " %8i %7.4f %7.4f %5.1f")
;
;Note that because you are writing the time dimension out first and then the lat/lon dimension, I had to transpose the GPM_fod array before writing it.
;Also, the faster way does require the use of more memory, because I'm now creating extra arrays to hold all the time, lat, and lon arrays in a long 1D arrays.
;However, the speed up for writing the two files without do loops was significant, mainly for the first double do loop.

; load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"  ; This library is automatically loaded 
                                                                    ; from NCL V6.5.0 onward.
                                                                    ; No need for user to explicitly load.
	begin

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Part1: Reading CSV files(stations_NW.csv)
;https://www.ncl.ucar.edu/Applications/Scripts/csv_2.ncl
  print("Part1: Reading .csv files (stations_NW.csv)") ;(rbscReport_edited5_NW_R24h_Meteomanz.csv)
  fname_NW     = "stations_NW.csv"
  lines_NW     = asciiread(fname_NW,-1,"string")

  delim = ","

;---Read fields 1-5 from fname_Report
  station_id_NW      = tointeger(str_get_field(lines_NW(1:),1,delim))
  icao_code_NW       =           str_get_field(lines_NW(1:),2,delim)
  st_name_NW         =           str_get_field(lines_NW(1:),3,delim)
  lat_NW             =   tofloat(str_get_field(lines_NW(1:),4,delim))
  lon_NW             =   tofloat(str_get_field(lines_NW(1:),5,delim))
  
;---Print the information
  ;print("station_id_NW is '" + station_id_NW + "', icao_code_NW is " + icao_code_NW + ", st_name_NW is " + st_name_NW)
  print("End Part1 Reading .csv files")

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Part2: Interpolating GPM data at stations which have been got in previous section
;Example 1: https://www.ncl.ucar.edu/Document/Functions/Built-in/linint2_points.shtml
;Named subscripting: https://www.ncl.ucar.edu/Document/Manuals/Ref_Manual/NclVariables.shtml
;a sample to read NetCDF file: https://www.ncl.ucar.edu/Applications/netcdf4.shtml

  print("Part2: Interpolating GPM files")
  GPM_fnames = systemfunc("ls GPM_IMERGV04DE_20160401_20170228/3B-DAY-E.MS.MRG.3IMERG.*.nc")
;  print(GPM_fnames)
;	Following line should change in case of changing of input paths.
  GPM_fnamesc = str_get_cols(GPM_fnames,56,63)
  print(GPM_fnamesc)
  print("dimsizes(GPM_fnames) = "+dimsizes(GPM_fnames))
  print("dimsizes(lat_NW) = "+dimsizes(lat_NW))

  GPM_fod=new((/dimsizes(GPM_fnames),dimsizes(lat_NW)/),float)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;Creating file to write R24_GPM over grid points.
;Following 3 lines are commented for 2nd and later parts
  grid_output_file   = "GPMV04DER24Grid_ESMF_bilinear_20160401_20170228.txt"
  header_GPM_G		 =(/"   date     lat     lon   R24h_GPM_Grid"/)
  hlist_GPM_G		 =[/header_GPM_G/]
  write_table(grid_output_file,"w",hlist_GPM_G,"%s ")      ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  interp_method		  = "patch"   ; "neareststod"  ; "bilinear"   , "patch"
  interp_output_file      = "GPMV04DLR24Interpolated_ESMF_"+interp_method+"_20160401_20170228.txt"
  Opt                     = True
  Opt@DstGridType	  = "unstructured"
  Opt@SrcGridType	  = "rectilinear"
;  Opt@Overwrite       = True
  Opt@DstGridLat      = lat_NW    ; unstructured latitudes  
  Opt@DstGridLon      = lon_NW    ;      "       longitudes
  Opt@InterpMethod    = interp_method
; WgtFileName should be "nc"?
  Opt@WgtFileName     = interp_output_file

  Opt@ForceOverwrite  = True

  do_loop_start_time = get_cpu_time()

  do i=0, dimsizes(GPM_fnames)-1
 
;djs setfileoption("nc", "FileStructure", "Advanced")
    GPM_fi = addfile(GPM_fnames(i), "r")
  ;;setfileoption("nc", "Format",  "NetCDF4Classic")
    printVarSummary(GPM_fi);;;
    
    lat_GPM  = GPM_fi->lat
    lon_GPM  = GPM_fi->lon
    R24h_GPM = GPM_fi->precipitationCal
    ;print("  --dimsizes(R24h_GPM)= "+dimsizes(R24h_GPM))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ;;Writing R24h_GPM over grid points in the file
      start_time = get_cpu_time()
      nlat = dimsizes(lat_GPM)
      nlon = dimsizes(lon_GPM)
      dims = (/nlat,nlon/)
      gpm_fnamesc_tmp := ndtooned(conform_dims(dims,GPM_fnamesc(i),-1))
      lat1d           := ndtooned(conform_dims(dims,lat_GPM,0))
      lon1d           := ndtooned(conform_dims(dims,lon_GPM,1))
      write_table(grid_output_file,"a",[/gpm_fnamesc_tmp,lat1d,lon1d, transpose(R24h_GPM)/], " %s %7.4f %7.4f %5.1f")
      end_time = get_cpu_time()
      print("Elapsed time for writing GPM_R24_Grid file (FAST) = " + (end_time-start_time) + " CPU seconds")
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;    GPM_fo = linint2_points(lon_GPM,lat_GPM,R24h_GPM(lat | :, lon | :), False , lon_NW, lat_NW, 0)
;           = dspnt2(xi, yi, zi, xo1d, yo1d)   ; Interpolate using dspnt2.
;	GPM_fo = dspnt2(lon_GPM,lat_GPM,R24h_GPM(lat | :, lon | :) , lon_NW, lat_NW)
	Opt@SrcGridLat      = lat_GPM      ; rectilinear lat[*]
	Opt@SrcGridLon      = lon_GPM      ;      "      lon[*]
	GPM_fo = ESMF_regrid(R24h_GPM,Opt) ; Do the regridding
    printVarSummary(GPM_fo)
    printMinMax(GPM_fo,0)

    GPM_fod(i,:)=GPM_fo

  end do

   do_loop_end_time = get_cpu_time()
   print("Elapsed time for do loop = " + (do_loop_end_time-do_loop_start_time) + " CPU seconds")
  print("Finish Writing R24h_GPM over grid points in GPM_R24_Grid.txt")

  ;Generating dates 20160401-20170228
;	Following 2 lines should be changed in case of changing times.
  TIME = yyyymmdd_time(2016, 2017, "integer")
  time = TIME({20160401:20170228})    ; coordinate subscripting
 
  header_GPM=(/"   date     lat     lon   R24h_GPM"/)
  hlist_GPM=[/header_GPM/]
  write_table(interp_output_file,"w",hlist_GPM,"%s ")

  start_time = get_cpu_time()
    ntim = dimsizes(time)
    nlat = dimsizes(lat_NW)
    dims = (/nlat,ntim/)
    time1d := ndtooned(conform_dims(dims,time,1))
    lat1d  := ndtooned(conform_dims(dims,lat_NW,0))
    lon1d  := ndtooned(conform_dims(dims,lon_NW,0))
    write_table(interp_output_file,"a",[/time1d,lat1d,lon1d,transpose(GPM_fod)/], \
                 " %8i %7.4f %7.4f %5.1f")
    end_time = get_cpu_time()
    print("Elapsed time for writing GPM_R24_Interpolated file (FAST) = " + (end_time-start_time) + " CPU seconds")
	print("End Part2 Interpolating GPM files")

	end