Index: trunk/tools/module_ForInterpolate.F90
===================================================================
--- trunk/tools/module_ForInterpolate.F90	(revision 1172)
+++ trunk/tools/module_ForInterpolate.F90	(revision 1173)
@@ -1,5 +1,11 @@
 ! Module to interpolate values from a giving projection and pressure interpolation
 ! To be included in a python
-! f2py -m module_ForInterpolate --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForInterpolate.F90  >& run_f2py.log
+! Follow compilation instructions from Makefile
+! Content
+!  LlInterpolateProjection: Subroutine which provides the indices for a given interpolation of a projection
+!  var2D_IntProj: Subroutine to interpolate a 2D variable
+!  var3D_IntProj: Subroutine to interpolate a 3D variable
+!  var4D_IntProj: Subroutine to interpolate a 4D variable
+!  var5D_IntProj: Subroutine to interpolate a 5D variable
 MODULE module_ForInterpolate
 
@@ -714,4 +720,399 @@
 END SUBROUTINE CoarseInterpolateExact
 
+SUBROUTINE LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy,   &
+  pdimx, pdimy)
+! Subroutine which provides the indices for a given interpolation of a projection
+
+  USE module_generic
+
+  IMPLICIT NONE
+
+  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
+  INTEGER, INTENT(in)                                    :: idimx, idimy, pdimx, pdimy
+  REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in)          :: projlon, projlat
+  REAL(r_k), DIMENSION(idimx,idimy), INTENT(in)          :: inlonv, inlatv
+  CHARACTER(LEN=50), INTENT(in)                          :: intkind
+  REAL(r_k), DIMENSION(3,16,pdimx,pdimy), INTENT(out)    :: outLlw
+
+! Local
+  INTEGER                                                :: i,j,iv, ix, iy
+  REAL(r_k)                                              :: mindiffLl, dist
+  REAL(r_k), DIMENSION(idimx,idimy)                      :: difflonlat
+  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
+  INTEGER, DIMENSION(2)                                  :: iLl
+  CHARACTER(LEN=50)                                      :: fname
+
+!!!!!!! Variables
+! idimx, idimy: dimension length of the input projection
+! pdimx, pdimy: dimension length of the target projection
+! in[lon/lat]: longitudes and latitudes of the target projection
+! proj[lon/lat]: longitudes and latitudes of the target projection
+! intkind: kind of interpolation
+!   'npp': nearest neighbourgh
+!   'dis': weighted distance, grid-output for SW, NW, NE, SE
+! outLlw: output interpolation result
+!    for point pi,pj; up to 16 different values of
+!       1st: i-index in input projection
+!       2nd: j-index in input projection
+!       3rd: weight for i-index, j-index to use for ponderation (0<1.)
+  fname = 'LlInterpolateProjection'
+
+  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
+  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
+
+  ! Using case outside loop to be more efficient
+  SELECT CASE(TRIM(intkind))
+    CASE ('dis')
+      DO i=1, pdimx
+        DO j=1, pdimy      
+          ! Find point
+          difflonlat = SQRT((projlon(i,j)-inlonv)**2. + (projlat(i,j)-inlatv)**2.)
+          mindiffLl = MINVAL(difflonlat)
+          ! Getting the four surrounding values
+          iLl = index2DArrayR(difflonlat, idimx, idimy, mindiffLl)
+          IF ( projlon(i,j) < inlonv(iLl(1),iLl(2)) ) THEN
+            outLlw(1,1,i,j) = MAX(iLl(1)-1,1)
+            outLlw(1,2,i,j) = MAX(iLl(1)-1,1)
+            outLlw(1,3,i,j) = MIN(iLl(1),idimx)
+            outLlw(1,4,i,j) = MIN(iLl(1),idimx)
+          ELSE
+            outLlw(1,1,i,j) = MIN(iLl(1),1)
+            outLlw(1,2,i,j) = MIN(iLl(1),1)
+            outLlw(1,3,i,j) = MAX(iLl(1)+1,idimx)
+            outLlw(1,4,i,j) = MAX(iLl(1)+1,idimx)
+          END IF
+          IF ( projlat(i,j) < inlatv(iLl(1),iLl(2)) ) THEN
+            outLlw(2,1,i,j) = MIN(iLl(1)-1,1)
+            outLlw(2,2,i,j) = MAX(iLl(1),idimy)
+            outLlw(2,3,i,j) = MAX(iLl(1),idimy)
+            outLlw(2,4,i,j) = MIN(iLl(1)-1,1)
+          ELSE
+            outLlw(2,1,i,j) = MIN(iLl(1),1)
+            outLlw(2,2,i,j) = MAX(iLl(1)+1,idimy)
+            outLlw(2,3,i,j) = MAX(iLl(1)+1,idimy)
+            outLlw(2,4,i,j) = MIN(iLl(1),1)
+          END IF
+          ! Computing distances
+          DO iv=1,4
+            ix = INT(outLlw(1,iv,i,j))
+            iy = INT(outLlw(2,iv,i,j))
+            dist = SQRT( (projlon(i,j)-inlonv(ix,iy))**2 + (projlat(i,j)-inlatv(ix,iy))**2 )
+            IF ( dist /= 0.) THEN
+              outLlw(3,iv,i,j) = 1./dist
+            ELSE
+              outLlw(3,iv,i,j) = 1.
+            END IF
+          END DO
+          ! Normalizing
+          outLlw(3,:,i,j) = outLlw(3,:,i,j)/(SUM(outLlw(3,:,i,j)))
+        END DO
+      END DO
+    CASE ('npp')
+      DO i=1, pdimx
+        DO j=1, pdimy      
+          ! Find point
+          difflonlat = SQRT((projlon(i,j)-inlonv)**2. + (projlat(i,j)-inlatv)**2.)
+          mindiffLl = MINVAL(difflonlat)
+          outLlw(1:2,1,i,j) = index2DArrayR(difflonlat, idimx, idimy, mindiffLl)
+          outLlw(3,1,i,j) = 1.
+        END DO
+      END DO
+  END SELECT
+
+END SUBROUTINE LlInterpolateProjection
+
+SUBROUTINE var2D_IntProj(var2Din, inlonv, inlatv, projlon, projlat, intkind, varout, idimx, idimy,    &
+  pdimx, pdimy)
+! Subroutine to interpolate a 2D variable
+
+  USE module_generic
+
+  IMPLICIT NONE
+
+  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
+  INTEGER, INTENT(in)                                    :: idimx, idimy, pdimx, pdimy
+  REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in)          :: projlon, projlat
+  REAL(r_k), DIMENSION(idimx,idimy), INTENT(in)          :: inlonv, inlatv
+  CHARACTER(LEN=50), INTENT(in)                          :: intkind
+  REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in)          :: var2Din
+  REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(out)         :: varout
+
+! Local
+  INTEGER                                                :: i, j, iv, ix, iy
+  REAL(r_k)                                              :: w
+  REAL(r_k), DIMENSION(3,16,pdimx,pdimy)                 :: outLlw
+  CHARACTER(LEN=50)                                      :: fname
+
+!!!!!!! Variables
+! idimx, idimy: dimension length of the input projection
+! pdimx, pdimy: dimension length of the target projection
+! in[lon/lat]: longitudes and latitudes of the target projection
+! proj[lon/lat]: longitudes and latitudes of the target projection
+! intkind: kind of interpolation
+!   'npp': nearest neighbourgh
+!   'dis': weighted distance, grid-output for SW, NW, NE, SE
+! outLlw: output interpolation result
+!    for point pi,pj; up to 16 different values of
+!       1st: i-index in input projection
+!       2nd: j-index in input projection
+!       3rd: weight for i-index, j-index to use for ponderation (0<1.)
+! var2Din: 2D variable to interpolate
+! varout: variable interpolated on the target projection
+  fname = 'var2D_IntProj'
+
+  CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,&
+    pdimy)
+
+  SELECT CASE (intkind) 
+    CASE('dis')
+      DO i=1, pdimx
+        DO j=1, pdimx
+          varout(i,j) = 0.
+          DO iv=1, 4
+            ix = INT(outLlw(1,iv,i,j))
+            iy = INT(outLlw(2,iv,i,j))
+            w =  outLlw(3,iv,i,j)
+            varout(i,j) = varout(i,j) + w*var2Din(ix,iy)
+          END DO
+        END DO
+      END DO
+    CASE('npp')
+      DO i=1, pdimx
+        DO j=1, pdimx
+          ix = INT(outLlw(1,1,i,j))
+          iy = INT(outLlw(2,1,i,j))
+          varout(i,j) = var2Din(ix,iy)
+        END DO
+      END DO
+  END SELECT
+
+END SUBROUTINE var2D_IntProj
+
+SUBROUTINE var3D_IntProj(var3Din, inlonv, inlatv, projlon, projlat, intkind, varout, idimx, idimy,    &
+  pdimx, pdimy, d3)
+! Subroutine to interpolate a 3D variable
+
+  USE module_generic
+
+  IMPLICIT NONE
+
+  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
+  INTEGER, INTENT(in)                                    :: idimx, idimy, pdimx, pdimy, d3
+  REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in)          :: projlon, projlat
+  REAL(r_k), DIMENSION(idimx,idimy), INTENT(in)          :: inlonv, inlatv
+  CHARACTER(LEN=50), INTENT(in)                          :: intkind
+  REAL(r_k), DIMENSION(pdimx,pdimy,d3), INTENT(in)       :: var3Din
+  REAL(r_k), DIMENSION(pdimx,pdimy,d3), INTENT(out)      :: varout
+
+! Local
+  INTEGER                                                :: i, j, k, iv, ix, iy
+  REAL(r_k)                                              :: w
+  REAL(r_k), DIMENSION(3,16,pdimx,pdimy)                 :: outLlw
+  CHARACTER(LEN=50)                                      :: fname
+
+!!!!!!! Variables
+! idimx, idimy: dimension length of the input projection
+! pdimx, pdimy: dimension length of the target projection
+! in[lon/lat]: longitudes and latitudes of the target projection
+! proj[lon/lat]: longitudes and latitudes of the target projection
+! intkind: kind of interpolation
+!   'npp': nearest neighbourgh
+!   'dis': weighted distance, grid-output for SW, NW, NE, SE
+! outLlw: output interpolation result
+!    for point pi,pj; up to 16 different values of
+!       1st: i-index in input projection
+!       2nd: j-index in input projection
+!       3rd: weight for i-index, j-index to use for ponderation (0<1.)
+! var3Din: 3D variable to interpolate
+! varout: variable interpolated on the target projection
+  fname = 'var3D_IntProj'
+
+  CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,&
+    pdimy)
+
+  SELECT CASE (intkind) 
+    CASE('dis')
+      DO i=1, pdimx
+        DO j=1, pdimx
+          DO k=1, d3
+            varout(i,j,k) = 0.
+            DO iv=1, 4
+              ix = INT(outLlw(1,iv,i,j))
+              iy = INT(outLlw(2,iv,i,j))
+              w =  outLlw(3,iv,i,j)
+              varout(i,j,k) = varout(i,j,k) + w*var3Din(ix,iy,k)
+            END DO
+          END DO
+        END DO
+      END DO
+    CASE('npp')
+      DO i=1, pdimx
+        DO j=1, pdimx
+          DO k=1, d3
+            ix = INT(outLlw(1,1,i,j))
+            iy = INT(outLlw(2,1,i,j))
+            varout(i,j,k) = var3Din(ix,iy,k)
+          END DO
+        END DO
+      END DO
+  END SELECT
+
+END SUBROUTINE var3D_IntProj
+
+SUBROUTINE var4D_IntProj(var4Din, inlonv, inlatv, projlon, projlat, intkind, varout, idimx, idimy,    &
+  pdimx, pdimy, d3, d4)
+! Subroutine to interpolate a 4D variable
+
+  USE module_generic
+
+  IMPLICIT NONE
+
+  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
+  INTEGER, INTENT(in)                                    :: idimx, idimy, pdimx, pdimy, d3, d4
+  REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in)          :: projlon, projlat
+  REAL(r_k), DIMENSION(idimx,idimy), INTENT(in)          :: inlonv, inlatv
+  CHARACTER(LEN=50), INTENT(in)                          :: intkind
+  REAL(r_k), DIMENSION(pdimx,pdimy,d3,d4), INTENT(in)    :: var4Din
+  REAL(r_k), DIMENSION(pdimx,pdimy,d3,d4), INTENT(out)   :: varout
+
+! Local
+  INTEGER                                                :: i, j, k, l, iv, ix, iy
+  REAL(r_k)                                              :: w
+  REAL(r_k), DIMENSION(3,16,pdimx,pdimy)                 :: outLlw
+  CHARACTER(LEN=50)                                      :: fname
+
+!!!!!!! Variables
+! idimx, idimy: dimension length of the input projection
+! pdimx, pdimy: dimension length of the target projection
+! in[lon/lat]: longitudes and latitudes of the target projection
+! proj[lon/lat]: longitudes and latitudes of the target projection
+! intkind: kind of interpolation
+!   'npp': nearest neighbourgh
+!   'dis': weighted distance, grid-output for SW, NW, NE, SE
+! outLlw: output interpolation result
+!    for point pi,pj; up to 16 different values of
+!       1st: i-index in input projection
+!       2nd: j-index in input projection
+!       3rd: weight for i-index, j-index to use for ponderation (0<1.)
+! var4Din: 4D variable to interpolate
+! varout: variable interpolated on the target projection
+  fname = 'var4D_IntProj'
+
+  CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,&
+    pdimy)
+
+  SELECT CASE (intkind) 
+    CASE('dis')
+      DO i=1, pdimx
+        DO j=1, pdimx
+          DO k=1, d3
+            DO l=1, d4
+              varout(i,j,k,l) = 0.
+              DO iv=1, 4
+                ix = INT(outLlw(1,iv,i,j))
+                iy = INT(outLlw(2,iv,i,j))
+                w =  outLlw(3,iv,i,j)
+                varout(i,j,k,l) = varout(i,j,k,l) + w*var4Din(ix,iy,k,l)
+              END DO
+            END DO
+          END DO
+        END DO
+      END DO
+    CASE('npp')
+      DO i=1, pdimx
+        DO j=1, pdimx
+          ix = INT(outLlw(1,1,i,j))
+          iy = INT(outLlw(2,1,i,j))
+          DO k=1, d3
+            DO l=1, d4
+              varout(i,j,k,l) = var4Din(ix,iy,k,l)
+            END DO
+          END DO
+        END DO
+      END DO
+  END SELECT
+
+END SUBROUTINE var4D_IntProj
+
+
+SUBROUTINE var5D_IntProj(var5Din, inlonv, inlatv, projlon, projlat, intkind, varout, idimx, idimy,    &
+  pdimx, pdimy, d3, d4, d5)
+! Subroutine to interpolate a 5D variable
+
+  USE module_generic
+
+  IMPLICIT NONE
+
+  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
+  INTEGER, INTENT(in)                                    :: idimx, idimy, pdimx, pdimy, d3, d4, d5
+  REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in)          :: projlon, projlat
+  REAL(r_k), DIMENSION(idimx,idimy), INTENT(in)          :: inlonv, inlatv
+  CHARACTER(LEN=50), INTENT(in)                          :: intkind
+  REAL(r_k), DIMENSION(pdimx,pdimy,d3,d4,d5), INTENT(in) :: var5Din
+  REAL(r_k), DIMENSION(pdimx,pdimy,d3,d4,d5), INTENT(out):: varout
+
+! Local
+  INTEGER                                                :: i, j, k, l, m, iv, ix, iy
+  REAL(r_k)                                              :: w
+  REAL(r_k), DIMENSION(3,16,pdimx,pdimy)                 :: outLlw
+  CHARACTER(LEN=50)                                      :: fname
+
+!!!!!!! Variables
+! idimx, idimy: dimension length of the input projection
+! pdimx, pdimy: dimension length of the target projection
+! in[lon/lat]: longitudes and latitudes of the target projection
+! proj[lon/lat]: longitudes and latitudes of the target projection
+! intkind: kind of interpolation
+!   'npp': nearest neighbourgh
+!   'dis': weighted distance, grid-output for SW, NW, NE, SE
+! outLlw: output interpolation result
+!    for point pi,pj; up to 16 different values of
+!       1st: i-index in input projection
+!       2nd: j-index in input projection
+!       3rd: weight for i-index, j-index to use for ponderation (0<1.)
+! var5Din: 5D variable to interpolate
+! varout: variable interpolated on the target projection
+  fname = 'var5D_IntProj'
+
+  CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,&
+    pdimy)
+
+  SELECT CASE (intkind) 
+    CASE('dis')
+      DO i=1, pdimx
+        DO j=1, pdimx
+          DO k=1, d3
+            DO l=1, d4
+              DO m=1, d5
+              varout(i,j,k,l,m) = 0.
+                DO iv=1, 4
+                  ix = INT(outLlw(1,iv,i,j))
+                  iy = INT(outLlw(2,iv,i,j))
+                  w =  outLlw(3,iv,i,j)
+                  varout(i,j,k,l,m) = varout(i,j,k,l,m) + w*var5Din(ix,iy,k,l,m)
+                END DO
+              END DO
+            END DO
+          END DO
+        END DO
+      END DO
+    CASE('npp')
+      DO i=1, pdimx
+        DO j=1, pdimx
+          ix = INT(outLlw(1,1,i,j))
+          iy = INT(outLlw(2,1,i,j))
+          DO k=1, d3
+            DO l=1, d4
+              DO m=1, d5
+                varout(i,j,k,l,m) = var5Din(ix,iy,k,l,m)
+              END DO
+            END DO
+          END DO
+        END DO
+      END DO
+  END SELECT
+
+END SUBROUTINE var5D_IntProj
+
 SUBROUTINE Interpolate(projlon, projlat, lonvs, latvs, mindiff, inpt, diffs, ilonlat, dimx, dimy,     &
   Ninpts)
