docs: def interp2d2d(x,y,Z,xout,yout,split_factor=1,**kwargs):
      """    
      INTERP2D2D : Interpolate a 2D matrix into another 2D matrix
      
      @param x: 1st dimension vector of size NX
      @param y: 2nd dimension vector of size NY
      @param Z: Array to interpolate (NXxNY)
      @param xout: 1st dimension vector of size NXout
      @param yout: 2nd dimension vector of size NYout
      
      @keyword split_factor:Nummber of times to split arrays. Nb of threads is equal to split_factor**2. 
      
      @return: Interpolated array (NXoutxNYout)
      
      @author: Renaud DUSSURGET, LER/PAC, Ifremer La Seyne
      """
      
      #Queued thread
      ##############
      class ThreadClass(threading.Thread):
          
          #We override the __init__ method
          def __init__(self, input_q,indices_q,result_q):
              threading.Thread.__init__(self)
              self.input = input_q
              self.indices = indices_q
              self.result = result_q
          
          def task(self,NaN=True):
              
              #grabs host from queue
              inargs = self.input.get()
              ind = inargs[0] #thread index
              x = inargs[1]
              y = inargs[2]
              Z = inargs[3]
              xout = inargs[4]
              yout = inargs[5]
              
              Zout=_interp2d2d(x,y,Z,xout,yout,**kwargs)
              
          #    Zout = sc.interpolate.griddata(points, gz, xi, **kwargs)
              
              self.indices.put(ind)
              self.result.put(Zout)
      
              del Zout
              
  #            if verbose > 0 : print "%s ended at time: %s" % (self.getName(), datetime.datetime.now())
          
          def run(self):
              
                  #Starts the queue
                  self.task()
                  #signals to queue job is done
                  self.input.task_done()
         
  
      #Setup input and output queues
      input_q = Queue.Queue()
      indices_q = Queue.Queue()
      result_q = Queue.Queue()
      
      #Map the data along X axis
      N_threads=split_factor**2
  #    over=np.ceil(cut/dx) #Overlay between each time series processed in parallel
      
      #Get dimensions to split matrix
      nxin=x.size
      nyin=y.size
      nxout=xout.size
      nyout=yout.size
      
      gx = np.reshape(np.repeat(x,nyin),(nxin,nyin))
      gy = np.reshape(np.repeat(y,nxin),(nyin,nxin)).transpose((1,0))
  
      gxout = np.reshape(np.repeat(xout,nyout),(nxout,nyout))
      gyout = np.reshape(np.repeat(yout,nxout),(nyout,nxout)).transpose((1,0))
      
      #Map output coordinates
      ind=[]
      
      #Map output coordinates
      xsplit=[]
      ysplit=[]
      for i in np.arange(split_factor) : 
          for j in np.arange(split_factor) :
              xsplit.append(i*(nxout/float(split_factor)))
              xsplit.append((i+1)*(nxout/float(split_factor))-1)
              ysplit.append(j*(nyout/float(split_factor)))
              ysplit.append((j+1)*(nyout/float(split_factor))-1)
      
      #Round
      xsplit=np.round(xsplit).astype(int)
      ysplit=np.round(ysplit).astype(int)
      
      N_threads = len(xsplit)/2
      
      
      
  #    th_xout=gxout[0:nxout/split_factor,0:2]
  #    th_yout=
  
  #    gz = Z.flatten()
       
  #    points = zip(*(gx.flatten(),gy.flatten())) 
      
      #spawn a pool of threads, and pass them queue instance 
      for i in np.arange(N_threads):
          t = ThreadClass(input_q,indices_q,result_q)
          t.setDaemon(True)
          t.start()
      
      #Feed threads with data
      for i in range(N_threads):
          xoind=xsplit[[i*2,(i*2)+1]]
          yoind=ysplit[[i*2,(i*2)+1]]
          xiind=x[xout[xoind].astype(int)].astype(int)
          yiind=y[yout[yoind].astype(int)].astype(int)
          th_x=x[xiind[0]:xiind[1]+1]
          th_y=y[yiind[0]:yiind[1]+1]
          th_xo=xout[xoind[0]:xoind[1]+1]
          th_yo=yout[yoind[0]:yoind[1]+1]
          th_z=Z[xiind[0]:xiind[1]+1,yiind[0]:yiind[1]+1]
  #        input_q.put((i,np.copy(points),gz.copy(),xout[iind[i]].copy(),yout[iind[i]].copy()))
          input_q.put((i,th_x,th_y,th_z,th_xo,th_yo))
  
     
  #    outvar=Z.copy()
  #    outvar.data[:]=outvar.fill_value
      
  #    for i in range(N_threads):
  #        print "%s launched time: %s" % (i, datetime.datetime.now())
  #        input_q.put((i,np.copy(points),gz.copy(),xout[iind[i]].copy(),yout[iind[i]].copy()))
  
      
  
      #wait on the queue until everything has been processed     
      input_q.join()
      
      #Sort threads
      tnb=[]
      for i in range(N_threads) : tnb.append(indices_q.get(i))
      tsort=np.argsort(tnb)
      
      #Get back the results for each thread in a list of results
      for i in np.arange(N_threads) :
          r=result_q.get(i)
          if i == 0 : dum = [r]
          else : dum.append(r)
  
      #Reorder data from each thread into output matrix
      for i in tsort :
          if i == tsort[0] : outmat = dum[i]
          else : outmat=np.ma.concatenate((outmat,dum[i]),0) if isinstance(outmat,np.ma.masked_array) else np.concatenate((outmat,dum[i]),0)
      
      if len(outmat) != len(outmat) :
          raise '[ERROR]Output array is not coherent with input array - check array reconstruction'
      return(outmat)