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

;;;
; 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
 T2m_orig  = (f->air)

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

;;; 
; 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_orig(dy,:,:)

;;;
;  Hey, interpolate!
   T2m_interp  = rgrid2rcm(lat,lon,T2m_orig_2d,lat2d,lon2d,1)

;;;
; !!!!!!! 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
   end if
       
   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"
 out_nc = addfile(fileout,"c")
 
 out_nc->time  = time
 out_nc->lat   = lat2d
 out_nc->lon   = lon2d
 out_nc->air   = T2mem_interp_3D

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