Index: /trunk/tools/module_ForInterpolate.F90
===================================================================
--- /trunk/tools/module_ForInterpolate.F90	(revision 715)
+++ /trunk/tools/module_ForInterpolate.F90	(revision 715)
@@ -0,0 +1,297 @@
+! Module to interpolate values from a giving projection
+! To be included in a python
+! f2py -m module_ForInterpolate --f90exec=/usr/bin/gfortran-4.7 -c module_ForInterpolate.F90 module_generic.F90 >& run_f2py.log
+MODULE module_ForInterpolate
+
+  CONTAINS
+
+SUBROUTINE CoarselonlatFind(dx, dy, ilon, ilat, nxlon, nxlat, fraclon, fraclat, lonv, latv, per,      &
+  Nperx, Npery, ilonlat, mindiffLl)
+! Function to search a given value from a coarser version of the data
+
+  USE module_generic
+
+  IMPLICIT NONE
+
+  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
+  INTEGER, INTENT(in)                                    :: dx, dy
+  REAL(r_k), DIMENSION(dx,dy), INTENT(in)                :: ilon, ilat
+  REAL(r_k), DIMENSION(Nperx,Npery), INTENT(in)          :: fraclon, fraclat
+  REAL(r_k), INTENT(in)                                  :: lonv, latv, per
+  REAL(r_k), DIMENSION(2), INTENT(in)                    :: nxlon, nxlat
+  INTEGER, INTENT(in)                                    :: Nperx, Npery
+  INTEGER, DIMENSION(2), INTENT(out)                     :: ilonlat
+  REAL(r_k), INTENT(out)                                 :: mindiffLl
+! Local
+  REAL(r_k), DIMENSION(Nperx,Npery)                      :: difffraclonlat
+  REAL(r_k)                                              :: mindifffracLl
+  INTEGER, DIMENSION(2)                                  :: ilonlatfrac
+  INTEGER                                                :: ixbeg, ixend, iybeg, iyend
+  INTEGER                                                :: fracx, fracy
+  REAL(r_k)                                              :: fraclonv, fraclatv
+  REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: difflonlat, lon, lat
+  CHARACTER(LEN=50)                                      :: fname
+  
+! Variables
+! ilon, ilat: original 2D matrices with the longitudes and the latitudes
+! lonv, latv: longitude and latitude to find
+! nxlon, nxlat: minimum and maximum longitude and latitude of the target lon,lat
+! per: fraction of the whole domain (as percentage)
+! Nper[x/y]: period (as fraction over 1) of the fractions of the original grid to use to explore
+! fraclon, fraclat: longitude and latitude fractional matricies to perform the first guess
+
+  fname = 'CoarselonlatFind'
+
+  IF (lonv < nxlon(1) .OR. lonv > nxlon(2)) THEN
+    PRINT *, TRIM(ErrWarnMsg('err'))
+    PRINT *,'  ' // TRIM(fname) // ': longitude outside data range!!'
+    PRINT *,'    given value:', lonv,' outside (',nxlon(1),' ,',nxlon(2),' )'
+    STOP
+  END IF
+  IF (latv < nxlat(1) .OR. latv > nxlat(2)) THEN
+    PRINT *, TRIM(ErrWarnMsg('err'))
+    PRINT *,'  ' // TRIM(fname) // ': latitude outside data range!!'
+    PRINT *,'    given value:', latv,' outside (',nxlat(1),' ,',nxlat(2),' )'
+    STOP
+  END IF
+
+! Initializing variables
+  ixbeg = 0
+  ixend = 0
+  iybeg = 0
+  iyend = 0
+
+  fracx = int(dx*per)
+  fracy = int(dy*per)
+
+!  PRINT *,'fraclon _______'
+!  PRINT *,fraclon
+
+!  PRINT *,'fraclat _______'
+!  PRINT *,fraclat
+
+! Fraction point
+  difffraclonlat = SQRT((fraclon-lonv)**2. + (fraclat-latv)**2.)
+  mindifffracLl = MINVAL(difffraclonlat)
+  ilonlatfrac = index2DArrayR(difffraclonlat, Nperx, Npery, mindifffracLl)
+
+!  PRINT *, 'mindifffracLl:', mindifffracLl, ' ilonlatfrac:', ilonlatfrac
+!  PRINT *, 'frac lon, lat:', fraclon(ilonlatfrac(1),ilonlatfrac(2)), ' ,',                            &
+!    fraclat(ilonlatfrac(1),ilonlatfrac(2))
+!  PRINT *, 'values lon, lat:', lonv, latv
+      
+! Providing fraction range
+  fraclonv = fraclon(ilonlatfrac(1),ilonlatfrac(2))
+  fraclatv = fraclat(ilonlatfrac(1),ilonlatfrac(2))
+
+  IF (fraclonv >= lonv .AND. fraclatv >= latv) THEN
+    IF (ilonlatfrac(1) > 0) THEN
+      ixbeg = (ilonlatfrac(1)-1)*fracx
+      ixend = ilonlatfrac(1)*fracx+1
+    ELSE
+      ixbeg = 0
+      ixend = fracx+1
+    END IF
+    IF (ilonlatfrac(2) > 0) THEN
+      iybeg = (ilonlatfrac(2)-1)*fracy
+      iyend = ilonlatfrac(2)*fracy+1
+    ELSE
+      iybeg = 0
+      iyend = fracy+1
+    END IF
+  ELSE IF (fraclonv < lonv .AND. fraclatv >= latv) THEN
+    IF (ilonlatfrac(1) < Nperx) THEN
+      IF (ilonlatfrac(1) /= 0) THEN
+        ixbeg = (ilonlatfrac(1)-1)*fracx
+        ixend = ilonlatfrac(1)*fracx+1
+      ELSE
+        ixbeg = 0
+        ixend = fracx+1
+      END IF
+    ELSE
+      ixbeg = Nperx*fracx
+      ixend = dx+1
+    END IF
+    IF (ilonlatfrac(2) > 0) THEN
+      iybeg = (ilonlatfrac(2)-1)*fracy
+      iyend = ilonlatfrac(2)*fracy+1
+    ELSE
+      iybeg = 0
+      iyend = fracy+1
+    END IF    
+  ELSE IF (fraclonv < lonv .AND. fraclatv < latv) THEN
+    IF (ilonlatfrac(1) < Nperx) THEN
+      IF (ilonlatfrac(1) /= 0) THEN
+        ixbeg = (ilonlatfrac(1)-1)*fracx
+        ixend = ilonlatfrac(1)*fracx+1
+      ELSE
+        ixbeg = 0
+        ixend = fracx+1
+      END IF
+    ELSE
+      ixbeg = Nperx*fracx
+      ixend = dx+1
+    ENDIF
+    IF (ilonlatfrac(2) < Npery) THEN
+      IF (ilonlatfrac(2) /= 0) THEN
+        iybeg = (ilonlatfrac(2)-1)*fracy
+        iyend = ilonlatfrac(2)*fracy+1
+      ELSE
+        iybeg = 0
+        iyend = fracy+1
+      END IF
+    ELSE
+      iybeg = Npery*fracy
+      iyend = dy+1
+    END IF
+  ELSE IF (fraclonv >= lonv .AND. fraclatv < latv) THEN
+    IF (ilonlatfrac(1) > 0) THEN
+      ixbeg = (ilonlatfrac(1)-1)*fracx
+      ixend = ilonlatfrac(1)*fracx+1
+    ELSE
+      ixbeg = 0
+      ixend = fracx+1
+    END IF
+    IF (ilonlatfrac(2) < Npery) THEN
+      IF (ilonlatfrac(2) /= 0) THEN
+        iybeg = (ilonlatfrac(2)-1)*fracy
+        iyend = ilonlatfrac(2)*fracy+1
+      ELSE
+        iybeg = 0
+        iyend = fracy+1
+      END IF
+    ELSE
+      iybeg = Npery*fracy
+      iyend = dy+1
+    END IF
+  END IF
+
+  IF (ALLOCATED(lon)) DEALLOCATE(lon)
+  ALLOCATE(lon(ixend-ixbeg+1, iyend-iybeg+1))
+  IF (ALLOCATED(lat)) DEALLOCATE(lat)
+  ALLOCATE(lat(ixend-ixbeg+1, iyend-iybeg+1))
+  IF (ALLOCATED(difflonlat)) DEALLOCATE(difflonlat)
+  ALLOCATE(difflonlat(ixend-ixbeg+1, iyend-iybeg+1))
+
+  lon = ilon(ixbeg:ixend,iybeg:iyend)
+  lat = ilat(ixbeg:ixend,iybeg:iyend)
+
+!  print *,'lon _______'
+!  print *,lon
+!  print *,'lat _______'
+!  print *,lat
+
+! Find point
+  difflonlat = SQRT((lon-lonv)**2. + (lat-latv)**2.)
+  mindiffLl = MINVAL(difflonlat)
+  ilonlat = index2DArrayR(difflonlat, ixend-ixbeg+1, iyend-iybeg+1, mindiffLl)
+
+  ilonlat(1) = ilonlat(1) + ixbeg
+  ilonlat(2) = ilonlat(2) + iybeg
+
+!  PRINT *,'mindiffLl:', mindiffLl, ' ilatlon:', ilatlon
+!  PRINT *,'lon, lat:', lon(ilonlat(1),ilonlat(2)), ' ,', lat(ilonlat(1),ilonlat(2))
+
+  RETURN
+
+END SUBROUTINE CoarselonlatFind
+
+SUBROUTINE CoarseInterpolate(dimx, dimy, projlon, projlat, Ninpts, lonvs, latvs, percen, mindiff,     &
+  ivar, newvar, newvarin, newvarinpt, newvarindiff)
+! Subroutine which finds the closest grid point within a projection throughout a first guest 
+!   approche from percentages of the whole domain
+
+  USE module_generic
+
+  IMPLICIT NONE
+
+  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
+  INTEGER, INTENT(in)                                    :: dimx, dimy
+  REAL(r_k), DIMENSION(dimx,dimy), INTENT(in)            :: projlon, projlat
+  INTEGER, INTENT(in)                                    :: Ninpts
+  REAL(r_k), DIMENSION(Ninpts), INTENT(in)               :: ivar, lonvs, latvs
+  REAL(r_k)                                              :: mindiff, percen
+  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
+! Local
+  INTEGER                                                :: iv
+  INTEGER, DIMENSION(2)                                  :: ilonlat
+  REAL(r_k)                                              :: mindiffLl
+  INTEGER                                                :: Ninpts1
+  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
+  REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: fractionlon, fractionlat
+  INTEGER                                                :: fracdx, fracdy
+  CHARACTER(LEN=50)                                      :: fname
+
+!!!!!!! Variables
+! dimx, dimy: dimension length of the target interpolation
+! proj[lon/lat]: longitudes and latitudes of the target interpolation
+! Ninpts: number of points to interpolate
+! [lon/lat]vs: longitudes and latitudes of the points to interpolate
+! mindiff: minimal accepted distance to the target point
+! percen: size (as percentage of the total domain) of the first guess portions to provide the first gues
+! ivar: values to localize in the target projection
+! newvar: localisation of the [lon/lat]vs point in the target projection
+! newvarin: number of point from the input data
+! 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
+! ncid: netCDF output file id
+
+  fname = 'CoarseInterpolate'
+  Ninpts1 = Ninpts/100
+
+  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
+  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
+
+  fracdx = INT(dimx*percen)
+  fracdy = INT(dimy*percen)
+
+  IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon)
+  ALLOCATE(fractionlon(fracdx, fracdy))
+  IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat)
+  ALLOCATE(fractionlat(fracdx, fracdy))
+
+  fractionlon = projlon(::fracdx,::fracdy)
+  fractionlat = projlat(::fracdx,::fracdy)
+
+  DO iv=1,Ninpts
+    IF (newvarinpt(iv) == 0) THEN
+      CALL CoarselonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, fractionlon,        &
+        fractionlat, lonvs(iv), latvs(iv), percen, fracdx, fracdy, ilonlat, 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
+          newvarinpt(iv) = 1
+          newvarindiff(iv) = mindiffLl
+          PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv),        &
+            ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:',           &
+            newvarindiff(iv), ' point:',ilonlat    
+        ELSE
+          PRINT *,TRIM(ErrWarnMsg('err'))
+          PRINT *,'  ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv),        &
+            ' not relocated !!'
+          PRINT *,'    mindiffl:', mindiffLl, ' ilon:', ilonlat(1), ' ilat:', ilonlat(2)
+          STOP
+        END IF
+
+!        IF (MOD(iv,Ninpts1) == 0) newnc.sync()
+      ELSE
+        PRINT *,TRIM(ErrWarnMsg('err'))
+        PRINT *,'  ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv),   &
+          ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ',       &
+          mindiff, ' !!'
+        PRINT *,'    found minimum difference:', mindiffLl
+        STOP
+      END IF
+    END IF
+  END DO
+
+END SUBROUTINE CoarseInterpolate
+
+END MODULE module_ForInterpolate
Index: /trunk/tools/module_generic.F90
===================================================================
--- /trunk/tools/module_generic.F90	(revision 715)
+++ /trunk/tools/module_generic.F90	(revision 715)
@@ -0,0 +1,59 @@
+MODULE module_generic
+! Module with generic functions
+
+  CONTAINS
+
+  FUNCTION Index2DArrayR(array2D, d1, d2, val)
+! Function to provide the first index of a given value inside a 2D real array
+
+    IMPLICIT NONE
+
+    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
+    INTEGER, INTENT(in)                                  :: d1, d2
+    REAL(r_k), INTENT(in)                                :: val
+    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: array2D
+    INTEGER, DIMENSION(2)                                :: Index2DArrayR
+
+! Local
+    INTEGER                                              :: i, j
+    CHARACTER(LEN=50)                                    :: fname
+
+    fname = 'Index2DArrayR'
+
+    Index2DArrayR = -1
+
+    DO i=1,d1
+      DO j=1,d2
+        IF (array2d(i,j) == val) THEN
+          Index2DArrayR(1) = i
+          Index2DArrayR(2) = j
+          EXIT
+        END IF
+      END DO
+    END DO
+
+  END FUNCTION Index2DArrayR
+
+  CHARACTER(LEN=50) FUNCTION ErrWarnMsg(msg)
+! Function to print error/warning message
+
+    IMPLICIT NONE
+
+    CHARACTER(LEN=3), INTENT(in)                         :: msg
+! Local
+    CHARACTER(LEN=50)                                    :: fname
+
+    fname = 'ErrWarnMsg'
+
+    IF (msg == 'err') THEN
+      ErrWarnMsg = 'ERROR -- error -- ERROR -- error'
+    ELSE IF (msg == 'wrn') THEN
+      ErrWarnMsg = 'WARNING -- warning -- WARNING -- warning'
+    ELSE
+      PRINT *,'ERROR -- error -- ERROR -- error'
+      PRINT *,'  ' // TRIM(fname) // ": '" // TRIM(msg) // "' does not exist!!" 
+      STOP
+    END IF
+  END FUNCTION ErrWarnMsg
+
+END MODULE module_generic
