;::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
;
;    This file is part of ICTP RegCM.
;
;    ICTP RegCM is free software: you can redistribute it and/or modify
;    it under the terms of the GNU General Public License as published by
;    the Free Software Foundation, either version 3 of the License, or
;    (at your option) any later version.
;
;    ICTP RegCM is distributed in the hope that it will be useful,
;    but WITHOUT ANY WARRANTY; without even the implied warranty of
;    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;    GNU General Public License for more details.
;
;    You should have received a copy of the GNU General Public License
;    along with ICTP RegCM.  If not, see <http://www.gnu.org/licenses/>.
;
;::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

;**************************************************
;************ Load necessary ncl code *************
;**************************************************

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/csm/shea_util.ncl"

begin

  ;User-defined vars
  ;Open the regcm header, where we'll grab lat2d and lon2d
  fho    =    addfile("/home/samy/RegCM-4.4.5.3_run/input/AFRICA_DOMAIN000.nc","r")
  ;Open the CRU Temperature data file
  f1    =    addfile("/home/samy/RegCM-4.4.5.3_run/output/clm_tg_2.nc","r")
  ;Open the CRU Precip data file
  ;f2    =    addfile(cruprefilename,"r")

  ;Grab lat2d and lon2d
  lat2df        =    fho->xlat(:,:) 
  lon2df        =    fho->xlon(:,:)

  ;Get the number of latitudes and longitudes
  dvar        =    dimsizes(lat2df)
  nlatf        =    dvar(0)
  nlonf        =    dvar(1)

  ;remove the border points, since the RCM output files
  ;do the same
  lat2d = lat2df(1:nlatf-3,1:nlonf-3)
  lon2d = lon2df(1:nlatf-3,1:nlonf-3)

  nlat        =    nlatf-3
  nlon        =    nlonf-3

  ;read the CRU lat, lon, and time
  lat_cru     = f1->lat
  lon_cru     = f1->lon
  time_cru    = f1->time

  ;Find the start and end times
  ;uttime = ut_calendar(time_cru,-1)
  istart = 0
  ;ind(uttime.eq.istartdate)
  iend = 4
   ;ind(uttime.eq.ienddate)
  ;if(ismissing(istart).or.ismissing(iend))then
   ; print("Error: couldn't find corresponding start and end dates in CRU file.  Abort.")
    ;return
  ;end if

   ntime = dimsizes(time_cru(istart:iend))

  print("Reading groundtemp from " + time_cru(istart) + " to " + time_cru(iend))
  ;Read in the CRU data and give it proper units
  v1_tg       =   f1->tg(istart:iend,:,:)
  ;v1_tg        =  doubletofloat(f1->tg(istart:iend,:,:))
  v1_tg        =    v1_tg + 273.15
  v1_tg@units    = "K"
  
  ;Remove the CRU file sitting in this directory
  sCRUFileName ="cru." + istart + "_" + iend + ".nc"
  system("/bin/rm -f "+ sCRUFileName)
  ;Open the CRU file for outputting the re-gridded CRU data
  fout    =    addfile(sCRUFileName,"c")

  ;Regrid the CRU temperature to the RCM grid
  print("Re-gridding tg")
  ;printVarSummary(tmp) 
  tg_regrid    = rgrid2rcm(lat_cru,lon_cru,v1_tg,lat2df,lon2df,1)
  tg = new((/ntime,1,nlat,nlon/),float)
  printVarSummary(tg)
  printVarSummary(tg_regrid)
  tg(:,:,:) = tg_regrid
  ;printVarSummary(tmp)
  ;(:,0,:,:) = ta_regrid
         
  ;Give ta the RCM coordinates
  tg!0 = "time"
  ;tmp!1 = "K"
  ;tg!1 = "iy"
  tg!1 = "iy"
  tg!2 = "jx"
  ;tg&time = v1_tg&time
  ;ta&m2 = 2
  tg&iy = lat2d&iy
  tg&jx = lat2d&jx
  tg@_FillValue = v1_tg@_FillValue
  tg@_FillValue = 1e-30
  tg@missing_value = tg@_FillValue
  ;printVarSummary(clt)
  ;Write out ta
  print("Writing tg")
  fout->tg = tg

  ;Read in the CRU precip data
  ;print("Reading pr (TPR) from " + time_cru(istart) + " to " + time_cru(iend))
  ;v2_PRE        =    short2flt(f2->PRE(istart:iend,:,:))

  ;Regrid the CRU precip data to the RCM grid
  ;print("Re-gridding pr")
  ;pr_regrid    = rgrid2rcm(lat_cru,lon_cru,v2_PRE,lat2d,lon2d,1)
  ;pr = new((/ntime,nlat,nlon/),float)
  ;pr(:,:,:) = pr_regrid
  ;Give it the RCM coordinates
  ;pr!0 = "time"
  ;pr!1 = "iy"
  ;pr!2 = "jx"
  ;pr&time = v2_PRE&time
  ;pr&iy = lat2d&iy
  ;pr&jx = lat2d&jx
  ;pr@_FillValue = v2_PRE@_FillValue
  ;pr@_FillValue = 1e-30
  ;pr@missing_value = pr@_FillValue

  ;Output the precip data
  ;print("Writing pr")
  ;fout->pr = pr(:,:,:)

end