; SMOPS: soil Moisture Operational Products System 
; https://www.ospo.noaa.gov/Products/land/smops/
; Values over land only
;
; This file still has to be loaded manually
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"  ; not needed for 6.5.0
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Part1: Reading .dat file(2_stations.dat)

; print("Part1: Reading .dat files (2_stations.dat)")
  
  drst    = "./"
  flst    = "2_stations.dat"
  ptst    = drst + flst
  lnst    = asciiread(ptst,-1,"string")
  delim   = ","

;---Read fields 1-5 from fname_Report
  stid    = tointeger(str_get_field(lnst(1:),1,delim))
  stcd    = str_get_field(lnst(1:),2,delim)
  stnm    = str_get_field(lnst(1:),3,delim)
  stlon   = tofloat(str_get_field(lnst(1:),4,delim))
  stlat   = tofloat(str_get_field(lnst(1:),5,delim))
  
;---Print the information
  print("stid is '" + stid + "', stnm is " + stnm +": lon/lat: "+stlon+"     "+stlat)
  print("---")

  latS    = min(stlat)
  latN    = max(stlat)
  lonL    = min(stlon)
  lonR    = max(stlon)
  print("latS="+latS+"  latN="+latN+"   lonL="+lonL+"  "+lonR)
  print("---")

; print("End Part1 Reading 2_stations.dat files")

  nsta    = dimsizes(stid)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Creating Coordinates for SMOPS data
; They are not in the netCDF file
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  nlat			    = 720
  latitude		    = latGlobeFo(nlat, "lat", "latitude", "degrees_north")
  latitude = latitude(::-1)                           ; data are North-to-South
  printVarSummary(latitude)
  printMinMax(latitude,0)
	
  mlon	                    = 1440
  longitude		    = lonGlobeFo(mlon, "longitude", "longitude", "degrees_east")  ; 0->360
  longitude		    = (/ longitude - 180. /)  ; subtract 180 from all values 
  longitude&longitude	    = longitude               ; update coordinates       printVarSummary(lon)
  printVarSummary(longitude)
  printMinMax(longitude,0)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; poisson_grid_fill parameters
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  nscan = 2000             ; usually *much* fewer
  eps   = 0.001            ; variable depended
  gtype = True             ; cyclic
  guess = 0                ; use zonal means
  relc  = 0.6              ; standard relaxation coef
  iopt  = 0

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Read/open NCL'ls crude (1x1) land-sea mask
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  a    = addfile("$NCARG_ROOT/lib/ncarg/data/cdf/landsea.nc","r")
  lsdata = a->LSMASK
  lsm  = landsea_mask(lsdata, latitude, longitude)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Loop and Open each file(s)
; latitude/longitude are not in the netCDF file: must be added
; non-CF convention FillValue must be added: _FillValue
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  SMOPS_dir                 = "./"
  SMOPS_fnames              = systemfunc("cd "+SMOPS_dir+ " ; ls NPR_SMOPS_CMAP_*.nc")
  nSMOPS                    = dimsizes(SMOPS_fnames)
  print(SMOPS_fnames)
  print("---")

  do i=0,nSMOPS-1                                   ; loop over all files
     print("=======================================================")
     SMOPS_fi		= addfile(SMOPS_dir+SMOPS_fnames(i), "r")
     sm	               := SMOPS_fi->Blended_SM      ; type short 
     sm@_FillValue     := sm@FillValue              ; add correct attribute 
     sm                := short2flt( sm )
     sm!0		= "latitude"
     sm!1		= "longitude"
     sm&latitude	=  latitude
     sm&longitude	=  longitude 
     printVarSummary(sm)
     printMinMax(sm, 0)
     print("---")

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Use "poisson_grid_fill" to fill in the mosaic gaps
; .    https://www.ncl.ucar.edu/Applications/grid_fill.shtml
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

     sm_pois      := sm
     poisson_grid_fill(sm_pois, gtype, guess, nscan, eps, relc, iopt)
    ;printVarSummary(sm_pois)
    ;printMinMax(sm_pois, 0)
    ;print("---")

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Use NCl's coarse builtin function to mask out water
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
     sm_new  = sm_pois
     sm_new  = where(lsm.eq.0 , sm_new@_FillValue , sm_new)
     sm_new  = where(lsm.eq.2 , sm_new@_FillValue , sm_new)
     sm_new  = where(lsm.eq.4 , sm_new@_FillValue , sm_new)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; linint2_points_Wrap requires South-to-North ordering  
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
     sm_new  = sm_new(::-1,:)   ; linint2_points requires South-to-North order
     printVarSummary(sm_new)
     printMinMax(sm_new, 0)
     print("---")

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Limit the resion for efficiency
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

     extra = 3              ; extra space
     latS  = latS-extra
     latN  = latN+extra
     lonL  = lonL-extra
     lonR  = lonR+extra
     sm_region = sm_new({latS:latN},{lonL:lonR})
     printVarSummary(sm_region)
     printMinMax(sm_new, 0)
     print("---")

     SMOPS_fo = linint2_points(sm_region&longitude, sm_region&latitude, sm_region, True , stlon, stlat, 0)
     print("linint2_points: "+stlon+"   "+stlat+"  "+SMOPS_fo)
  end do