Index: trunk/tools/nc_var.py
===================================================================
--- trunk/tools/nc_var.py	(revision 1172)
+++ trunk/tools/nc_var.py	(revision 1173)
@@ -44,8 +44,8 @@
   'get_namelist_vars', 'get_Variables',                                              \
   'getvalues_lonlat', 'grattr',                                                      \
-  'grmattr', 'idims', 'igattrs', 'increaseDimvar', 'isgattrs', 'isvattrs', 'ivars',  \
-  'ivattrs', 'LMDZ_toCF', 'maskvar', 'model_characteristics',                        \
+  'grmattr', 'idims', 'igattrs', 'increaseDimvar', 'isgattrs',                       \
+  'isvattrs', 'ivars', 'ivattrs', 'LMDZ_toCF', 'maskvar', 'model_characteristics',   \
   'ncreplace', 'ncstepdiff', 'netcdf_concatenation', 'netcdf_fold_concatenation',    \
-  'netcdf_fold_concatenation_HMT', 'Partialmap_Entiremap',                           \
+  'netcdf_fold_concatenation_HMT', 'reproject', 'Partialmap_Entiremap',              \
   'Partialmap_EntiremapFor', 'Partialmap_EntiremapForExact',                         \
   'pinterp', 'remapnn',                                                              \
