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

data_path  = "/Users/linx/Data/NSIDC/NOAA/G02202_v2/south/monthly/"

year1      = 1979
year2      = 2015
ny         = year2-year1+1
year       = fspan(1979.,2013.999,ny*12)

method       = "bilinear" ; "neareststod"   ;"patch"  ;"bilinear"

; read the original nc files
y   = new((/ny*12,332,316/),"float")
index = "n07"

do i = 0, year2-year1
do j = 1, 12

   yymm  = ((i+1979)*100+j)
   if(yymm.eq.198707)then
      index = "f08"
   end if
   if(yymm.eq.199112)then
      index = "f11"
   end if
   if(yymm.eq.199510)then 
      index = "f13"
   end if
   if(yymm.eq.200801)then
      index = "f17"
   end if
   data_filename = "seaice_conc_monthly_sh_" + index + "_" + ((i+1979)*100+j) + "_v02r00.nc"
   print(" " + data_filename)

   fi     = addfile(data_path + data_filename,"r")
   lon    = fi->longitude
   lat    = fi->latitude    

   data   = fi->goddard_merged_seaice_conc_monthly    ; byte type in default
   x      = data(0,:,:)*data@scale_factor                           

   ;flag_values :  ( -5, -4, -3, -2, -1 )
   ;flag_meanings : pole_hole unused coastal land_mask missing_data  
   n1   = num(data.eq.-1)
   n2   = num(data.eq.-2)
   n3   = num(data.eq.-3)
   n4   = num(data.eq.-4)
   n5   = num(data.eq.-5)
   print(n1+" " + n2 + " " + n3 + " " + n4 + " " + n5)

   x      = where(data(0,:,:).lt.0, -999., x)
   x@_FillValue = -999.

   y(i*12+j-1,:,:)   = x

   nmiss   = num(ismissing(x))
   n0      = num(x.eq.0)
   ntotal  = num(x.gt.0)
   print(i+1979 + "  " + j + "  missing "  + nmiss + "  zero values " + n0 + "  positive " + ntotal)
 
end do
end do
  
  y@lat2d  = lat
  ;lon      = where(lon.lt.0,lon+360,lon)     ; -180 ~ 180 longitude order
  ;lon@units= "degrees_east"
  y@lon2d  = lon                     

   print(min(lon) + "   " + max(lon))
;  Loop over , read over

;  read the destination grid file (HadISST)
     fi   = addfile("/Users/linx/Data/Hadley/HadISST.2.2.0.0_sea_ice_concentration.nc","r")
     hsic = fi->sic(:,{-90.:-30.},:)
     hlon = fi->longitude                   ; 0 ~ 360
     hlat = fi->latitude({-90.:-30.})
   

; ************* regridding *********************.
    Opt                  = True
    Opt@SrcFileName      = "source_grid_file.nc"       ; source file name
    Opt@DstFileName      = "destination_grid_file.nc"   ; destination file name
    Opt@ForceOverwrite   = True

    ;Opt@SrcGridMask      = where(.not.ismissing(y(0,:,:)),1,0)
    ;Opt@DstGridMask      = where(.not.ismissing(hsic(0,:,:)),1,0)

    Opt@DstGridLat       = hlat     
    Opt@DstGridLon       = hlon
    printVarSummary(hlon)

      print("Generating interpolation weights from CMIP5 to")
      print("World 1 degree grid using the " + method + " method.")

      Opt@WgtFileName  = "ESMF_regrid_wgt_nsidc.sic.nc" 
      Opt@InterpMethod = method
      
      printVarSummary(y)
      thetao_regrid = ESMF_regrid(y,Opt)
      printVarSummary(thetao_regrid)

      time = flt2dble(fspan(0., ny*12-1., ny*12) )   ;185001-201212:1956
      time@units      = "months since 1979-01-01 00:00"
      time@long_name  = "Time"
      thetao_regrid!0 = "time"
      thetao_regrid&time  = time

; ******* output netcdf file ******************
       system("rm -f " + "nsidc_bt.sic_regrid.nc")   ; remove any pre-existing file
       ncdf = addfile("nsidc_bt.sic_regrid.nc" ,"c")  ; open output netCDF file

    ;===================================================================
    ; create global attributes of the file (optional)
    ;===================================================================
       fAtt               = True            ; assign file attributes
       fAtt@title         = "NSIDC Bootstrap sea ice concentration"
       ;fAtt@source_file   =  "original-file.nc"
       ;fAtt@Conventions   = "None"
       fAtt@creation_date = systemfunc ("date")
       fileattdef( ncdf, fAtt )            ; copy file attributes

    ;===================================================================
    ; make time an UNLIMITED dimension; recommended  for most applications
    ;===================================================================
       filedimdef(ncdf,"time",-1,True) 

    ;===================================================================
    ; output variables directly; NCL will call appropriate functions
    ; to write the meta data associated with each variable
    ;===================================================================
       ncdf->sic  = thetao_regrid
       ncdf->time = time

      delete(thetao_regrid)
                  

end