; Loops thru a 1D variable [eg, time] and finds
; the indices which match up with the values of the cvWant 1D array
; e.g. time_want = (/1948, 1957, 1964, 1965, 1989/)
; indWant   = get1DindexMultiple (time, time_want)
; note that the values of cvWant must EXIST in cv

undef("get1DindexMultiple")
function get1DindexMultiple (cv[*],  cvWant[*])
local nWant, indWant, n, kn, kind, nkk
begin
  ncv     = dimsizes (cv)
  nWant   = dimsizes (cvWant)
  indWant = new (nWant*ncv, integer)


  kn = -1
  do n=0,nWant-1
     kind = ind( cv.eq.cvWant(n) )
     if (.not.ismissing(kind(0))) then
         nkk = dimsizes(kind)
         do km=0,nkk-1
            kn = kn+1
            indWant(kn) = kind(km)
         end do
         delete( nkk )
      end if
      delete(kind)
  end do

  return (indWant(0:kn))
end

undef("sort_xy")
procedure sort_xy(x:numeric, y:numeric, opt:numeric)
local x_old, y_old, x_ind, y_ind
begin
   ;--- sort x ---
   if (opt .eq. 1) then
     x_old = x
     qsort(x)
     x_ind = get1DindexMultiple(x_old, x)
print(x_ind)
     y_old = y
print("===> dimsizes(y) <===")
print(dimsizes(y))
print("===> dimsizes(y(x_ind)) <===")
print(dimsizes(y(x_ind)))
     y     = y(x_ind)
   ;--- sort y ---
   else
     y_old = y
     qsort(y)
     y_ind = get1DindexMultiple(y_old, y)
print(y_ind)
     x_old = x
print(dimsizes(x))
print(dimsizes(x(y_ind)))
     x     = x(y_ind)
   end if
end

undef("SORT_XY")
function SORT_XY(x:numeric, y:numeric, opt:numeric)
local x_old, y_old, x_ind, y_ind
begin
   ;--- sort x ---
   if (opt .eq. 1) then
     x_old = x
     qsort(x)
     x_ind = get1DindexMultiple(x_old, x)
     y_old = y
     return( y(x_ind) )
   ;--- sort y ---
   else
     y_old = y
     qsort(y)
     y_ind = get1DindexMultiple(y_old, y)
     x_old = x
     return( x(y_ind) )
   end if
end


       xxx = (/ 3, 5, 2, 10, 6, 4, 3 /)
       yyy = (/ 1, 8, 4, 3 , 7, 6, 2 /)

       qqq = SORT_XY(xxx, yyy, 1)
       print(qqq)
       print("=======================") 
       sort_xy(xxx, yyy, 1)