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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"

begin

;;;
; First, grab the latitude and longitude to interpolate FROM)
 f = addfile("air.1950.nc","r")
 lat=f->lat
 lon=f->lon
 yr = "1950"

;;;
; for the function rgrid2rcm, the latitudes need to be monotonically increasing, so do that here
 lat=lat(::-1)

;;;
; Now grab the latitude/longitude to interpolate TO
 fwrf = addfile("geo_em.d01.nc","r")

 lat2d=fwrf->XLAT_M(0,:,:)
 lon2d=fwrf->XLONG_M(0,:,:)
 dims=dimsizes(lat2d)

;;;
; Make all longitudes positive for interpolating purposes
 lon2d = where(lon2d.lt.0,lon2d+360.,lon2d)

;;; 
; get data for interpolating
 dims1d = dimsizes(f->air)
 print(dims1d)
 dims1d(2) = dims1d(2) + 1 ; increase the longitude dim
 T2m_orig  = f->air
 T2m_mod = new(dims1d,typeof(T2m_orig))
 print(dims1d) 
 T2m_mod(:,:,0:dims1d(2)-2) = T2m_orig  ; fill all except the new last longitude column with the original array
 T2m_mod(:,:,dims1d(2)-1) = T2m_mod(:,:,0) ; fill the last longitude column with the values from the first longitude column
 T2m_mod&lon(dims1d(2)-1) = lon(0) + 360  ; fill in the last longitude coodinate value with the first value (0) + 360 (maintaining monotonicity)
 printVarSummary(T2m_mod)
 lonnew = T2m_mod&lon
 
 

;;;
; Switch the order of the latitudes, since I had to switch the order of latitudes above
 T2m_orig  = T2m_orig(:,::-1,:)
 printVarSummary(T2m_orig)

;;; 
; Grab time variable which has how many days in each year
 time=f ->time
 dimtimes = dimsizes(time)

 T2mem_interp_3D = new((/dimtimes,dims(0),dims(1)/),float)

;;;
; copy attributes
 T2mem_interp_3D!0 = "time"
 T2mem_interp_3D!1 = "south_north"
 T2mem_interp_3D!2 = "west_east"

 do dy=0,dimtimes-1

   if( mod(dy, 10).eq.0 ) then
     print("Interpolating day "+(dy+1))
   end if
     
;;;
; Grab data for just this day
   T2m_orig_2d  = T2m_mod(dy,:,:)

;;;
;  Hey, interpolate!
   T2m_interp  = rgrid2rcm(lat,lonnew,T2m_orig_2d,lat2d,lon2d,1)
;   print(num(ismissing(T2m_interp)))
;   print(num(ismissing(lat2d)))
;   print(num(ismissing(lon2d)))
;   printMinMax(lon2d,0)
;   printMinMax(lat2d,0)
;  printMinMax(lon,0)
;   printMinMax(lat,0)
;  print(num(ismissing(T2m_orig_2d)))

  wks = gsn_open_wks("png","test-interp")
   res = True
  res@gsnMaximize           = True
  res@tiMainString          = "Lat/Lon interpolated to WRF grid"
  res@tiMainFont            = 22

;  res@cnFillColors          = ispan(2,200,10)
  res@cnFillPalette   =  "rainbow"
  res@cnFillOn              = True
  res@cnLineThicknessF      = 2.0
 ; res@cnLinesOn = False
  res@pmLabelBarDisplayMode = "Always"
  res@lbPerimOn             = False
;  plot = gsn_contour(wks,T2m_orig_2d,res)
;  plot = gsn_contour(wks,lat2d,res)
;  plot = gsn_contour(wks,lon2d,res)
;  plot = gsn_contour(wks,T2m_interp,res)
  T2m_interp@lat2d = lat2d
  T2m_interp@lon2d = lon2d
  printVarSummary(T2m_interp)
  ;poisson_grid_fill(T2m_interp,Flase
  res@mpProjection = "orthographic"
  res@mpLimitMode = "latlon"
  res@mpMinLatF = 58
  ;res@trGridType = "triangularmesh"
  res@mpCenterLatF = 90
  res@mpCenterLonF = -45 
  res@mpGridAndLimbOn = True
  plot = gsn_csm_contour_map(wks,T2m_interp,res)
 
  print("finished plots")
  
;;;
; !!!!!!! this is where the error happens !!!!!!!!!!
   if( any(ismissing(T2m_interp)) ) then
     print("There was missing data after interpolation on day "+dy+", exiting, 1")
     exit
   else   
     print("There are " + num(ismissing(T2m_interp)) + " missing values")
   end if
   exit
       
   copy_VarAtts(T2m_orig_2d  , T2m_interp)
   T2m_interp!0 = "south_north"
   T2m_interp!1 = "west_east"

   T2mem_interp_3D(dy,:,:) = T2m_interp


   delete(T2m_orig_2d)
   delete(T2m_interp)
 end do
    
 fileout = yr+".Greenland_domain.NCEP.nc"
 system("rm -f " + fileout)
 out_nc = addfile(fileout,"c")
 
 out_nc->time  = time
 out_nc->lat   = lat2d
 out_nc->lon   = lon2d
 print (out_nc)
 printVarSummary(T2mem_interp_3D)
 out_nc->air   = T2mem_interp_3D

 delete( [/T2mem_interp_3D, time/] )
 delete( [/T2m_orig, dimtimes/] )
 print("********************************************************")
  
end