Index: trunk/tools/module_ForInterpolate.F90
===================================================================
--- trunk/tools/module_ForInterpolate.F90	(revision 739)
+++ trunk/tools/module_ForInterpolate.F90	(revision 740)
@@ -743,6 +743,6 @@
 END SUBROUTINE CoarseInterpolateExact
 
-SUBROUTINE Interpolate(projlon, projlat, lonvs, latvs, mindiff, ivar, newvar, newvarin, newvarinpt,   &
-  newvarindiff, dimx, dimy, Ninpts)
+SUBROUTINE Interpolate(projlon, projlat, lonvs, latvs, mindiff, ivar, newvarinpt,                     &
+  newvarindiff, ilonlat, dimx, dimy, Ninpts)
 ! Subroutine which finds the closest grid point within a projection
 
@@ -757,16 +757,14 @@
   REAL(r_k), DIMENSION(Ninpts), INTENT(in)               :: ivar, lonvs, latvs
   REAL(r_k), INTENT(in)                                  :: mindiff
-  REAL(r_k), DIMENSION(dimx,dimy), INTENT(out)           :: newvar
-  INTEGER, DIMENSION(dimx,dimy), INTENT(out)             :: newvarin
   INTEGER, DIMENSION(Ninpts), INTENT(out)                :: newvarinpt
   REAL(r_k), DIMENSION(Ninpts), INTENT(out)              :: newvarindiff
+  INTEGER, DIMENSION(Ninpts,2), INTENT(out)              :: ilonlat
 
 ! Local
   INTEGER                                                :: iv,i,j
   INTEGER                                                :: ierr
-  INTEGER, DIMENSION(2)                                  :: ilonlat
-  REAL(r_k), DIMENSION(dimx,dimy)                        :: difflonlat
   REAL(r_k)                                              :: mindiffLl
   INTEGER                                                :: Ninpts1
+  REAL(r_k), DIMENSION(dimx,dimy)                        :: difflonlat
   REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
   CHARACTER(LEN=50)                                      :: fname
@@ -783,4 +781,5 @@
 ! newvarinpt: integer value indicating if the value has been already located (0: no, 1: yes)
 ! newvarindiff: distance of point from the input data to the closest target point
+! ilonlat: longitude and latitude of the point
 ! ncid: netCDF output file id
 
@@ -813,12 +812,10 @@
       difflonlat = SQRT((projlon-lonvs(iv))**2. + (projlat-latvs(iv))**2.)
       mindiffLl = MINVAL(difflonlat)
-      ilonlat = index2DArrayR(difflonlat, dimx, dimy, mindiffLl)
+      ilonlat(iv,:) = index2DArrayR(difflonlat, dimx, dimy, mindiffLl)
 
       IF (mindiffLl <= mindiff) THEN 
 !        percendone(iv,Ninpts,0.5,'done:')
 
-        IF (ilonlat(1) >= 0 .AND. ilonlat(1) >= 0) THEN
-          newvar(ilonlat(1),ilonlat(2)) = ivar(iv)
-          newvarin(ilonlat(1),ilonlat(2)) = iv
+        IF (ilonlat(iv,1) >= 0 .AND. ilonlat(iv,2) >= 0) THEN
           newvarinpt(iv) = 1
           newvarindiff(iv) = mindiffLl
@@ -830,5 +827,5 @@
           PRINT *,'  ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv),        &
             ' not relocated !!'
-          PRINT *,'    mindiffl:', mindiffLl, ' ilon:', ilonlat(1), ' ilat:', ilonlat(2)
+          PRINT *,'    mindiffl:', mindiffLl, ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2)
           STOP
         END IF
Index: trunk/tools/nc_var_tools.py
===================================================================
--- trunk/tools/nc_var_tools.py	(revision 739)
+++ trunk/tools/nc_var_tools.py	(revision 740)
@@ -20405,16 +20405,12 @@
 #              lonvss, latvss, fracd, mindiff, ovars)
 # Slow way
-            newvar,newvarin,newvarinpt[pts],newvarindiff[pts] =                      \
+            newvarinpt[pts],newvarindiff[pts],ilonlatv =                             \
               fin.module_forinterpolate.interpolate(projlon, projlat,                \
               lonvss, latvss, np.float64(mindiff), ovars)
+            for i in range(fracs):
+#                print 'Lluis: ilonlatv:', ilonlatv[i,:]
+                newvar[ilonlatv[i,1],ilonlatv[i,0]] = ovars[i]
+                newvarin[ilonlatv[i,1],ilonlatv[i,0]] = i
             newnc.sync()
-            print newvar 
-            print 'newvar:',type(newvar), newvar.shape, np.min(newvar), np.max(newvar)
-            print 'newvarin:',type(newvarin), newvarin.shape, np.min(newvarin), np.max(newvarin)
-            newnc.close()
-            for i in range(iri,ire,1):
-                print lonvss[i], latvss[i], ovars[i] 
-    
-            quit()
 
     elif kind == 'lonlat':
