load "./shapefile_mask_data.ncl"
load "./shapefile_utils.ncl"

;--- Open WRF output file.

  filename         = ("ref.nc")
  a                = addfile(filename,"r")
  var              = a->pr 
  var1             = dim_avg_n_Wrap(var, 0) 
  var@_FillValue   = -9999
                     

  opt              = True
  opt@return_mask  = True
  var_mask_2d      = shapefile_mask_data(var1,"cuencasx.shp",opt)  ;usando var1 ou var


;---Promote to a 3D array so we can mask all of "var"

  var_mask_3d       = conform_dims(dimsizes(var),var_mask_2d,(/1,2/))
                      copy_VarMeta(var(:,:,:),var_mask_3d)    
                      printMinMax( var_mask_3d, 0)

;---Apply the 3D mask to original variable

  var_mask          = where(var_mask_3d.eq.1,var,var@_FillValue)
  print("Min/max after masking: " + min(var_mask) + "/" + max(var_mask))
  copy_VarMeta(var(:,:,:),var_mask) 

WRITE_NETCDF = True
if(WRITE_NETCDF) then
  fout = "ref_mask.nc"    
  system("rm -f " + fout)
  ncdf = addfile(fout, "c")
  ncdf->pr =  var_mask
  print("")
  print(" data saved to " + fout)
end if