@@ -247,4 +247,6 @@
 elif oper == 'Partialmap_EntiremapForExact':
     ncvar.Partialmap_EntiremapForExact(opts.values, opts.ncfile, opts.varname)
+elif oper == 'reproject':
+    ncvar.reproject(opts.values, opts.ncfile, opts.varname)
 elif oper == 'seasmean':
     ncvar.seasmean(timename, opts.ncfile, opts.varname)
Index: trunk/tools/nc_var_tools.py
===================================================================
--- trunk/tools/nc_var_tools.py	(revision 1172)
+++ trunk/tools/nc_var_tools.py	(revision 1173)
@@ -90,4 +90,5 @@
 # remapnn: Function to remap to the nearest neightbor a variable using projection from another file
 # remapnn_old: Function to remap to the nearest neightbor a variable using projection from another fil
+# reproject: Function to reproject values to another one
 # seasmean: Function to compute the seasonal mean of a variable
 # sellonlatbox: Function to select a lotlan box from a data-set
@@ -15782,4 +15783,250 @@
 #Partialmap_Entiremap('longitude,latitude,std,5000.,lonlat', '/home/lluis/etudes/DYNAMICO/ORCHIDEE/interpolation/data/carteveg5km.nc', 'vegetation_map')
 
+def reproject(values, filen, variable):
+    """ Function to reproject values to another one
+      reproject(values, filen, variable)
+      values=[inlonn],[inlatn],[projfile],[projlonn],[projlatn],[kind],[dimvars]
+        inlonn= name of the longitude values in the file to interpolate
+        inlatn= name of the latitude values in the file to interpolate
+        projlonn= name of the longitude values in the file with the target projection
+        projlatn= name of the latitude values in the file with the target projection
+        kind= kind of interpolation
+          'npp': nearest point
+          'dis': distance weighted (4 closest)
+        dimvars= ':', separated list of couples of dimension, and variable with its values [dimn]@[variabledim]
+      filen= name of file to interpolate
+      variable = ',' list of variable to interpolate ('all' for all variables)
+    """
+    import module_ForInt as fin
+    fname = 'reproject'
+
+    ofile = fname + '.nc'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print reproject.__doc__
+        quit()
+
+    arguments = '[inlonn],[inlatn],[projfile],[projlonn],[projlatn],[kind],[dimvars]'
+    gen.check_arguments(fname, values, arguments, ',')
+    
+    inlonn = values.split(',')[0]
+    inlatn = values.split(',')[1]
+    projfile = values.split(',')[2]
+    projlonn = values.split(',')[3]
+    projlatn = values.split(',')[4]
+    kind = values.split(',')[5]
+    dimvars = gen.Str_list(values.split(',')[6],':')
+
+    varns = gen.Str_list(variable,',')
+
+    inf = NetCDFFile(filen, 'r')
+    if not os.path.isfile(projfile):
+        print errormsg
+        print '  ' + fname + ": file with target projection '" + projfile +          \
+          "' does not exist !!"
+        quit(-1)
+    projf = NetCDFFile(projfile, 'r')
+
+    if not inf.variables.has_key(inlonn):
+        print errormsg
+        print '  ' +fname+ ": input file '" + filen + "' does not have longitude '" +\
+          inlonn + "' !!"
+        quit(-1)
+    if not inf.variables.has_key(inlatn):
+        print errormsg
+        print '  ' +fname+ ": input file '" + filen + "' does not have latitude '" + \
+          inlatn + "' !!"
+        quit(-1)
+    if not projf.variables.has_key(projlonn):
+        print errormsg
+        print '  ' + fname + ": file with target projection '" + projfile +          \
+          "' does not have longitude '" + projlonn + "' !!"
+        quit(-1)
+    if not projf.variables.has_key(projlatn):
+        print errormsg
+        print '  ' + fname + ": file with target projection '" + projfile +          \
+          "' does not have latitude '" + projlatn + "' !!"
+        quit(-1)
+
+    oinlon = inf.variables[inlonn]
+    oinlat = inf.variables[inlatn]
+    oprojlon = projf.variables[projlonn]
+    oprojlat = projf.variables[projlatn]
+
+    inlonvals0 = oinlon[:]
+    inlatvals0 = oinlat[:]
+    olonvals0 = oprojlon[:]
+    olatvals0 = oprojlat[:]
+
+    # Getting 2D longitudes and latitudes
+    ilonvals, ilatvals = gen.lonlat2D(ilonvals0, ilatvals0)
+    olonvals, olatvals = gen.lonlat2D(olonvals0, olatvals0)
+
+    if len(ilonvals.shape) != 2:
+        print errormsg
+        print '  ' + fname + ': wrong input longitude rank:',  ilonvals.shape,       \
+          'it must be 2!!'
+        quit(-1)
+
+    if len(ilatvals.shape) != 2:
+        print errormsg
+        print '  ' + fname + ': wrong input latitude rank:',  ilatvals.shape,        \
+          'it must be 2!!'
+        quit(-1)
+
+    if len(olonvals.shape) != 2:
+        print errormsg
+        print '  ' + fname + ': wrong output longitude rank:', olonvals.shape,       \
+          'it must be 2!!'
+        quit(-1)
+
+    if len(olatvals.shape) != 2:
+        print errormsg
+        print '  ' + fname + ': wrong output latitude rank:', olatvals.shape,        \
+          'it must be 2!!'
+        quit(-1)
+
+    for i in range(2):
+        if ilonvals[i] != ilatvals[i]:
+            print errormsg
+            print '  ' + fname + ': dimension',i,'of input longitude:', ilonvals[i], \
+              'differs from latitude:', ilatvals[i], 'they must be the same!!'
+            quit(-1)
+        if rlonvals[i] != rlatvals[i]:
+            print errormsg
+            print '  ' + fname + ': dimension',i,'of output longitude:', ilonvals[i],\
+              'differs from latitude:', ilatvals[i], 'they must be the same!!'
+            quit(-1)
+    
+    idx = ilonvals.shape[1]
+    idy = ilonvals.shape[0]
+    rdx = rlonvals.shape[1]
+    rdy = rlonvals.shape[0]
+
+    ilon = ilonvals.astype('float64')
+    ilat = ilatvals.astype('float64')
+    rlon = rlonvals.astype('float64')
+    rlat = rlatvals.astype('float64')
+
+    ilont = ilon.transpose()
+    ilatt = ilat.transpose()
+    rlont = rlon.transpose()
+    rlatt = rlat.transpose()
+
+    # dimension-vardims
+    vardimns = {}
+    for vardn in dimvars:
+        dimn = vardn.split('@')[0]
+        dimvn = vardn.split('@')[1]
+        vardimns[dimn] = dimvn
+
+    # File creation
+    onewnc = NetCDFFile(ofile, 'w')
+
+    # Dimensions
+    newdim = onewnc.createDimension('lon',rdx)
+    newdim = onewnc.createDimension('lat',rdy)
+
+    # Which variables to re-project
+    if varns[0] == 'all':
+        varns = inf.variables.keys()
+        # No interpolation of lon, lat
+        varns.remove(inlonn)
+        varns.remove(inlatn)
+
+    for varn in varns:
+        if not inf.variables.has_key(inlatn):
+            print errormsg
+            print '  ' + fname + ": input file '" + filen + "' does not have '" +    \
+              varn + "' !!"
+            quit(-1)
+        print '  ' + fname + ": re-projecting '" + varn + "' ..."
+        ovarin = inf.variables[varn]
+        varin = ovarin[:]
+
+        Nvarind = len(varin.shape)
+        if Nvarind < 2:
+            print errormsg
+            print '  ' + fname + ': wrong input data rank:', varin.shape,            \
+              'it must be at least 2!!'
+            quit(-1)
+
+        varin = varin.astype('float64')
+        varint = varin.transpose()
+
+        if Nvarind == 2:
+            varout = fin.var2D_IntProj(var2Din=varint, inlonv=inlont, inlatv=inlatt, \
+              projlon=rlont, projlat=rlatt, intkind=kind, idimx=idx, idimy=idy,      \
+              pdimx=rdx, pdimy=rdy)
+        elif Nvarind == 3:
+            varout = fin.var3D_IntProj(var3Din=varint, inlonv=inlont, inlatv=inlatt, \
+              projlon=rlont, projlat=rlatt, intkind=kind, idimx=idx, idimy=idy,      \
+              pdimx=rdx, pdimy=rdy, d3=varin.shape[0])
+        elif Nvarind == 4:
+            varout = fin.var4D_IntProj(var4Din=varint, inlonv=inlont, inlatv=inlatt, \
+              projlon=rlont, projlat=rlatt, intkind=kind, idimx=idx, idimy=idy,      \
+              pdimx=rdx, pdimy=rdy, d3=varin.shape[0], d4=varin.shape[1])
+        elif Nvarind == 5:
+            varout = fin.var4D_IntProj(var4Din=varint, inlonv=inlont, inlatv=inlatt, \
+              projlon=rlont, projlat=rlatt, intkind=kind, idimx=idx, idimy=idy,      \
+              pdimx=rdx, pdimy=rdy, d3=varin.shape[0], d4=varin.shape[1],            \
+              d5=varin.shape[2])
+        else:
+            print errormsg
+            print '  ' + fname + ': rank:', Nvarind, 'for input data not ready !!'
+            quit(-1)
+    
+        # Variable's dimensions and their related variables
+        dimsfile = onewnc.dimensions()
+        for dimn in ovarin.dimensions():
+            if not gen.searchInlist(dimsfile, dimn):
+                odimn = inf.dimensions[dimn]
+                if odimn.isunlimited():
+                    newdim = onewnc.createDimension(dimn, None)
+                else:
+                    newdim = onewnc.createDimension(dimn, len(odimn))
+                if not vardimsn.has_key(dimn):
+                    print errormsg
+                    print '  ' + fname + ": dimension '" + dimn + "' does not have "+\
+                      'assigned any variable !!'
+                    print '    assigned pairs (via [dimvars]) _______'
+                    gen.print_dictionary(vardimsn)
+                    quit(-1)
+                else:
+                    add_vars(inf,onewnc,vardimsn[dimn])
+
+        # Interpolated variable
+        vardims = ovarin.dimensions()
+        newvardims = gen.replace_list(vardims,inlonn,'lon')
+        newvardims = gen.replace_list(newvardims,inlatn,'lat')
+        varattrs = ovarin.ncattrs()
+        if gen.searchInlist(varattrs,'_FillValue'):        
+            fillv = ovarin.getncattr('_FillValue')
+            newvar = onewnc.createVariable(varn, 'f4', tuple(newvardims),            \
+              fill_value=fillv)
+            varattrs.remove('_FillValue')
+        else:
+            newvar = onewnc.createVariable(varn, 'f4', tuple(newvardims))
+
+        newvar[:] = varout[:]
+        for attrn in varattrs:
+            attrv = ovarin.getncattr(attrn)
+            set_attribute(newvar, attrn, attrv)
+
+        onewnc.sync()
+
+    # Global attributes
+    add_globalattrs(inf,onewnc,'all')
+    onewnc.sync()
+
+    inf.close()
+    projf.close()
+    onewnc.close()
+
+    print fname + ": successful re-projection saved in '" + ofile + "' !!"
+
+    return
+
 def lonlatvarsFile(lonvals,latvals,dnx,dny,oLlnc):
     """ Function to CF-define the longitudes and latitudes variables within a file
