Index: unk/tools/interpolate.F90
===================================================================
--- /trunk/tools/interpolate.F90	(revision 1605)
+++ 	(revision )
@@ -1,296 +1,0 @@
-PROGRAM interpolate
-! Program to interpolate values from a giving projection
-! To be included in a python
-! f2py -m ForInterpolate --f90exec=/usr/bin/gfortran-4.7 -c interpolate.F90 module_generic.F90
-
-  IMPLICIT NONE
-
-  CHARACTER(LEN=50)                                      :: main, ErrWarnMsg
-
-  main='interpolate'
-
-END PROGRAM interpolate
-
-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
-
-  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, ncid)
-! 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
-  INTEGER, INTENT(in)                                    :: ncid
-  REAL(r_k), DIMENSION(dimx,dimy), INTENT(inout)         :: newvar
-  INTEGER, DIMENSION(dimx,dimy), INTENT(inout)           :: newvarin
-  INTEGER, DIMENSION(Ninpts), INTENT(inout)              :: newvarinpt
-  REAL(r_k), DIMENSION(Ninpts), INTENT(inout)            :: 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
Index: /trunk/tools/interpolate.f90
===================================================================
--- /trunk/tools/interpolate.f90	(revision 1606)
+++ /trunk/tools/interpolate.f90	(revision 1606)
@@ -0,0 +1,296 @@
+PROGRAM interpolate
+! Program to interpolate values from a giving projection
+! To be included in a python
+! f2py -m ForInterpolate --f90exec=/usr/bin/gfortran-4.7 -c interpolate.F90 module_generic.F90
+
+  IMPLICIT NONE
+
+  CHARACTER(LEN=50)                                      :: main, ErrWarnMsg
+
+  main='interpolate'
+
+END PROGRAM interpolate
+
+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
+
+  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, ncid)
+! 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
+  INTEGER, INTENT(in)                                    :: ncid
+  REAL(r_k), DIMENSION(dimx,dimy), INTENT(inout)         :: newvar
+  INTEGER, DIMENSION(dimx,dimy), INTENT(inout)           :: newvarin
+  INTEGER, DIMENSION(Ninpts), INTENT(inout)              :: newvarinpt
+  REAL(r_k), DIMENSION(Ninpts), INTENT(inout)            :: 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
Index: unk/tools/module_ForDiagnostics.F90
===================================================================
--- /trunk/tools/module_ForDiagnostics.F90	(revision 1605)
+++ 	(revision )
@@ -1,462 +1,0 @@
-!! Fortran version of different diagnostics
-! L. Fita. LMD May 2016
-! gfortran module_generic.o module_ForDiagnosticsVars.o -c module_ForDiagnostics.F90
-!
-! f2py -m module_ForDiagnostics --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForDiagnosticsVars.F90 module_ForDiagnostics.F90
-
-MODULE module_ForDiagnostics
-
-  USE module_generic
-  USE module_ForDiagnosticsVars
-
-  CONTAINS
-
-!!!!!!! Calculations
-! compute_cllmh4D3: Computation of low, medium and high cloudiness from a 4D CLDFRA and pressure  being 3rd dimension the z-dim 
-! compute_cllmh3D3: Computation of low, medium and high cloudiness from a 3D CLDFRA and pressure  being 3rd dimension the z-dim 
-! compute_cllmh: Computation of low, medium and high cloudiness
-! compute_clt4D3: Computation of total cloudiness from a 4D CLDFRA being 3rd dimension the z-dim
-! compute_clt3D3: Computation of total cloudiness from a 3D CLDFRA being 3rd dimension the z-dim
-! compute_clt: Computation of total cloudiness
-
-!!!
-! Calculations
-!!!
-
-  SUBROUTINE compute_cllmh4D2(cldfra4D, pres4D, cllmh4D2, d1, d2, d3, d4)
-! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ from a 4D CLDFRA and pressure 
-!   where zdim is the 2nd dimension (thus, cldfra4D(d1,d2,d3,d4) --> cllmh(3,d1,d3,d4) 1: low, 2: medium, 3: high
-! It should be properly done via an 'INTERFACE', but...
-
-    IMPLICIT NONE
-
-    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
-    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
-    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: cldfra4D, pres4D
-    REAL(r_k), DIMENSION(3,d1,d3,d4), INTENT(out)        :: cllmh4D2
-
-! Local
-    INTEGER                                              :: i,j,k, zdim, Ndim
-    CHARACTER(LEN=50)                                    :: fname
-
-!!!!!!! Variables
-! cldfra4D: 4D cloud fraction values [1]
-! pres4D: 4D pressure values [Pa]
-! Ndim: number of dimensions of the input data
-! d[1-4]: dimensions of 'cldfra4D'
-! zdim: number of the vertical-dimension within the matrix
-! cltlmh4D2: low, medium, high cloudiness for the 4D cldfra and d2 being 'zdim'
-
-    fname = 'compute_cllmh4D2'
-    zdim = 2
-    Ndim = 4
-
-    DO i=1, d1
-      DO j=1, d3
-        DO k=1, d4
-          cllmh4D2(:,i,j,k) = var_cllmh(cldfra4D(i,:,j,k), pres4D(i,:,j,k), d2)
-        END DO
-      END DO
-    END DO
-    
-    RETURN 
-
-  END SUBROUTINE compute_cllmh4D2
-
-  SUBROUTINE compute_cllmh3D1(cldfra3D, pres3D, cllmh3D1, d1, d2, d3)
-! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ from a 3D CLDFRA and pressure 
-!   where zdim is the 1st dimension (thus, cldfra3D(d1,d2,d3) --> cllmh(3,d2,d3) 1: low, 2: medium, 3: high
-! It should be properly done via an 'INTERFACE', but...
-
-    IMPLICIT NONE
-
-    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
-    INTEGER, INTENT(in)                                  :: d1, d2, d3
-    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: cldfra3D, pres3D
-    REAL(r_k), DIMENSION(3,d2,d3), INTENT(out)           :: cllmh3D1
-
-! Local
-    INTEGER                                              :: i,j,k, zdim, Ndim
-    CHARACTER(LEN=50)                                    :: fname
-
-!!!!!!! Variables
-! cldfra3D: 3D cloud fraction values [1]
-! pres3D: 3D pressure values [Pa]
-! Ndim: number of dimensions of the input data
-! d[1-3]: dimensions of 'cldfra3D'
-! zdim: number of the vertical-dimension within the matrix
-! cltlmh3D1: low, medium, high cloudiness for the 3D cldfra and d1 being 'zdim'
-
-    fname = 'compute_cllmh3D1'
-    zdim = 1
-    Ndim = 3
-
-    DO i=1, d1
-      DO j=1, d2
-        cllmh3D1(:,i,j) = var_cllmh(cldfra3D(:,i,j), pres3D(:,i,j), d1)
-      END DO
-    END DO
-    
-    RETURN 
-
-  END SUBROUTINE compute_cllmh3D1
-
-  SUBROUTINE compute_cllmh(cldfra1D, cldfra2D, cldfra3D, cldfra4D, pres1D, pres2D, pres3D, pres4D,    &
-    Ndim, zdim, cllmh1D, cllmh2D1, cllmh2D2, cllmh3D1, cllmh3D2, cllmh3D3, cllmh4D1, cllmh4D2,        &
-    cllmh4D3, cllmh4D4, d1, d2, d3, d4)
-! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ
-
-    IMPLICIT NONE
-
-    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
-    INTEGER, INTENT(in)                                  :: Ndim, d1, d2, d3, d4, zdim
-    REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in)       :: cldfra1D, pres1D
-    REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(in)    :: cldfra2D, pres2D
-    REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL, INTENT(in) :: cldfra3D, pres3D
-    REAL(r_k), DIMENSION(d1,d2,d3,d4), OPTIONAL,                                                      &
-      INTENT(in)                                         :: cldfra4D, pres4D
-    REAL(r_k), DIMENSION(3), OPTIONAL, INTENT(out)       :: cllmh1D
-    REAL(r_k), DIMENSION(d1,3), OPTIONAL, INTENT(out)    :: cllmh2D1
-    REAL(r_k), DIMENSION(d2,3), OPTIONAL, INTENT(out)    :: cllmh2D2
-    REAL(r_k), DIMENSION(d2,d3,3), OPTIONAL, INTENT(out) :: cllmh3D1
-    REAL(r_k), DIMENSION(d1,d3,3), OPTIONAL, INTENT(out) :: cllmh3D2
-    REAL(r_k), DIMENSION(d1,d2,3), OPTIONAL, INTENT(out) :: cllmh3D3
-    REAL(r_k), DIMENSION(d2,d3,d4,3), OPTIONAL,                                                       &
-      INTENT(out)                                        :: cllmh4D1
-    REAL(r_k), DIMENSION(d1,d3,d4,3), OPTIONAL,                                                       &
-      INTENT(out)                                        :: cllmh4D2
-    REAL(r_k), DIMENSION(d1,d2,d4,3), OPTIONAL,                                                       &
-      INTENT(out)                                        :: cllmh4D3
-    REAL(r_k), DIMENSION(d1,d2,d3,3), OPTIONAL,                                                       &
-      INTENT(out)                                        :: cllmh4D4
-
-! Local
-    INTEGER                                              :: i,j,k
-    CHARACTER(LEN=50)                                    :: fname
-
-!!!!!!! Variables
-! cldfra[1-4]D: cloud fraction values [1]
-! pres[1-4]D: pressure values [Pa]
-! Ndim: number of dimensions of the input data
-! d[1-4]: dimensions of 'cldfra'
-! zdim: number of the vertical-dimension within the matrix
-! cllmh1D: low, medium and high cloudiness for the 1D cldfra
-! cllmh2D1: low, medium and high cloudiness for the 2D cldfra and d1 being 'zdim'
-! cllmh2D2: low, medium and high cloudiness for the 2D cldfra and d2 being 'zdim'
-! cllmh3D1: low, medium and high cloudiness for the 3D cldfra and d1 being 'zdim'
-! cllmh3D2: low, medium and high cloudiness for the 3D cldfra and d2 being 'zdim'
-! cllmh3D3: low, medium and high cloudiness for the 3D cldfra and d3 being 'zdim'
-! cllmh4D1: low, medium and high cloudiness for the 4D cldfra and d1 being 'zdim'
-! cllmh4D2: low, medium and high cloudiness for the 4D cldfra and d2 being 'zdim'
-! cllmh4D3: low, medium and high cloudiness for the 4D cldfra and d3 being 'zdim'
-! cllmh4D4: low, medium and high cloudiness for the 4D cldfra and d4 being 'zdim'
-
-    fname = 'compute_cllmh'
-
-    SELECT CASE (Ndim)
-      CASE (1)
-        cllmh1D = var_cllmh(cldfra1D, pres1D, d1)
-      CASE (2)
-        IF (zdim == 1) THEN
-          DO i=1, d2
-            cllmh2D1(i,:) = var_cllmh(cldfra2D(:,i), pres2D(:,i), d1)
-          END DO
-        ELSE IF (zdim == 2) THEN
-          DO i=1, d1
-            cllmh2D2(i,:) = var_cllmh(cldfra2D(:,i), pres2D(i,:), d2)
-          END DO
-        ELSE
-          PRINT *,TRIM(ErrWarnMsg('err'))
-          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
-          PRINT *,'    accepted values: 1,2'
-          STOP
-        END IF
-      CASE (3)
-        IF (zdim == 1) THEN
-          DO i=1, d2
-            DO j=1, d3
-              cllmh3D1(i,j,:) = var_cllmh(cldfra3D(:,i,j), pres3D(:,i,j), d1)
-            END DO
-          END DO
-        ELSE IF (zdim == 2) THEN
-          DO i=1, d1
-            DO j=1, d3
-              cllmh3D2(i,j,:) = var_cllmh(cldfra3D(i,:,j), pres3D(i,:,j), d2)
-            END DO
-          END DO
-        ELSE IF (zdim == 3) THEN
-          DO i=1, d1
-            DO j=1, d2
-              cllmh3D3(i,j,:) = var_cllmh(cldfra3D(i,j,:), pres3D(i,j,:), d3)
-            END DO
-          END DO
-        ELSE
-          PRINT *,TRIM(ErrWarnMsg('err'))
-          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
-          PRINT *,'    accepted values: 1,2,3'
-          STOP
-        END IF
-      CASE (4)
-        IF (zdim == 1) THEN
-          DO i=1, d2
-            DO j=1, d3
-              DO k=1, d4
-                cllmh4D1(i,j,k,:) = var_cllmh(cldfra4D(:,i,j,k), pres4D(:,i,j,k), d1)
-              END DO
-            END DO
-          END DO
-        ELSE IF (zdim == 2) THEN
-          DO i=1, d1
-            DO j=1, d3
-              DO k=1, d4
-                cllmh4D2(i,j,k,:) = var_cllmh(cldfra4D(i,:,j,k), pres4D(i,:,j,k), d2)
-              END DO
-            END DO
-          END DO
-        ELSE IF (zdim == 3) THEN
-          DO i=1, d2
-            DO j=1, d3
-              DO k=1, d4
-                cllmh4D3(i,j,k,:) = var_cllmh(cldfra4D(i,j,:,k), pres4D(i,j,:,k), d3)
-              END DO
-            END DO
-          END DO
-        ELSE IF (zdim == 4) THEN
-          DO i=1, d1
-            DO j=1, d2
-              DO k=1, d3
-                cllmh4D4(i,j,k,:) = var_cllmh(cldfra4D(i,j,k,:), pres4D(i,j,k,:), d4)
-              END DO
-            END DO
-          END DO
-        ELSE
-          PRINT *,TRIM(ErrWarnMsg('err'))
-          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
-          PRINT *,'    accepted values: 1,2,3,4'
-          STOP
-        END IF
-      CASE DEFAULT
-        PRINT *,TRIM(ErrWarnMsg('err'))
-        PRINT *,'  ' // TRIM(fname) // ': Ndim:', Ndim,' not ready !!'
-        STOP
-      END SELECT
-
-    RETURN 
-
-  END SUBROUTINE compute_cllmh
-
-  SUBROUTINE compute_clt4D2(cldfra4D, clt4D2, d1, d2, d3, d4)
-! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ from a 4D CLDFRA
-!   where zdim is the 2nd dimension (thus, cldfra4D(d1,d2,d3,d4) --> clt(d1,d3,d4)
-! It should be properly done via an 'INTERFACE', but...
-
-    IMPLICIT NONE
-
-    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
-    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
-    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: cldfra4D
-    REAL(r_k), DIMENSION(d1,d3,d4), INTENT(out)          :: clt4D2
-
-! Local
-    INTEGER                                              :: i,j,k, zdim, Ndim
-    CHARACTER(LEN=50)                                    :: fname
-
-!!!!!!! Variables
-! cldfra4D: 4D cloud fraction values [1]
-! Ndim: number of dimensions of the input data
-! d[1-4]: dimensions of 'cldfra4D'
-! zdim: number of the vertical-dimension within the matrix
-! clt4D2: total cloudiness for the 4D cldfra and d2 being 'zdim'
-
-    fname = 'compute_clt4D2'
-    zdim = 2
-    Ndim = 4
-
-    DO i=1, d1
-      DO j=1, d3
-        DO k=1, d4
-          clt4D2(i,j,k) = var_clt(cldfra4D(i,:,j,k), d2)
-        END DO
-      END DO
-    END DO
-    
-    RETURN 
-
-  END SUBROUTINE compute_clt4D2
-
-  SUBROUTINE compute_clt3D1(cldfra3D, clt3D1, d1, d2, d3)
-! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ from a 3D CLDFRA
-!   where zdim is the 1st dimension (thus, cldfra4D(d1,d2,d3) --> clt(d2,d3)
-! It should be properly done via an 'INTERFACE', but...
-
-    IMPLICIT NONE
-
-    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
-    INTEGER, INTENT(in)                                  :: d1, d2, d3
-    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: cldfra3D
-    REAL(r_k), DIMENSION(d2,d3), INTENT(out)             :: clt3D1
-
-! Local
-    INTEGER                                              :: i,j,k, zdim, Ndim
-    CHARACTER(LEN=50)                                    :: fname
-
-!!!!!!! Variables
-! cldfra3D: 3D cloud fraction values [1]
-! Ndim: number of dimensions of the input data
-! d[1-3]: dimensions of 'cldfra3D'
-! zdim: number of the vertical-dimension within the matrix
-! clt3D1: total cloudiness for the 3D cldfra and d1 being 'zdim'
-
-    fname = 'compute_clt3D1'
-    zdim = 1
-    Ndim = 3
-
-    DO i=1, d2
-      DO j=1, d3
-        clt3D1(i,j) = var_clt(cldfra3D(:,i,j), d1)
-      END DO
-    END DO
-    
-    RETURN 
-
-  END SUBROUTINE compute_clt3D1
-
-  SUBROUTINE compute_clt(cldfra1D, cldfra2D, cldfra3D, cldfra4D, Ndim, zdim, clt1D, clt2D1, clt2D2,   &
-    clt3D1, clt3D2, clt3D3, clt4D1, clt4D2, clt4D3, clt4D4, d1, d2, d3, d4)
-! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ
-
-    IMPLICIT NONE
-
-    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
-    INTEGER, INTENT(in)                                  :: Ndim, d1, d2, d3, d4, zdim
-    REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in)       :: cldfra1D
-    REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(in)    :: cldfra2D
-    REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL, INTENT(in) :: cldfra3D
-    REAL(r_k), DIMENSION(d1,d2,d3,d4), OPTIONAL,                                                      &
-      INTENT(in)                                         :: cldfra4D
-    REAL(r_k), OPTIONAL, INTENT(out)                     :: clt1D
-    REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(out)      :: clt2D1
-    REAL(r_k), DIMENSION(d2), OPTIONAL, INTENT(out)      :: clt2D2
-    REAL(r_k), DIMENSION(d2,d3), OPTIONAL, INTENT(out)   :: clt3D1
-    REAL(r_k), DIMENSION(d1,d3), OPTIONAL, INTENT(out)   :: clt3D2
-    REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(out)   :: clt3D3
-    REAL(r_k), DIMENSION(d2,d3,d4), OPTIONAL,INTENT(out) :: clt4D1
-    REAL(r_k), DIMENSION(d1,d3,d4), OPTIONAL,INTENT(out) :: clt4D2
-    REAL(r_k), DIMENSION(d1,d2,d4), OPTIONAL,INTENT(out) :: clt4D3
-    REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL,INTENT(out) :: clt4D4
-
-! Local
-    INTEGER                                              :: i,j,k
-    CHARACTER(LEN=50)                                    :: fname
-
-!!!!!!! Variables
-! cldfra[1-4]D: cloud fraction values [1]
-! Ndim: number of dimensions of the input data
-! d[1-4]: dimensions of 'cldfra'
-! zdim: number of the vertical-dimension within the matrix
-! clt1D: total cloudiness for the 1D cldfra
-! clt2D1: total cloudiness for the 2D cldfra and d1 being 'zdim'
-! clt2D2: total cloudiness for the 2D cldfra and d2 being 'zdim'
-! clt3D1: total cloudiness for the 3D cldfra and d1 being 'zdim'
-! clt3D2: total cloudiness for the 3D cldfra and d2 being 'zdim'
-! clt3D3: total cloudiness for the 3D cldfra and d3 being 'zdim'
-! clt4D1: total cloudiness for the 4D cldfra and d1 being 'zdim'
-! clt4D2: total cloudiness for the 4D cldfra and d2 being 'zdim'
-! clt4D3: total cloudiness for the 4D cldfra and d3 being 'zdim'
-! clt4D4: total cloudiness for the 4D cldfra and d4 being 'zdim'
-
-    fname = 'compute_clt'
-
-    SELECT CASE (Ndim)
-      CASE (1)
-        clt1D = var_clt(cldfra1D, d1)
-      CASE (2)
-        IF (zdim == 1) THEN
-          DO i=1, d2
-            clt2D1(i) = var_clt(cldfra2D(:,i), d1)
-          END DO
-        ELSE IF (zdim == 2) THEN
-          DO i=1, d1
-            clt2D2(i) = var_clt(cldfra2D(:,i), d2)
-          END DO
-        ELSE
-          PRINT *,TRIM(ErrWarnMsg('err'))
-          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
-          PRINT *,'    accepted values: 1,2'
-          STOP
-        END IF
-      CASE (3)
-        IF (zdim == 1) THEN
-          DO i=1, d2
-            DO j=1, d3
-              clt3D1(i,j) = var_clt(cldfra3D(:,i,j), d1)
-            END DO
-          END DO
-        ELSE IF (zdim == 2) THEN
-          DO i=1, d1
-            DO j=1, d3
-              clt3D2(i,j) = var_clt(cldfra3D(i,:,j), d2)
-            END DO
-          END DO
-        ELSE IF (zdim == 3) THEN
-          DO i=1, d1
-            DO j=1, d2
-              clt3D3(i,j) = var_clt(cldfra3D(i,j,:), d3)
-            END DO
-          END DO
-        ELSE
-          PRINT *,TRIM(ErrWarnMsg('err'))
-          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
-          PRINT *,'    accepted values: 1,2,3'
-          STOP
-        END IF
-      CASE (4)
-        IF (zdim == 1) THEN
-          DO i=1, d2
-            DO j=1, d3
-              DO k=1, d4
-                clt4D1(i,j,k) = var_clt(cldfra4D(:,i,j,k), d1)
-              END DO
-            END DO
-          END DO
-        ELSE IF (zdim == 2) THEN
-          DO i=1, d1
-            DO j=1, d3
-              DO k=1, d4
-                clt4D2(i,j,k) = var_clt(cldfra4D(i,:,j,k), d2)
-              END DO
-            END DO
-          END DO
-        ELSE IF (zdim == 3) THEN
-          DO i=1, d2
-            DO j=1, d3
-              DO k=1, d4
-                clt4D3(i,j,k) = var_clt(cldfra4D(i,j,:,k), d3)
-              END DO
-            END DO
-          END DO
-        ELSE IF (zdim == 4) THEN
-          DO i=1, d1
-            DO j=1, d2
-              DO k=1, d3
-                clt4D4(i,j,k) = var_clt(cldfra4D(i,j,k,:), d4)
-              END DO
-            END DO
-          END DO
-        ELSE
-          PRINT *,TRIM(ErrWarnMsg('err'))
-          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
-          PRINT *,'    accepted values: 1,2,3,4'
-          STOP
-        END IF
-      CASE DEFAULT
-        PRINT *,TRIM(ErrWarnMsg('err'))
-        PRINT *,'  ' // TRIM(fname) // ': Ndim:', Ndim,' not ready !!'
-        STOP
-      END SELECT
-
-    RETURN 
-
-  END SUBROUTINE compute_clt
-
-END MODULE module_ForDiagnostics
Index: /trunk/tools/module_ForDiagnostics.f90
===================================================================
--- /trunk/tools/module_ForDiagnostics.f90	(revision 1606)
+++ /trunk/tools/module_ForDiagnostics.f90	(revision 1606)
@@ -0,0 +1,462 @@
+!! Fortran version of different diagnostics
+! L. Fita. LMD May 2016
+! gfortran module_generic.o module_ForDiagnosticsVars.o -c module_ForDiagnostics.F90
+!
+! f2py -m module_ForDiagnostics --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForDiagnosticsVars.F90 module_ForDiagnostics.F90
+
+MODULE module_ForDiagnostics
+
+  USE module_generic
+  USE module_ForDiagnosticsVars
+
+  CONTAINS
+
+!!!!!!! Calculations
+! compute_cllmh4D3: Computation of low, medium and high cloudiness from a 4D CLDFRA and pressure  being 3rd dimension the z-dim 
+! compute_cllmh3D3: Computation of low, medium and high cloudiness from a 3D CLDFRA and pressure  being 3rd dimension the z-dim 
+! compute_cllmh: Computation of low, medium and high cloudiness
+! compute_clt4D3: Computation of total cloudiness from a 4D CLDFRA being 3rd dimension the z-dim
+! compute_clt3D3: Computation of total cloudiness from a 3D CLDFRA being 3rd dimension the z-dim
+! compute_clt: Computation of total cloudiness
+
+!!!
+! Calculations
+!!!
+
+  SUBROUTINE compute_cllmh4D2(cldfra4D, pres4D, cllmh4D2, d1, d2, d3, d4)
+! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ from a 4D CLDFRA and pressure 
+!   where zdim is the 2nd dimension (thus, cldfra4D(d1,d2,d3,d4) --> cllmh(3,d1,d3,d4) 1: low, 2: medium, 3: high
+! It should be properly done via an 'INTERFACE', but...
+
+    IMPLICIT NONE
+
+    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
+    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
+    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: cldfra4D, pres4D
+    REAL(r_k), DIMENSION(3,d1,d3,d4), INTENT(out)        :: cllmh4D2
+
+! Local
+    INTEGER                                              :: i,j,k, zdim, Ndim
+    CHARACTER(LEN=50)                                    :: fname
+
+!!!!!!! Variables
+! cldfra4D: 4D cloud fraction values [1]
+! pres4D: 4D pressure values [Pa]
+! Ndim: number of dimensions of the input data
+! d[1-4]: dimensions of 'cldfra4D'
+! zdim: number of the vertical-dimension within the matrix
+! cltlmh4D2: low, medium, high cloudiness for the 4D cldfra and d2 being 'zdim'
+
+    fname = 'compute_cllmh4D2'
+    zdim = 2
+    Ndim = 4
+
+    DO i=1, d1
+      DO j=1, d3
+        DO k=1, d4
+          cllmh4D2(:,i,j,k) = var_cllmh(cldfra4D(i,:,j,k), pres4D(i,:,j,k), d2)
+        END DO
+      END DO
+    END DO
+    
+    RETURN 
+
+  END SUBROUTINE compute_cllmh4D2
+
+  SUBROUTINE compute_cllmh3D1(cldfra3D, pres3D, cllmh3D1, d1, d2, d3)
+! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ from a 3D CLDFRA and pressure 
+!   where zdim is the 1st dimension (thus, cldfra3D(d1,d2,d3) --> cllmh(3,d2,d3) 1: low, 2: medium, 3: high
+! It should be properly done via an 'INTERFACE', but...
+
+    IMPLICIT NONE
+
+    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
+    INTEGER, INTENT(in)                                  :: d1, d2, d3
+    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: cldfra3D, pres3D
+    REAL(r_k), DIMENSION(3,d2,d3), INTENT(out)           :: cllmh3D1
+
+! Local
+    INTEGER                                              :: i,j,k, zdim, Ndim
+    CHARACTER(LEN=50)                                    :: fname
+
+!!!!!!! Variables
+! cldfra3D: 3D cloud fraction values [1]
+! pres3D: 3D pressure values [Pa]
+! Ndim: number of dimensions of the input data
+! d[1-3]: dimensions of 'cldfra3D'
+! zdim: number of the vertical-dimension within the matrix
+! cltlmh3D1: low, medium, high cloudiness for the 3D cldfra and d1 being 'zdim'
+
+    fname = 'compute_cllmh3D1'
+    zdim = 1
+    Ndim = 3
+
+    DO i=1, d1
+      DO j=1, d2
+        cllmh3D1(:,i,j) = var_cllmh(cldfra3D(:,i,j), pres3D(:,i,j), d1)
+      END DO
+    END DO
+    
+    RETURN 
+
+  END SUBROUTINE compute_cllmh3D1
+
+  SUBROUTINE compute_cllmh(cldfra1D, cldfra2D, cldfra3D, cldfra4D, pres1D, pres2D, pres3D, pres4D,    &
+    Ndim, zdim, cllmh1D, cllmh2D1, cllmh2D2, cllmh3D1, cllmh3D2, cllmh3D3, cllmh4D1, cllmh4D2,        &
+    cllmh4D3, cllmh4D4, d1, d2, d3, d4)
+! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ
+
+    IMPLICIT NONE
+
+    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
+    INTEGER, INTENT(in)                                  :: Ndim, d1, d2, d3, d4, zdim
+    REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in)       :: cldfra1D, pres1D
+    REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(in)    :: cldfra2D, pres2D
+    REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL, INTENT(in) :: cldfra3D, pres3D
+    REAL(r_k), DIMENSION(d1,d2,d3,d4), OPTIONAL,                                                      &
+      INTENT(in)                                         :: cldfra4D, pres4D
+    REAL(r_k), DIMENSION(3), OPTIONAL, INTENT(out)       :: cllmh1D
+    REAL(r_k), DIMENSION(d1,3), OPTIONAL, INTENT(out)    :: cllmh2D1
+    REAL(r_k), DIMENSION(d2,3), OPTIONAL, INTENT(out)    :: cllmh2D2
+    REAL(r_k), DIMENSION(d2,d3,3), OPTIONAL, INTENT(out) :: cllmh3D1
+    REAL(r_k), DIMENSION(d1,d3,3), OPTIONAL, INTENT(out) :: cllmh3D2
+    REAL(r_k), DIMENSION(d1,d2,3), OPTIONAL, INTENT(out) :: cllmh3D3
+    REAL(r_k), DIMENSION(d2,d3,d4,3), OPTIONAL,                                                       &
+      INTENT(out)                                        :: cllmh4D1
+    REAL(r_k), DIMENSION(d1,d3,d4,3), OPTIONAL,                                                       &
+      INTENT(out)                                        :: cllmh4D2
+    REAL(r_k), DIMENSION(d1,d2,d4,3), OPTIONAL,                                                       &
+      INTENT(out)                                        :: cllmh4D3
+    REAL(r_k), DIMENSION(d1,d2,d3,3), OPTIONAL,                                                       &
+      INTENT(out)                                        :: cllmh4D4
+
+! Local
+    INTEGER                                              :: i,j,k
+    CHARACTER(LEN=50)                                    :: fname
+
+!!!!!!! Variables
+! cldfra[1-4]D: cloud fraction values [1]
+! pres[1-4]D: pressure values [Pa]
+! Ndim: number of dimensions of the input data
+! d[1-4]: dimensions of 'cldfra'
+! zdim: number of the vertical-dimension within the matrix
+! cllmh1D: low, medium and high cloudiness for the 1D cldfra
+! cllmh2D1: low, medium and high cloudiness for the 2D cldfra and d1 being 'zdim'
+! cllmh2D2: low, medium and high cloudiness for the 2D cldfra and d2 being 'zdim'
+! cllmh3D1: low, medium and high cloudiness for the 3D cldfra and d1 being 'zdim'
+! cllmh3D2: low, medium and high cloudiness for the 3D cldfra and d2 being 'zdim'
+! cllmh3D3: low, medium and high cloudiness for the 3D cldfra and d3 being 'zdim'
+! cllmh4D1: low, medium and high cloudiness for the 4D cldfra and d1 being 'zdim'
+! cllmh4D2: low, medium and high cloudiness for the 4D cldfra and d2 being 'zdim'
+! cllmh4D3: low, medium and high cloudiness for the 4D cldfra and d3 being 'zdim'
+! cllmh4D4: low, medium and high cloudiness for the 4D cldfra and d4 being 'zdim'
+
+    fname = 'compute_cllmh'
+
+    SELECT CASE (Ndim)
+      CASE (1)
+        cllmh1D = var_cllmh(cldfra1D, pres1D, d1)
+      CASE (2)
+        IF (zdim == 1) THEN
+          DO i=1, d2
+            cllmh2D1(i,:) = var_cllmh(cldfra2D(:,i), pres2D(:,i), d1)
+          END DO
+        ELSE IF (zdim == 2) THEN
+          DO i=1, d1
+            cllmh2D2(i,:) = var_cllmh(cldfra2D(:,i), pres2D(i,:), d2)
+          END DO
+        ELSE
+          PRINT *,TRIM(ErrWarnMsg('err'))
+          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
+          PRINT *,'    accepted values: 1,2'
+          STOP
+        END IF
+      CASE (3)
+        IF (zdim == 1) THEN
+          DO i=1, d2
+            DO j=1, d3
+              cllmh3D1(i,j,:) = var_cllmh(cldfra3D(:,i,j), pres3D(:,i,j), d1)
+            END DO
+          END DO
+        ELSE IF (zdim == 2) THEN
+          DO i=1, d1
+            DO j=1, d3
+              cllmh3D2(i,j,:) = var_cllmh(cldfra3D(i,:,j), pres3D(i,:,j), d2)
+            END DO
+          END DO
+        ELSE IF (zdim == 3) THEN
+          DO i=1, d1
+            DO j=1, d2
+              cllmh3D3(i,j,:) = var_cllmh(cldfra3D(i,j,:), pres3D(i,j,:), d3)
+            END DO
+          END DO
+        ELSE
+          PRINT *,TRIM(ErrWarnMsg('err'))
+          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
+          PRINT *,'    accepted values: 1,2,3'
+          STOP
+        END IF
+      CASE (4)
+        IF (zdim == 1) THEN
+          DO i=1, d2
+            DO j=1, d3
+              DO k=1, d4
+                cllmh4D1(i,j,k,:) = var_cllmh(cldfra4D(:,i,j,k), pres4D(:,i,j,k), d1)
+              END DO
+            END DO
+          END DO
+        ELSE IF (zdim == 2) THEN
+          DO i=1, d1
+            DO j=1, d3
+              DO k=1, d4
+                cllmh4D2(i,j,k,:) = var_cllmh(cldfra4D(i,:,j,k), pres4D(i,:,j,k), d2)
+              END DO
+            END DO
+          END DO
+        ELSE IF (zdim == 3) THEN
+          DO i=1, d2
+            DO j=1, d3
+              DO k=1, d4
+                cllmh4D3(i,j,k,:) = var_cllmh(cldfra4D(i,j,:,k), pres4D(i,j,:,k), d3)
+              END DO
+            END DO
+          END DO
+        ELSE IF (zdim == 4) THEN
+          DO i=1, d1
+            DO j=1, d2
+              DO k=1, d3
+                cllmh4D4(i,j,k,:) = var_cllmh(cldfra4D(i,j,k,:), pres4D(i,j,k,:), d4)
+              END DO
+            END DO
+          END DO
+        ELSE
+          PRINT *,TRIM(ErrWarnMsg('err'))
+          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
+          PRINT *,'    accepted values: 1,2,3,4'
+          STOP
+        END IF
+      CASE DEFAULT
+        PRINT *,TRIM(ErrWarnMsg('err'))
+        PRINT *,'  ' // TRIM(fname) // ': Ndim:', Ndim,' not ready !!'
+        STOP
+      END SELECT
+
+    RETURN 
+
+  END SUBROUTINE compute_cllmh
+
+  SUBROUTINE compute_clt4D2(cldfra4D, clt4D2, d1, d2, d3, d4)
+! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ from a 4D CLDFRA
+!   where zdim is the 2nd dimension (thus, cldfra4D(d1,d2,d3,d4) --> clt(d1,d3,d4)
+! It should be properly done via an 'INTERFACE', but...
+
+    IMPLICIT NONE
+
+    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
+    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
+    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: cldfra4D
+    REAL(r_k), DIMENSION(d1,d3,d4), INTENT(out)          :: clt4D2
+
+! Local
+    INTEGER                                              :: i,j,k, zdim, Ndim
+    CHARACTER(LEN=50)                                    :: fname
+
+!!!!!!! Variables
+! cldfra4D: 4D cloud fraction values [1]
+! Ndim: number of dimensions of the input data
+! d[1-4]: dimensions of 'cldfra4D'
+! zdim: number of the vertical-dimension within the matrix
+! clt4D2: total cloudiness for the 4D cldfra and d2 being 'zdim'
+
+    fname = 'compute_clt4D2'
+    zdim = 2
+    Ndim = 4
+
+    DO i=1, d1
+      DO j=1, d3
+        DO k=1, d4
+          clt4D2(i,j,k) = var_clt(cldfra4D(i,:,j,k), d2)
+        END DO
+      END DO
+    END DO
+    
+    RETURN 
+
+  END SUBROUTINE compute_clt4D2
+
+  SUBROUTINE compute_clt3D1(cldfra3D, clt3D1, d1, d2, d3)
+! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ from a 3D CLDFRA
+!   where zdim is the 1st dimension (thus, cldfra4D(d1,d2,d3) --> clt(d2,d3)
+! It should be properly done via an 'INTERFACE', but...
+
+    IMPLICIT NONE
+
+    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
+    INTEGER, INTENT(in)                                  :: d1, d2, d3
+    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: cldfra3D
+    REAL(r_k), DIMENSION(d2,d3), INTENT(out)             :: clt3D1
+
+! Local
+    INTEGER                                              :: i,j,k, zdim, Ndim
+    CHARACTER(LEN=50)                                    :: fname
+
+!!!!!!! Variables
+! cldfra3D: 3D cloud fraction values [1]
+! Ndim: number of dimensions of the input data
+! d[1-3]: dimensions of 'cldfra3D'
+! zdim: number of the vertical-dimension within the matrix
+! clt3D1: total cloudiness for the 3D cldfra and d1 being 'zdim'
+
+    fname = 'compute_clt3D1'
+    zdim = 1
+    Ndim = 3
+
+    DO i=1, d2
+      DO j=1, d3
+        clt3D1(i,j) = var_clt(cldfra3D(:,i,j), d1)
+      END DO
+    END DO
+    
+    RETURN 
+
+  END SUBROUTINE compute_clt3D1
+
+  SUBROUTINE compute_clt(cldfra1D, cldfra2D, cldfra3D, cldfra4D, Ndim, zdim, clt1D, clt2D1, clt2D2,   &
+    clt3D1, clt3D2, clt3D3, clt4D1, clt4D2, clt4D3, clt4D4, d1, d2, d3, d4)
+! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ
+
+    IMPLICIT NONE
+
+    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
+    INTEGER, INTENT(in)                                  :: Ndim, d1, d2, d3, d4, zdim
+    REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in)       :: cldfra1D
+    REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(in)    :: cldfra2D
+    REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL, INTENT(in) :: cldfra3D
+    REAL(r_k), DIMENSION(d1,d2,d3,d4), OPTIONAL,                                                      &
+      INTENT(in)                                         :: cldfra4D
+    REAL(r_k), OPTIONAL, INTENT(out)                     :: clt1D
+    REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(out)      :: clt2D1
+    REAL(r_k), DIMENSION(d2), OPTIONAL, INTENT(out)      :: clt2D2
+    REAL(r_k), DIMENSION(d2,d3), OPTIONAL, INTENT(out)   :: clt3D1
+    REAL(r_k), DIMENSION(d1,d3), OPTIONAL, INTENT(out)   :: clt3D2
+    REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(out)   :: clt3D3
+    REAL(r_k), DIMENSION(d2,d3,d4), OPTIONAL,INTENT(out) :: clt4D1
+    REAL(r_k), DIMENSION(d1,d3,d4), OPTIONAL,INTENT(out) :: clt4D2
+    REAL(r_k), DIMENSION(d1,d2,d4), OPTIONAL,INTENT(out) :: clt4D3
+    REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL,INTENT(out) :: clt4D4
+
+! Local
+    INTEGER                                              :: i,j,k
+    CHARACTER(LEN=50)                                    :: fname
+
+!!!!!!! Variables
+! cldfra[1-4]D: cloud fraction values [1]
+! Ndim: number of dimensions of the input data
+! d[1-4]: dimensions of 'cldfra'
+! zdim: number of the vertical-dimension within the matrix
+! clt1D: total cloudiness for the 1D cldfra
+! clt2D1: total cloudiness for the 2D cldfra and d1 being 'zdim'
+! clt2D2: total cloudiness for the 2D cldfra and d2 being 'zdim'
+! clt3D1: total cloudiness for the 3D cldfra and d1 being 'zdim'
+! clt3D2: total cloudiness for the 3D cldfra and d2 being 'zdim'
+! clt3D3: total cloudiness for the 3D cldfra and d3 being 'zdim'
+! clt4D1: total cloudiness for the 4D cldfra and d1 being 'zdim'
+! clt4D2: total cloudiness for the 4D cldfra and d2 being 'zdim'
+! clt4D3: total cloudiness for the 4D cldfra and d3 being 'zdim'
+! clt4D4: total cloudiness for the 4D cldfra and d4 being 'zdim'
+
+    fname = 'compute_clt'
+
+    SELECT CASE (Ndim)
+      CASE (1)
+        clt1D = var_clt(cldfra1D, d1)
+      CASE (2)
+        IF (zdim == 1) THEN
+          DO i=1, d2
+            clt2D1(i) = var_clt(cldfra2D(:,i), d1)
+          END DO
+        ELSE IF (zdim == 2) THEN
+          DO i=1, d1
+            clt2D2(i) = var_clt(cldfra2D(:,i), d2)
+          END DO
+        ELSE
+          PRINT *,TRIM(ErrWarnMsg('err'))
+          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
+          PRINT *,'    accepted values: 1,2'
+          STOP
+        END IF
+      CASE (3)
+        IF (zdim == 1) THEN
+          DO i=1, d2
+            DO j=1, d3
+              clt3D1(i,j) = var_clt(cldfra3D(:,i,j), d1)
+            END DO
+          END DO
+        ELSE IF (zdim == 2) THEN
+          DO i=1, d1
+            DO j=1, d3
+              clt3D2(i,j) = var_clt(cldfra3D(i,:,j), d2)
+            END DO
+          END DO
+        ELSE IF (zdim == 3) THEN
+          DO i=1, d1
+            DO j=1, d2
+              clt3D3(i,j) = var_clt(cldfra3D(i,j,:), d3)
+            END DO
+          END DO
+        ELSE
+          PRINT *,TRIM(ErrWarnMsg('err'))
+          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
+          PRINT *,'    accepted values: 1,2,3'
+          STOP
+        END IF
+      CASE (4)
+        IF (zdim == 1) THEN
+          DO i=1, d2
+            DO j=1, d3
+              DO k=1, d4
+                clt4D1(i,j,k) = var_clt(cldfra4D(:,i,j,k), d1)
+              END DO
+            END DO
+          END DO
+        ELSE IF (zdim == 2) THEN
+          DO i=1, d1
+            DO j=1, d3
+              DO k=1, d4
+                clt4D2(i,j,k) = var_clt(cldfra4D(i,:,j,k), d2)
+              END DO
+            END DO
+          END DO
+        ELSE IF (zdim == 3) THEN
+          DO i=1, d2
+            DO j=1, d3
+              DO k=1, d4
+                clt4D3(i,j,k) = var_clt(cldfra4D(i,j,:,k), d3)
+              END DO
+            END DO
+          END DO
+        ELSE IF (zdim == 4) THEN
+          DO i=1, d1
+            DO j=1, d2
+              DO k=1, d3
+                clt4D4(i,j,k) = var_clt(cldfra4D(i,j,k,:), d4)
+              END DO
+            END DO
+          END DO
+        ELSE
+          PRINT *,TRIM(ErrWarnMsg('err'))
+          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
+          PRINT *,'    accepted values: 1,2,3,4'
+          STOP
+        END IF
+      CASE DEFAULT
+        PRINT *,TRIM(ErrWarnMsg('err'))
+        PRINT *,'  ' // TRIM(fname) // ': Ndim:', Ndim,' not ready !!'
+        STOP
+      END SELECT
+
+    RETURN 
+
+  END SUBROUTINE compute_clt
+
+END MODULE module_ForDiagnostics
Index: unk/tools/module_ForDiagnosticsVars.F90
===================================================================
--- /trunk/tools/module_ForDiagnosticsVars.F90	(revision 1605)
+++ 	(revision )
@@ -1,124 +1,0 @@
-!! Fortran version of different diagnostics
-! L. Fita. LMD May 2016
-! gfortran module_generic.o -c module_ForDiagnostics.F90
-!
-! f2py -m module_ForDiagnostics --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForDiagnostics.F90
-
-MODULE module_ForDiagnosticsVars
-
-  USE module_generic
-
-  IMPLICIT NONE
-
-!  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
-  REAL(r_k), PARAMETER                                   :: ZEPSEC=1.0D-12
-! Low limit pressure for medium clouds [Pa]
-  REAL(r_k), PARAMETER                                   :: prmhc = 44000.d0
-! Low limit pressure for High clouds [Pa]
-  REAL(r_k), PARAMETER                                   :: prmlc = 68000.d0
-
-  REAL(r_k), PARAMETER                                   :: zero=0.d0
-  REAL(r_k), PARAMETER                                   :: one=1.d0
-  REAL(r_k), PARAMETER                                   :: two=2.d0
-
-  CONTAINS
-
-!!!!!!! Variables
-! var_cllmh: low, medium, high-cloud [0,1]
-! var_clt: total cloudiness [0,1]
-
-!!!!!!! Calculations
-! compute_clt: Computation of total cloudiness
-
-!!!
-! Variables
-!!!
-
-  FUNCTION var_cllmh(clfra, p, dz)
-! Function to compute cllmh on a 1D column 1: low-cloud; 2: medium-cloud; 3: high-cloud [1]
-
-    IMPLICIT NONE
-
-!    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
-    INTEGER, INTENT(in)                                  :: dz
-    REAL(r_k), DIMENSION(dz), INTENT(in)                 :: clfra, p
-    REAL(r_k), DIMENSION(3)                              :: var_cllmh
-
-! Local
-    INTEGER                                              :: iz
-    CHARACTER(LEN=50)                                    :: fname
-    REAL(r_k)                                            :: zclearl, zcloudl, zclearm, zcloudm,       &
-      zclearh, zcloudh
-
-!!!!!!! Variables
-! clfra: cloudfraction as 1D verical-column [1]
-! p: pressure values of the column
-    fname = 'var_cllmh'
-
-    zclearl = one
-    zcloudl = zero
-    zclearm = one
-    zcloudm = zero
-    zclearh = one
-    zcloudh = zero
-
-    var_cllmh = one
-
-    DO iz=1, dz
-      IF (p(iz) < prmhc) THEN
-        var_cllmh(3) = var_cllmh(3)*(one-MAX(clfra(iz),zcloudh))/(one-MIN(zcloudh,one-ZEPSEC))
-        zcloudh = clfra(iz)
-      ELSE IF ( (p(iz) >= prmhc) .AND. (p(iz) < prmlc) ) THEN
-        var_cllmh(2) = var_cllmh(2)*(one-MAX(clfra(iz),zcloudm))/(one-MIN(zcloudm,one-ZEPSEC))
-        zcloudm = clfra(iz)
-      ELSE IF (p(iz) >= prmlc) THEN
-        var_cllmh(1) = var_cllmh(1)*(one-MAX(clfra(iz),zcloudl))/(one-MIN(zcloudl,one-ZEPSEC))
-        zcloudl = clfra(iz)
-      ELSE
-        PRINT *,'  ' // TRIM(fname) // ': This is weird, pressure:', p(iz), ' Pa fails out!!'
-        PRINT *,'    from high, low cloud pressures:', prmhc, ' ,', prmlc,' Pa at z-level:', iz
-        PRINT *,'    p_high > p:', prmhc,'> ',p(iz),' Pa'
-        PRINT *,'    p_low > p >= p_high:', prmlc,'> ',p(iz),' >=', prmhc,' Pa'
-        PRINT *,'    p_low >= p:', prmlc,'>= ',p(iz),' Pa'
-        STOP
-      END IF
-    END DO
-
-    var_cllmh = one - var_cllmh
-
-    RETURN 
-
-  END FUNCTION var_cllmh
-
-  REAL(r_k) FUNCTION var_clt(clfra, dz)
-! Function to compute the total cloud following 'newmicro.F90' from LMDZ using 1D vertical 
-!   column values
-
-    IMPLICIT NONE
-
-!    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
-    INTEGER, INTENT(in)                                  :: dz
-    REAL(r_k), DIMENSION(dz), INTENT(in)                 :: clfra
-! Local
-    INTEGER                                              :: iz
-    REAL(r_k)                                            :: zclear, zcloud
-    CHARACTER(LEN=50)                                    :: fname
-!!!!!!! Variables
-! cfra: 1-column cloud fraction values
-
-    fname = 'var_clt'
-
-    zclear = one
-    zcloud = zero
-
-    DO iz=1,dz
-      zclear = zclear*(one-MAX(clfra(iz),zcloud))/(one-MIN(zcloud,1.-ZEPSEC))
-      var_clt = one - zclear
-      zcloud = clfra(iz)
-    END DO
-
-    RETURN
-
-  END FUNCTION var_clt
-
-END MODULE module_ForDiagnosticsVars
Index: /trunk/tools/module_ForDiagnosticsVars.f90
===================================================================
--- /trunk/tools/module_ForDiagnosticsVars.f90	(revision 1606)
+++ /trunk/tools/module_ForDiagnosticsVars.f90	(revision 1606)
@@ -0,0 +1,124 @@
+!! Fortran version of different diagnostics
+! L. Fita. LMD May 2016
+! gfortran module_generic.o -c module_ForDiagnostics.F90
+!
+! f2py -m module_ForDiagnostics --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForDiagnostics.F90
+
+MODULE module_ForDiagnosticsVars
+
+  USE module_generic
+
+  IMPLICIT NONE
+
+!  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
+  REAL(r_k), PARAMETER                                   :: ZEPSEC=1.0D-12
+! Low limit pressure for medium clouds [Pa]
+  REAL(r_k), PARAMETER                                   :: prmhc = 44000.d0
+! Low limit pressure for High clouds [Pa]
+  REAL(r_k), PARAMETER                                   :: prmlc = 68000.d0
+
+  REAL(r_k), PARAMETER                                   :: zero=0.d0
+  REAL(r_k), PARAMETER                                   :: one=1.d0
+  REAL(r_k), PARAMETER                                   :: two=2.d0
+
+  CONTAINS
+
+!!!!!!! Variables
+! var_cllmh: low, medium, high-cloud [0,1]
+! var_clt: total cloudiness [0,1]
+
+!!!!!!! Calculations
+! compute_clt: Computation of total cloudiness
+
+!!!
+! Variables
+!!!
+
+  FUNCTION var_cllmh(clfra, p, dz)
+! Function to compute cllmh on a 1D column 1: low-cloud; 2: medium-cloud; 3: high-cloud [1]
+
+    IMPLICIT NONE
+
+!    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
+    INTEGER, INTENT(in)                                  :: dz
+    REAL(r_k), DIMENSION(dz), INTENT(in)                 :: clfra, p
+    REAL(r_k), DIMENSION(3)                              :: var_cllmh
+
+! Local
+    INTEGER                                              :: iz
+    CHARACTER(LEN=50)                                    :: fname
+    REAL(r_k)                                            :: zclearl, zcloudl, zclearm, zcloudm,       &
+      zclearh, zcloudh
+
+!!!!!!! Variables
+! clfra: cloudfraction as 1D verical-column [1]
+! p: pressure values of the column
+    fname = 'var_cllmh'
+
+    zclearl = one
+    zcloudl = zero
+    zclearm = one
+    zcloudm = zero
+    zclearh = one
+    zcloudh = zero
+
+    var_cllmh = one
+
+    DO iz=1, dz
+      IF (p(iz) < prmhc) THEN
+        var_cllmh(3) = var_cllmh(3)*(one-MAX(clfra(iz),zcloudh))/(one-MIN(zcloudh,one-ZEPSEC))
+        zcloudh = clfra(iz)
+      ELSE IF ( (p(iz) >= prmhc) .AND. (p(iz) < prmlc) ) THEN
+        var_cllmh(2) = var_cllmh(2)*(one-MAX(clfra(iz),zcloudm))/(one-MIN(zcloudm,one-ZEPSEC))
+        zcloudm = clfra(iz)
+      ELSE IF (p(iz) >= prmlc) THEN
+        var_cllmh(1) = var_cllmh(1)*(one-MAX(clfra(iz),zcloudl))/(one-MIN(zcloudl,one-ZEPSEC))
+        zcloudl = clfra(iz)
+      ELSE
+        PRINT *,'  ' // TRIM(fname) // ': This is weird, pressure:', p(iz), ' Pa fails out!!'
+        PRINT *,'    from high, low cloud pressures:', prmhc, ' ,', prmlc,' Pa at z-level:', iz
+        PRINT *,'    p_high > p:', prmhc,'> ',p(iz),' Pa'
+        PRINT *,'    p_low > p >= p_high:', prmlc,'> ',p(iz),' >=', prmhc,' Pa'
+        PRINT *,'    p_low >= p:', prmlc,'>= ',p(iz),' Pa'
+        STOP
+      END IF
+    END DO
+
+    var_cllmh = one - var_cllmh
+
+    RETURN 
+
+  END FUNCTION var_cllmh
+
+  REAL(r_k) FUNCTION var_clt(clfra, dz)
+! Function to compute the total cloud following 'newmicro.F90' from LMDZ using 1D vertical 
+!   column values
+
+    IMPLICIT NONE
+
+!    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
+    INTEGER, INTENT(in)                                  :: dz
+    REAL(r_k), DIMENSION(dz), INTENT(in)                 :: clfra
+! Local
+    INTEGER                                              :: iz
+    REAL(r_k)                                            :: zclear, zcloud
+    CHARACTER(LEN=50)                                    :: fname
+!!!!!!! Variables
+! cfra: 1-column cloud fraction values
+
+    fname = 'var_clt'
+
+    zclear = one
+    zcloud = zero
+
+    DO iz=1,dz
+      zclear = zclear*(one-MAX(clfra(iz),zcloud))/(one-MIN(zcloud,1.-ZEPSEC))
+      var_clt = one - zclear
+      zcloud = clfra(iz)
+    END DO
+
+    RETURN
+
+  END FUNCTION var_clt
+
+END MODULE module_ForDiagnosticsVars
Index: unk/tools/module_ForInterpolate.F90
===================================================================
--- /trunk/tools/module_ForInterpolate.F90	(revision 1605)
+++ 	(revision )
@@ -1,1671 +1,0 @@
-! Module to interpolate values from a giving projection and pressure interpolation
-! To be included in a python
-! 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
-
-  USE module_generic
-
-  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
-
-  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 CoarselonlatFindExact(dx, dy, ilon, ilat, nxlon, nxlat, fracx, fracy, fraclon, fraclat,    &
-  iv, lonv, latv, per, Nperx, Npery, mindiff, ilonlat, mindiffLl)
-! Function to search a given value from a coarser version of the data
-
-  IMPLICIT NONE
-
-!  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
-  INTEGER, INTENT(in)                                    :: dx, dy, iv
-  REAL(r_k), DIMENSION(dx,dy), INTENT(in)                :: ilon, ilat
-  INTEGER, INTENT(in)                                    :: fracx, fracy
-  REAL(r_k), DIMENSION(Nperx,Npery), INTENT(in)          :: fraclon, fraclat
-  REAL(r_k), INTENT(in)                                  :: lonv, latv, per, mindiff
-  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
-  INTEGER                                                :: i
-  REAL(r_k), DIMENSION(Nperx,Npery)                      :: difffraclonlat
-  REAL(r_k)                                              :: mindifffracLl
-  INTEGER, DIMENSION(2)                                  :: ilonlatfrac
-  INTEGER                                                :: ixbeg, ixend, iybeg, iyend
-  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
-! iv: point in the input data
-! 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
-! frac[x/y]: Number of grid points for each fraction
-! fraclon, fraclat: longitude and latitude fractional matricies to perform the first guess
-! mindiff: authorized minimal distance between input and interpolated point
-! ilonlat: grid point on the total lon,lat matrix
-! mindiffLl: distance between input and interpolated point
-
-  fname = 'CoarselonlatFindExact'
-
-  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
-
-! 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
-    PRINT *,'Lluis!',fraclonv, '>=', lonv,'&', fraclatv, '>=', latv
-    IF (ilonlatfrac(1) > 1) THEN
-      ixbeg = (ilonlatfrac(1)-1)*fracx
-      ixend = ilonlatfrac(1)*fracx+1
-    ELSE
-      PRINT *,'Lluis 2'
-      ixbeg = 1
-      ixend = fracx+1
-    END IF
-    IF (ilonlatfrac(2) > 1) THEN
-      iybeg = (ilonlatfrac(2)-1)*fracy
-      iyend = ilonlatfrac(2)*fracy+1
-    ELSE
-      iybeg = 1
-      iyend = fracy+1
-    END IF
-  ELSE IF (fraclonv < lonv .AND. fraclatv >= latv) THEN
-    PRINT *,'Lluis!',fraclonv, '<', lonv,'&', fraclatv, '>=', latv
-    IF (ilonlatfrac(1) < Nperx) THEN
-      PRINT *,'Lluis 2'
-      IF (ilonlatfrac(1) /= 1) THEN
-        ixbeg = (ilonlatfrac(1)-1)*fracx
-        ixend = ilonlatfrac(1)*fracx+1
-      ELSE
-        ixbeg = 1
-        ixend = fracx+1
-      END IF
-    ELSE
-      ixbeg = Nperx*fracx
-      ixend = dx+1
-    END IF
-    IF (ilonlatfrac(2) > 1) THEN
-      iybeg = (ilonlatfrac(2)-1)*fracy
-      iyend = ilonlatfrac(2)*fracy+1
-    ELSE
-      iybeg = 1
-      iyend = fracy+1
-    END IF    
-  ELSE IF (fraclonv < lonv .AND. fraclatv < latv) THEN
-    PRINT *,'Lluis!',fraclonv, '<', lonv,'&', fraclatv, '<', latv
-    IF (ilonlatfrac(1) < Nperx) THEN
-      IF (ilonlatfrac(1) /= 1) THEN
-        ixbeg = (ilonlatfrac(1)-1)*fracx
-        ixend = ilonlatfrac(1)*fracx+1
-      ELSE
-        ixbeg = 1
-        ixend = fracx+1
-      END IF
-    ELSE
-      ixbeg = Nperx*fracx
-      ixend = dx+1
-    ENDIF
-    IF (ilonlatfrac(2) < Npery) THEN
-      IF (ilonlatfrac(2) /= 1) THEN
-        iybeg = (ilonlatfrac(2)-1)*fracy
-        iyend = ilonlatfrac(2)*fracy+1
-      ELSE
-        iybeg = 1
-        iyend = fracy+1
-      END IF
-    ELSE
-      iybeg = Npery*fracy
-      iyend = dy+1
-    END IF
-  ELSE IF (fraclonv >= lonv .AND. fraclatv < latv) THEN
-    PRINT *,'Llui!',fraclonv, '>=', lonv,'&', fraclatv, '<', latv
-    IF (ilonlatfrac(1) > 1) THEN
-      ixbeg = (ilonlatfrac(1)-1)*fracx
-      ixend = ilonlatfrac(1)*fracx+1
-    ELSE
-      ixbeg = 1
-      ixend = fracx+1
-    END IF
-    IF (ilonlatfrac(2) < Npery) THEN
-      IF (ilonlatfrac(2) /= 1) THEN
-        iybeg = (ilonlatfrac(2)-1)*fracy
-        iyend = ilonlatfrac(2)*fracy+1
-      ELSE
-        iybeg = 1
-        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)
-
-  IF (mindiffLl > mindiff) THEN 
-    difflonlat = SQRT((lon-lonv)**2. + (lat-latv)**2.)
-    mindiffLl = MINVAL(difflonlat)
-  END IF
-
-  IF (mindiffLl > mindiff) THEN
-    PRINT *,TRIM(ErrWarnMsg('err'))
-    PRINT *,'  ' // TRIM(fname) // ': not equivalent point closer than:',mindiff,' found!!'
-    PRINT *,'    at input point iv:', iv,' lon/lat:', lonv,', ',latv,' distance:',mindiffLl
-    PRINT *,'    Fraction values _______ (',Nperx,', ',Npery ,')'
-    PRINT *,'      fraclon'
-    DO i=1, Nperx
-      PRINT *,'      ',fraclon(i,:)
-    END DO
-    PRINT *,'      fraclat'
-    DO i=1, Nperx
-      PRINT *,'      ',fraclat(i,:)
-    END DO
-    PRINT *,'      frac lon, lat:', fraclon(ilonlatfrac(1),ilonlatfrac(2)), ' ,',                     &
-     fraclat(ilonlatfrac(1),ilonlatfrac(2))
-    PRINT *,'      mindifffracLl:', mindifffracLl, ' ilonlatfrac:', ilonlatfrac
-    PRINT *,'    Coarse values _______'
-    PRINT *,'      indices. x:', ixbeg, ', ', ixend, ' y:', iybeg, ', ', iyend
-    PRINT *,'      lon range:', '(',ilon(ixbeg,iybeg),', ',ilon(ixend,iyend),')'
-    PRINT *,'      lat range:', '(',ilat(ixbeg,iybeg),', ',ilat(ixend,iyend),')'
-    PRINT *,'      lon', UBOUND(lon)
-    DO i=1, ixend-ixbeg+1
-      PRINT *,'      ',lon(i,:)
-    END DO
-    PRINT *,'      lat', UBOUND(lat)
-    DO i=1, ixend-ixbeg+1
-      PRINT *,'      ',lat(i,:)
-    END DO
-    STOP
-  END IF
-
-  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 CoarselonlatFindExact
-
-SUBROUTINE lonlatFind(dx, dy, ilon, ilat, nxlon, nxlat, lonv, latv, ilonlat, mindiffLl)
-! Function to search a given value from a coarser version of the data
-
-  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), INTENT(in)                                  :: lonv, latv
-  REAL(r_k), DIMENSION(2), INTENT(in)                    :: nxlon, nxlat
-  INTEGER, DIMENSION(2), INTENT(out)                     :: ilonlat
-  REAL(r_k), INTENT(out)                                 :: mindiffLl
-! Local
-  REAL(r_k), DIMENSION(dx,dy)                            :: difflonlat
-  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
-
-  fname = 'lonlatFind'
-
-  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
-
-! Find point
-  difflonlat = SQRT((ilon-lonv)**2. + (ilat-latv)**2.)
-  mindiffLl = MINVAL(difflonlat)
-  ilonlat = index2DArrayR(difflonlat, dx, dy, mindiffLl)
-
-!  PRINT *,'mindiffLl:', mindiffLl, ' ilatlon:', ilatlon
-!  PRINT *,'lon, lat:', lon(ilonlat(1),ilonlat(2)), ' ,', lat(ilonlat(1),ilonlat(2))
-
-  RETURN
-
-END SUBROUTINE lonlatFind
-
-SUBROUTINE CoarseInterpolate(projlon, projlat, lonvs, latvs, percen, mindiff, inpt, ilonlat,          &
-  mindiffLl, dimx, dimy, Ninpts)
-! Subroutine which finds the closest grid point within a projection throughout a first guest 
-!   approche from percentages of the whole domain
-
-  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)               :: inpt, lonvs, latvs
-  REAL(r_k), INTENT(in)                                  :: mindiff, percen
-  INTEGER, DIMENSION(Ninpts,2), INTENT(out)              :: ilonlat
-  REAL(r_k), DIMENSION(Ninpts), INTENT(out)              :: mindiffLl
-
-! Local
-  INTEGER                                                :: iv,i,j
-  INTEGER                                                :: ierr
-  INTEGER                                                :: Ninpts1
-  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
-  REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: fractionlon, fractionlat
-  INTEGER                                                :: dfracdx, dfracdy, 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
-! inpt: whether the point has already been localized (1) or not (0)
-! ilonlat: Longitude and Latitude of the input points
-! mindiffLl: minimum difference between target and source longitude/latitude (in degrees)
-
-  fname = 'CoarseInterpolate'
-  Ninpts1 = Ninpts/100
-
-  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
-  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
-
-  PRINT *,'  ' // TRIM(fname) //' total space:', dimx, ', ', dimy, ' %', percen
-
-  dfracdx = INT(1./percen+1)
-  dfracdy = INT(1./percen+1)
-  fracdx = INT(dimx*percen)
-  fracdy = INT(dimy*percen)
-  PRINT *,'  ' // TRIM(fname) //' fraction:', dfracdx, ', ', dfracdy, ' freq:', fracdx,', ',fracdy
-
-  IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon)
-  ALLOCATE(fractionlon(dfracdx, dfracdy), STAT=ierr)
-  IF (ierr /= 0) THEN
-    PRINT *,TRIM(ErrWarnMsg('err'))
-    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlon' !!"
-    STOP
-  END IF
-  IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat)
-  ALLOCATE(fractionlat(dfracdx, dfracdy), STAT=ierr)
-  IF (ierr /= 0) THEN
-    PRINT *,TRIM(ErrWarnMsg('err'))
-    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlat' !!"
-    STOP
-  END IF
-
-  DO i=1,dfracdx
-    DO j=1,dfracdy
-      fractionlon(i,j) = projlon(fracdx*(i-1)+1,fracdy*(j-1)+1)
-      fractionlat(i,j) = projlat(fracdx*(i-1)+1,fracdy*(j-1)+1)
-!      PRINT *,'i,j:',i,', ',j,' frac ij:',fracdx*(i-1),', ',fracdy*(j-1),' lonlat:', fractionlon(i,j),&
-!        ', ',fractionlat(i,j)
-    END DO
-  END DO
-
-!  PRINT *,'  ' // TRIM(fname) // ' fractions of:'
-!  PRINT *,' lon _______ (',dfracdx,', ',dfracdy,')'
-!  DO i=1,dfracdx
-!    PRINT *,fractionlon(i,:)
-!  END DO
-!  PRINT *,' lat_______'
-!  DO i=1,dfracdx
-!    PRINT *,fractionlat(i,:)
-!  END DO
-
-  DO iv=1,Ninpts
-    IF (inpt(iv) == 0) THEN
-      CALL CoarselonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, fractionlon,        &
-        fractionlat, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, ilonlat(iv,:), mindiffLl(iv))
-
-      IF ((mindiffLl(iv) <= mindiff) .AND. .NOT.(ilonlat(iv,1) >= 0 .AND. ilonlat(iv,1) >= 0)) THEN
-        PRINT *,TRIM(ErrWarnMsg('err'))
-        PRINT *,'  ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv),        &
-          ' not relocated !!'
-        PRINT *,'    mindiffl:', mindiffLl(iv), ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2)
-        STOP
-      END IF
-
-    END IF
-  END DO
-
-END SUBROUTINE CoarseInterpolate
-
-SUBROUTINE CoarseInterpolateExact(projlon, projlat, lonvs, latvs, percen, mindiff, ivar, newvar,      &
-  newvarin, newvarinpt, newvarindiff, dimx, dimy, Ninpts)
-! Subroutine which finds the closest grid point within a projection throughout a first guest 
-!   and then whole domain approche from percentages of the whole domain 
-
-  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), INTENT(in)                                  :: 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,i,j
-  INTEGER                                                :: ierr
-  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                                                :: dfracdx, dfracdy, 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 = 'CoarseInterpolateExact'
-  Ninpts1 = Ninpts/100
-
-  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
-  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
-
-  PRINT *,'  ' // TRIM(fname) //' total space:', dimx, ', ', dimy, ' %', percen
-
-  dfracdx = INT(1./percen+1)
-  dfracdy = INT(1./percen+1)
-  fracdx = INT(dimx*percen)
-  fracdy = INT(dimy*percen)
-  PRINT *,'  ' // TRIM(fname) //' fraction:', dfracdx, ', ', dfracdy, ' freq:', fracdx,', ',fracdy
-
-  IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon)
-  ALLOCATE(fractionlon(dfracdx, dfracdy), STAT=ierr)
-  IF (ierr /= 0) THEN
-    PRINT *,TRIM(ErrWarnMsg('err'))
-    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlon' !!"
-    STOP
-  END IF
-  IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat)
-  ALLOCATE(fractionlat(dfracdx, dfracdy), STAT=ierr)
-  IF (ierr /= 0) THEN
-    PRINT *,TRIM(ErrWarnMsg('err'))
-    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlat' !!"
-    STOP
-  END IF
-
-  DO i=1,dfracdx
-    DO j=1,dfracdy
-      fractionlon(i,j) = projlon(fracdx*(i-1)+1,fracdy*(j-1)+1)
-      fractionlat(i,j) = projlat(fracdx*(i-1)+1,fracdy*(j-1)+1)
-!      PRINT *,'i,j:',i,', ',j,' frac ij:',fracdx*(i-1),', ',fracdy*(j-1),' lonlat:', fractionlon(i,j),&
-!        ', ',fractionlat(i,j)
-    END DO
-  END DO
-
-!  PRINT *,'  ' // TRIM(fname) // ' fractions of:'
-!  PRINT *,' lon _______ (',dfracdx,', ',dfracdy,')'
-!  DO i=1,dfracdx
-!    PRINT *,fractionlon(i,:)
-!  END DO
-!  PRINT *,' lat_______'
-!  DO i=1,dfracdx
-!    PRINT *,fractionlat(i,:)
-!  END DO
-
-  DO iv=1,Ninpts
-    IF (newvarinpt(iv) == 0) THEN
-      CALL CoarselonlatFindExact(dimx, dimy, projlon, projlat, extremelon, extremelat, fracdx, fracdy,&
-        fractionlon, fractionlat, iv, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, mindiff,        &
-        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 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
-
-  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)                                              :: inclon, inclat, maxdiffprojlonlat,        &
-    maxdiffinlonlat
-  REAL(r_k), DIMENSION(idimx,idimy)                      :: difflonlat
-  REAL(r_k), DIMENSION(idimx,idimy)                      :: idifflon, idifflat
-  REAL(r_k), DIMENSION(pdimx,pdimy)                      :: difflon, difflat
-  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat, ipos
-  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) /)
- 
-  ! Maximum distance between grid points in input projection
-  idifflon = 0.
-  idifflat = 0.
-  idifflon(1:idimx-1,:) = inlonv(2:idimx,:)-inlonv(1:idimx-1,:)
-  idifflat(:,1:idimy-1) = inlatv(:,2:idimy)-inlatv(:,1:idimy-1)
-  maxdiffinlonlat = MAXVAL(SQRT(idifflon**2. + idifflat**2.))
-  ! Maximum distance between grid points in target projection
-  difflon = 0.
-  difflat = 0.
-  difflon(1:pdimx-1,:) = projlon(2:pdimx,:)-projlon(1:pdimx-1,:)
-  difflat(:,1:pdimy-1) = projlat(:,2:pdimy)-projlat(:,1:pdimy-1)
-  maxdiffprojlonlat = MAXVAL(SQRT(difflon**2. + difflat**2.))
-
-  IF (maxdiffinlonlat > maxdiffprojlonlat) THEN
-    PRINT *,TRIM(warnmsg)
-    PRINT *,'  ' //TRIM(fname)// '; input resolution: ', maxdiffinlonlat, ' is coarser than target:', &
-      maxdiffprojlonlat, ' !!'
-  END IF
-
-  ! Using case outside loop to be more efficient
-  SELECT CASE(TRIM(intkind))
-    CASE ('dis')
-      inclon = inlonv(2,1) - inlonv(1,1)
-      inclat = inlatv(1,2) - inlatv(1,1)
-
-      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)
-          IF ( (mindiffLl > maxdiffprojlonlat) .AND. (mindiffLl > maxdiffinlonlat)) THEN
-            outLlw(3,:,i,j) = 0.
-            outLlw(3,:,i,j) = -1.
-          ELSE
-            ! Getting the four surrounding values
-            iLl = index2DArrayR(difflonlat, idimx, idimy, mindiffLl)
-            IF ( (projlon(i,j) < inlonv(iLl(1),iLl(2)) .AND. inclon > 0.) .OR.                          &
-              (projlon(i,j) > inlonv(iLl(1),iLl(2)) .AND. inclon < 0.) ) 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) = MAX(iLl(1),1)
-              outLlw(1,2,i,j) = MAX(iLl(1),1)
-              outLlw(1,3,i,j) = MIN(iLl(1)+1,idimx)
-              outLlw(1,4,i,j) = MIN(iLl(1)+1,idimx)
-            END IF
-            IF ( (projlat(i,j) < inlatv(iLl(2),iLl(2)) .AND. inclat > 0.) .OR.                          &
-              (projlat(i,j) > inlatv(iLl(2),iLl(2)) .AND. inclat < 0.) )  THEN
-              outLlw(2,1,i,j) = MAX(iLl(2)-1,1)
-              outLlw(2,2,i,j) = MIN(iLl(2),idimy)
-              outLlw(2,3,i,j) = MIN(iLl(2),idimy)
-              outLlw(2,4,i,j) = MAX(iLl(2)-1,1)
-            ELSE
-              outLlw(2,1,i,j) = MAX(iLl(2),1)
-              outLlw(2,2,i,j) = MIN(iLl(2)+1,idimy)
-              outLlw(2,3,i,j) = MIN(iLl(2)+1,idimy)
-              outLlw(2,4,i,j) = MAX(iLl(2),1)
-            END IF
-            ! Computing distances
-!Keep the print for future checks?
-!            IF (MOD(i+j,200) == 0) THEN
-!              PRINT *,'center point:',i,j,'=', iLl,':',projlon(i,j),projlat(i,j),' incs',inclon,' ,',inclat
-!              PRINT *,'SW:', outLlw(1,1,i,j), outLlw(2,1,i,j),':',inlonv(outLlw(1,1,i,j), outLlw(2,1,i,j)),&
-!                inlatv(outLlw(1,1,i,j), outLlw(2,1,i,j))
-!              PRINT *,'NW:', outLlw(1,2,i,j), outLlw(2,2,i,j),':',inlonv(outLlw(1,2,i,j), outLlw(2,2,i,j)),&
-!                inlatv(outLlw(1,2,i,j), outLlw(2,2,i,j))
-!              PRINT *,'NE:', outLlw(1,3,i,j), outLlw(2,3,i,j),':',inlonv(outLlw(1,3,i,j), outLlw(2,3,i,j)),&
-!                inlatv(outLlw(1,3,i,j), outLlw(2,3,i,j))
-!              PRINT *,'SE:', outLlw(1,4,i,j), outLlw(2,4,i,j),':',inlonv(outLlw(1,4,i,j), outLlw(2,4,i,j)),&
-!                inlatv(outLlw(1,4,i,j), outLlw(2,4,i,j))
-!            END IF
-            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
-!              IF (i+j == 2) PRINT *,'iv:',iv,'dist:',dist,'weight:',outLlw(3,iv,i,j)
-            END DO
-            ! Normalizing
-            outLlw(3,:,i,j) = outLlw(3,:,i,j)/(SUM(outLlw(3,:,i,j)))
-!             IF (i+j == 2) PRINT *,'Normalized weights:',outLlw(3,:,i,j),':',SUM(outLlw(3,:,i,j))
-          END IF
-        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)
-          ipos = index2DArrayR(difflonlat, idimx, idimy, mindiffLl)
-          outLlw(1,1,i,j) = REAL(ipos(1))
-          outLlw(2,1,i,j) = REAL(ipos(2))
-          ! We do not want that values larger that the maximum distance between target grid points
-!          PRINT *,i,j,':',mindiffLl,'maxdiffLl:',maxdiffprojlonlat
-          IF ((mindiffLl .gt. maxdiffprojlonlat) .AND. (mindiffLl > maxdiffinlonlat)) THEN
-!            PRINT *,'  ' // TRIM(fname) // ': reprojected minimum distance to nearest grid point:',   &
-!              mindiffLl, ' larger than the maximum distance between grid points in target projection!!'
-            outLlw(3,1,i,j) = -1.
-          ELSE
-            outLlw(3,1,i,j) = 1.
-          END IF
-          ix = INT(outLlw(1,1,i,j))
-          iy = INT(outLlw(2,1,i,j))
-        END DO
-      END DO
-  END SELECT
-
-END SUBROUTINE LlInterpolateProjection
-
-SUBROUTINE var2D_IntProj(var2Din, inlonv, inlatv, projlon, projlat, intkind, mask, varout, idimx,     &
-  idimy, pdimx, pdimy)
-! Subroutine to interpolate a 2D variable
-
-  IMPLICIT NONE
-
-  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(idimx,idimy), INTENT(in)          :: var2Din
-  INTEGER, DIMENSION(idimx,idimy), INTENT(in)            :: mask
-  REAL(r_k), DIMENSION(pdimx,pdimy), 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.)
-! var2Din: 2D variable to interpolate
-! mask: mask of the intpu values (1: good, 0: none)
-! 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, pdimy
-          IF (outLlw(3,1,i,j) == -1.) THEN
-            varout(i,j) = fillVal64
-          ELSE 
-            varout(i,j) = 0.
-            DO iv=1, 4
-              ix = INT(outLlw(1,iv,i,j))
-              iy = INT(outLlw(2,iv,i,j))
-              IF (mask(ix,iy) == 1) THEN
-                w =  outLlw(3,iv,i,j)
-                varout(i,j) = varout(i,j) + w*var2Din(ix,iy)
-              END IF
-            END DO
-          END IF
-        END DO
-      END DO
-    CASE('npp')
-      DO i=1, pdimx
-        DO j=1, pdimy
-          ix = INT(outLlw(1,1,i,j))
-          iy = INT(outLlw(2,1,i,j))
-          IF ( (outLlw(3,1,i,j) == -1.) .OR. (mask(ix,iy) == 0) ) THEN
-            varout(i,j) = fillVal64
-          ELSE 
-            varout(i,j) = var2Din(ix,iy)*outLlw(3,1,i,j)
-          END IF
-        END DO
-      END DO
-  END SELECT
-
-END SUBROUTINE var2D_IntProj
-
-
-SUBROUTINE var3D_IntProj(var3Din, inlonv, inlatv, projlon, projlat, intkind, mask, varout, idimx,     &
-  idimy, pdimx, pdimy, d3)
-! Subroutine to interpolate a 3D variable
-
-  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(idimx,idimy,d3), INTENT(in)       :: var3Din
-  INTEGER, DIMENSION(idimx,idimy,d3), INTENT(in)         :: mask
-  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
-! mask: mask of the intpu values (1: good, 0: none)
-! 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, pdimy
-          IF (ALL(outLlw(3,:,i,j) == -1.)) THEN
-            varout(i,j,:) = fillVal64
-          ELSE 
-            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))
-                IF (mask(ix,iy,k) == 1) THEN
-                  w =  outLlw(3,iv,i,j)
-                  varout(i,j,k) = varout(i,j,k) + w*var3Din(ix,iy,k)
-                END IF
-              END DO
-            END DO
-          END IF
-        END DO
-      END DO
-    CASE('npp')
-      DO i=1, pdimx
-        DO j=1, pdimy
-          ix = INT(outLlw(1,1,i,j))
-          iy = INT(outLlw(2,1,i,j))
-          IF ( (outLlw(3,1,i,j) == -1.) .OR. (mask(ix,iy,1) == 0) ) THEN
-            varout(i,j,:) = fillVal64
-          ELSE 
-            DO k=1, d3
-              varout(i,j,k) = var3Din(ix,iy,k)*outLlw(3,1,i,j)
-            END DO
-          END IF
-        END DO
-      END DO
-  END SELECT
-
-END SUBROUTINE var3D_IntProj
-
-SUBROUTINE var4D_IntProj(var4Din, inlonv, inlatv, projlon, projlat, intkind, mask, varout, idimx,     &
-  idimy, pdimx, pdimy, d3, d4)
-! Subroutine to interpolate a 4D variable
-
-  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(idimx,idimy,d3,d4), INTENT(in)    :: var4Din
-  INTEGER, DIMENSION(idimx,idimy,d3,d4), INTENT(in)      :: mask
-  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
-! mask: mask of the intpu values (1: good, 0: none)
-! 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, pdimy
-          IF (ALL(outLlw(3,:,i,j) == -1.)) THEN
-            varout(i,j,:,:) = fillVal64
-          ELSE 
-            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))
-                  IF (mask(ix,iy,k,l) == 1) THEN
-                    w =  outLlw(3,iv,i,j)
-                    varout(i,j,k,l) = varout(i,j,k,l) + w*var4Din(ix,iy,k,l)
-                  END IF
-                END DO
-              END DO
-            END DO
-          END IF
-        END DO
-      END DO
-    CASE('npp')
-      DO i=1, pdimx
-        DO j=1, pdimy
-          ix = INT(outLlw(1,1,i,j))
-          iy = INT(outLlw(2,1,i,j))
-          IF ( (outLlw(3,1,i,j) == -1.) .OR. (mask(ix,iy,1,1) == 0) ) THEN
-            varout(i,j,:,:) = fillVal64
-          ELSE 
-            DO k=1, d3
-              DO l=1, d4
-                varout(i,j,k,l) = var4Din(ix,iy,k,l)*outLlw(3,1,i,j)
-              END DO
-            END DO
-          END IF
-        END DO
-      END DO
-  END SELECT
-
-END SUBROUTINE var4D_IntProj
-
-SUBROUTINE var5D_IntProj(var5Din, inlonv, inlatv, projlon, projlat, intkind, mask, varout, idimx,     &
-  idimy, pdimx, pdimy, d3, d4, d5)
-! Subroutine to interpolate a 5D variable
-
-  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(idimx,idimy,d3,d4,d5), INTENT(in) :: var5Din
-  INTEGER, DIMENSION(idimx,idimy,d3,d4,d5), INTENT(in)   :: mask
-  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
-! mask: mask of the intpu values (1: good, 0: none)
-! 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, pdimy
-          IF (ALL(outLlw(3,:,i,j) == -1.)) THEN
-            varout(i,j,:,:,:) = fillVal64
-          ELSE 
-            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))
-                    IF (mask(ix,iy,k,l,m) == 1) THEN
-                      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 IF
-                  END DO
-                END DO
-              END DO
-            END DO
-          END IF
-        END DO
-      END DO
-    CASE('npp')
-      DO i=1, pdimx
-        DO j=1, pdimy
-          ix = INT(outLlw(1,1,i,j))
-          iy = INT(outLlw(2,1,i,j))
-          IF ( (outLlw(3,1,i,j) == -1.) .OR. (mask(ix,iy,1,1,1) == 0) ) THEN
-            varout(i,j,:,:,:) = fillVal64
-          ELSE 
-            DO k=1, d3
-              DO l=1, d4
-                DO m=1, d5
-                  varout(i,j,k,l,m) = var5Din(ix,iy,k,l,m)*outLlw(3,1,i,j)
-                END DO
-              END DO
-            END DO
-          END IF
-        END DO
-      END DO
-  END SELECT
-
-END SUBROUTINE var5D_IntProj
-
-SUBROUTINE Interpolate(projlon, projlat, lonvs, latvs, mindiff, inpt, diffs, ilonlat, dimx, dimy,     &
-  Ninpts)
-! Subroutine which finds the closest grid point within a projection
-
-  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)               :: lonvs, latvs
-  REAL(r_k), INTENT(in)                                  :: mindiff
-  INTEGER, DIMENSION(Ninpts), INTENT(inout)              :: inpt
-  REAL(r_k), DIMENSION(Ninpts), INTENT(out)              :: diffs
-  INTEGER, DIMENSION(Ninpts,2), INTENT(out)              :: ilonlat
-
-! Local
-  INTEGER                                                :: iv
-  REAL(r_k)                                              :: mindiffLl
-  INTEGER                                                :: Ninpts1
-  REAL(r_k), DIMENSION(dimx,dimy)                        :: difflonlat
-  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
-  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
-! inpt: whether the point has already been localized
-! diffs: distance of point from the input data to the closest target point
-! ilonlat: longitude and latitude of the point
-! ncid: netCDF output file id
-
-  fname = 'Interpolate'
-  Ninpts1 = Ninpts/100
-
-  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
-  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
-
-  DO iv=1,Ninpts
-    IF (inpt(iv) <= 0) THEN
-! Not using the subroutine, not efficient!
-!      CALL lonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, lonvs(iv), latvs(iv),     &
-!        ilonlat, mindiffLl)
-
-      IF (lonvs(iv) < extremelon(1) .OR. lonvs(iv) > extremelon(2)) THEN
-        PRINT *, TRIM(ErrWarnMsg('err'))
-        PRINT *,'  ' // TRIM(fname) // ': longitude outside data range!!'
-        PRINT *,'    given value:', lonvs(iv),' outside (',extremelon(1),' ,',extremelon(2),' )'
-        STOP
-      END IF
-      IF (latvs(iv) < extremelat(1) .OR. latvs(iv) > extremelat(2)) THEN
-        PRINT *, TRIM(ErrWarnMsg('err'))
-        PRINT *,'  ' // TRIM(fname) // ': latitude outside data range!!'
-        PRINT *,'    given value:', latvs(iv),' outside (',extremelat(1),' ,',extremelat(2),' )'
-        STOP
-      END IF
-
-! Find point
-      difflonlat = SQRT((projlon-lonvs(iv))**2. + (projlat-latvs(iv))**2.)
-      mindiffLl = MINVAL(difflonlat)
-      ilonlat(iv,:) = index2DArrayR(difflonlat, dimx, dimy, mindiffLl)
-
-      IF (mindiffLl <= mindiff) THEN 
-!        percendone(iv,Ninpts,0.5,'done:')
-
-        IF (ilonlat(iv,1) >= 0 .AND. ilonlat(iv,2) >= 0) THEN
-          diffs(iv) = mindiffLl
-          inpt(iv) = 1
-!          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(iv,1), ' ilat:', ilonlat(iv,2)
-          STOP
-        END IF
-
-!        IF (MOD(iv,Ninpts1) == 0) newnc.sync()
-      ELSE
-! Because doing boxes and Goode is not conitnuos, we should jump this error message
-        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 Interpolate
-
-SUBROUTINE Interpolate1DLl(projlon, projlat, lonvs, latvs, mindiff, inpt, diffs, ilonlat, dimx, dimy, &
-  Ninpts)
-! Subroutine which finds the closest grid point within a projection with 1D longitudes and latitudes
-
-  IMPLICIT NONE
-
-!  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
-  INTEGER, INTENT(in)                                    :: dimx, dimy
-  REAL(r_k), DIMENSION(dimx), INTENT(in)                 :: projlon
-  REAL(r_k), DIMENSION(dimy), INTENT(in)                 :: projlat
-  INTEGER, INTENT(in)                                    :: Ninpts
-  REAL(r_k), DIMENSION(Ninpts), INTENT(in)               :: lonvs, latvs
-  REAL(r_k), INTENT(in)                                  :: mindiff
-  INTEGER, DIMENSION(Ninpts), INTENT(inout)              :: inpt
-  REAL(r_k), DIMENSION(Ninpts), INTENT(out)              :: diffs
-  INTEGER, DIMENSION(Ninpts,2), INTENT(out)              :: ilonlat
-
-! Local
-  INTEGER                                                :: iv
-  REAL(r_k)                                              :: mindifflo, mindiffLa, mindiffLl
-  INTEGER                                                :: Ninpts1
-  REAL(r_k), DIMENSION(dimx)                             :: difflon
-  REAL(r_k), DIMENSION(dimy)                             :: difflat
-  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
-  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
-! inpt: whether the point has already been localized
-! diffs: distance of point from the input data to the closest target point
-! ilonlat: longitude and latitude of the point
-! ncid: netCDF output file id
-
-  fname = 'Interpolate1DLl'
-  Ninpts1 = Ninpts/100
-
-  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
-  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
-
-  DO iv=1,Ninpts
-    IF (inpt(iv) <= 0) THEN
-! Not using the subroutine, not efficient!
-!      CALL lonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, lonvs(iv), latvs(iv),     &
-!        ilonlat, mindiffLl)
-
-      IF (lonvs(iv) < extremelon(1) .OR. lonvs(iv) > extremelon(2)) THEN
-        PRINT *, TRIM(ErrWarnMsg('err'))
-        PRINT *,'  ' // TRIM(fname) // ': longitude outside data range!!'
-        PRINT *,'    given value:', lonvs(iv),' outside (',extremelon(1),' ,',extremelon(2),' )'
-        STOP
-      END IF
-      IF (latvs(iv) < extremelat(1) .OR. latvs(iv) > extremelat(2)) THEN
-        PRINT *, TRIM(ErrWarnMsg('err'))
-        PRINT *,'  ' // TRIM(fname) // ': latitude outside data range!!'
-        PRINT *,'    given value:', latvs(iv),' outside (',extremelat(1),' ,',extremelat(2),' )'
-        STOP
-      END IF
-
-! Find point
-      difflon = SQRT((projlon-lonvs(iv))**2.)
-      difflat = SQRT((projlat-latvs(iv))**2.)
-      mindifflo = MINVAL(difflon)
-      mindiffLa = MINVAL(difflat)
-      mindifflL = SQRT(mindifflo*mindifflo + mindiffLa*mindiffLa)
-      ilonlat(iv,1) = index1DArrayR(difflon, dimx, mindifflo)
-      ilonlat(iv,2) = index1DArrayR(difflat, dimy, mindiffLa)
-!      PRINT *,'  Lluis: iv',iv,' lonvs:', lonvs(iv),' latvs:',latvs(iv)
-!      PRINT *,'  Lluis: mindifflo:', mindifflo,' ilonlat(1):',ilonlat(iv,1)
-!      PRINT *,'  Lluis: mindiffLa:', mindiffLa,' ilonlat(2):',ilonlat(iv,2)
-
-
-      IF (mindiffLl <= mindiff) THEN 
-!        percendone(iv,Ninpts,0.5,'done:')
-
-        IF (ilonlat(iv,1) >= 1 .AND. ilonlat(iv,2) >= 1) THEN
-          diffs(iv) = mindiffLl
-          inpt(iv) = 1
-!          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(iv,1), ' ilat:', ilonlat(iv,2)
-          STOP
-        END IF
-
-!        IF (MOD(iv,Ninpts1) == 0) newnc.sync()
-      ELSE
-! Because doing boxes and Goode is not conitnuos, we should jump this error message
-        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 Interpolate1DLl
-
-
- SUBROUTINE interp (data_in, pres_field, interp_levels, psfc, ter, tk, qv, LINLOG, extrapolate,       &
-                     GEOPT, MISSING, data_out, ix, iy, iz, it, num_metgrid_levels)
-! Interpolation subroutine from the p_interp.F90 NCAR program
-!  Program to read wrfout data and interpolate to pressure levels
-!  The program reads namelist.pinterp
-!  November 2007 - Cindy Bruyere
-!
-     INTEGER, INTENT(IN)                                 :: ix, iy, iz, it
-     INTEGER, INTENT(IN)                                 :: num_metgrid_levels, LINLOG
-     REAL, DIMENSION(ix,iy,iz,it), INTENT(IN)            :: data_in, pres_field, tk, qv
-     REAL, DIMENSION(ix,iy,it), INTENT(IN)               :: psfc
-     REAL, DIMENSION(ix,iy), INTENT(IN)                  :: ter
-     REAL, DIMENSION(num_metgrid_levels), INTENT(IN)     :: interp_levels
-     INTEGER, INTENT(IN)                                 :: extrapolate
-     REAL, INTENT(IN)                                    :: MISSING
-     LOGICAL, INTENT(IN)                                 :: GEOPT
-     REAL, DIMENSION(ix,iy,num_metgrid_levels,it),                                                    &
-       INTENT(OUT)                                       :: data_out
-
-! Local
-     INTEGER                                             :: i, j, itt, k, kk, kin
-     REAL, DIMENSION(num_metgrid_levels)                 :: data_out1D
-     REAL, DIMENSION(iz)                                 :: data_in1D, pres_field1D
-     REAL, DIMENSION(ix, iy, num_metgrid_levels, it)     :: N
-     REAL                                                :: sumA, sumN, AVE_geopt
-
-!!!!!!! Variables
-! data_out: interpolated field
-! data_in: field to interpolate
-! pres_field: pressure field [Pa]
-! interp_levels: pressure levels to interpolate [hPa]
-! psfc: surface pressure [Pa]
-! ter: terrein height [m]
-! tk: temperature [K]
-! qv: mositure mizing ratio [kg/kg]
-! i[x/y/z/t]: size of the matrices
-! num_metgrid_levels: number of pressure values to interpolate
-! LINLOG: if abs(linlog)=1 use linear interp in pressure
-!         if abs(linlog)=2 linear interp in ln(pressure)
-! extrapolate: whether to set to missing value below/above model ground and top (0), or extrapolate (1)
-! GEOPT: Wether the file is the geopotential file or not
-! MISSING: Missing value
-
-     N = 1.0
-
-     expon=287.04*.0065/9.81
-
-     do itt = 1, it
-        do j = 1, iy
-        do i = 1, ix
-           data_in1D(:)    = data_in(i,j,:,itt)
-           pres_field1D(:) = pres_field(i,j,:,itt)
-           CALL int1D (data_out1D, data_in1D, pres_field1D, interp_levels, iz, num_metgrid_levels, LINLOG, MISSING)
-           data_out(i,j,:,itt) = data_out1D(:)
-        end do
-        end do
-     end do
-
-
-     ! Fill in missing values
-     IF ( extrapolate == 0 ) RETURN       !! no extrapolation - we are out of here
-
-     ! First find where about 400 hPa is located
-     kk = 0
-     find_kk : do k = 1, num_metgrid_levels
-        kk = k
-        if ( interp_levels(k) <= 40000. ) exit find_kk
-     end do find_kk
-
-     
-     IF ( GEOPT ) THEN     !! geopt is treated different below ground
-
-        do itt = 1, it
-           do k = 1, kk
-              do j = 1, iy
-              do i = 1, ix
-                 IF ( data_out(i,j,k,itt) == MISSING .AND. interp_levels(k) < psfc(i,j,itt) ) THEN
-
-!                We are below the first model level, but above the ground 
-
-                    data_out(i,j,k,itt) = ((interp_levels(k) - pres_field(i,j,1,itt))*ter(i,j)*9.81 +  &
-                                           (psfc(i,j,itt) - interp_levels(k))*data_in(i,j,1,itt) ) /   &
-                                          (psfc(i,j,itt) - pres_field(i,j,1,itt))
-
-                 ELSEIF ( data_out(i,j,k,itt) == MISSING ) THEN
-
-!                We are below both the ground and the lowest data level.
-
-!                First, find the model level that is closest to a "target" pressure
-!                level, where the "target" pressure is delta-p less that the local
-!                value of a horizontally smoothed surface pressure field.  We use
-!                delta-p = 150 hPa here. A standard lapse rate temperature profile
-!                passing through the temperature at this model level will be used
-!                to define the temperature profile below ground.  This is similar
-!                to the Benjamin and Miller (1990) method, except that for
-!                simplicity, they used 700 hPa everywhere for the "target" pressure.
-!                Code similar to what is implemented in RIP4
-
-                    ptarget = (psfc(i,j,itt)*.01) - 150.
-                    dpmin=1.e4
-                    kupper = 0
-                    loop_kIN : do kin=iz,1,-1
-                       kupper = kin
-                       dp=abs( (pres_field(i,j,kin,itt)*.01) - ptarget )
-                       if (dp.gt.dpmin) exit loop_kIN
-                       dpmin=min(dpmin,dp)
-                    enddo loop_kIN
-
-                    pbot=max(pres_field(i,j,1,itt),psfc(i,j,itt))
-                    zbot=min(data_in(i,j,1,itt)/9.81,ter(i,j))
-
-                    tbotextrap=tk(i,j,kupper,itt)*(pbot/pres_field(i,j,kupper,itt))**expon
-                    tvbotextrap=virtual(tbotextrap,qv(i,j,1,itt))
-
-                    data_out(i,j,k,itt) = (zbot+tvbotextrap/.0065*(1.-(interp_levels(k)/pbot)**expon))*9.81
-               
-                 ENDIF
-              enddo
-              enddo
-           enddo
-        enddo
-
-
-        !!! Code for filling missing data with an average - we don't want to do this
-        !!do itt = 1, it
-           !!loop_levels : do k = 1, num_metgrid_levels
-              !!sumA = SUM(data_out(:,:,k,itt), MASK = data_out(:,:,k,itt) /= MISSING)
-              !!sumN = SUM(N(:,:,k,itt), MASK = data_out(:,:,k,itt) /= MISSING)
-              !!IF ( sumN == 0. ) CYCLE loop_levels
-              !!AVE_geopt = sumA/sumN
-              !!WHERE ( data_out(:,:,k,itt) == MISSING )
-                 !!data_out(:,:,k,itt) = AVE_geopt
-              !!END WHERE
-           !!end do loop_levels
-        !!end do
-
-     END IF
-     
-     !!! All other fields and geopt at higher levels come here
-     do itt = 1, it
-        do j = 1, iy
-        do i = 1, ix
-          do k = 1, kk
-             if ( data_out(i,j,k,itt) == MISSING ) data_out(i,j,k,itt) = data_in(i,j,1,itt)
-          end do
-          do k = kk+1, num_metgrid_levels
-             if ( data_out(i,j,k,itt) == MISSING ) data_out(i,j,k,itt) = data_in(i,j,iz,itt)
-          end do
-        end do
-        end do
-     end do
-
- END SUBROUTINE interp 
-
- SUBROUTINE int1D(xxout, xxin, ppin, ppout, npin, npout, LINLOG, MISSING)
-
-! Modified from int2p - NCL code
-! routine to interpolate from one set of pressure levels
-! .   to another set  using linear or ln(p) interpolation
-!
-! NCL: xout = int2p (pin,xin,pout,linlog)
-! This code was originally written for a specific purpose.
-! .   Several features were added for incorporation into NCL's
-! .   function suite including linear extrapolation.
-!
-! nomenclature:
-!
-! .   ppin   - input pressure levels. The pin can be
-! .            be in ascending or descending order
-! .   xxin   - data at corresponding input pressure levels
-! .   npin   - number of input pressure levels >= 2
-! .   ppout  - output pressure levels (input by user)
-! .            same (ascending or descending) order as pin
-! .   xxout  - data at corresponding output pressure levels
-! .   npout  - number of output pressure levels
-! .   linlog - if abs(linlog)=1 use linear interp in pressure
-! .            if abs(linlog)=2 linear interp in ln(pressure)
-! .   missing- missing data code. 
-
-!                                                ! input types
-      INTEGER   :: npin,npout,linlog,ier
-      real      :: ppin(npin),xxin(npin),ppout(npout)
-      real      :: MISSING       
-     logical                                          :: AVERAGE
-!                                                ! output
-      real      :: xxout(npout)
-      INTEGER   :: j1,np,nl,nin,nlmax,nplvl
-      INTEGER   :: nlsave,np1,no1,n1,n2,nlstrt
-      real      :: slope,pa,pb,pc
-
-! automatic arrays
-      real      :: pin(npin),xin(npin),p(npin),x(npin)
-      real      :: pout(npout),xout(npout)
-
-
-      xxout = MISSING
-      pout  = ppout
-      p     = ppin
-      x     = xxin
-      nlmax = npin
-
-! exact p-level matches
-      nlstrt = 1
-      nlsave = 1
-      do np = 1,npout
-          xout(np) = MISSING
-          do nl = nlstrt,nlmax
-              if (pout(np).eq.p(nl)) then
-                  xout(np) = x(nl)
-                  nlsave = nl + 1
-                  go to 10
-              end if
-          end do
-   10     nlstrt = nlsave
-      end do
-
-      if (LINLOG.eq.1) then
-          do np = 1,npout
-              do nl = 1,nlmax - 1
-                  if (pout(np).lt.p(nl) .and. pout(np).gt.p(nl+1)) then
-                      slope = (x(nl)-x(nl+1))/ (p(nl)-p(nl+1))
-                      xout(np) = x(nl+1) + slope* (pout(np)-p(nl+1))
-                  end if
-              end do
-          end do
-      elseif (LINLOG.eq.2) then
-          do np = 1,npout
-              do nl = 1,nlmax - 1
-                  if (pout(np).lt.p(nl) .and. pout(np).gt.p(nl+1)) then
-                      pa = log(p(nl))
-                      pb = log(pout(np))
-! special case: in case someone inadvertently enter p=0.
-                      if (p(nl+1).gt.0.d0) then
-                          pc = log(p(nl+1))
-                      else
-                          pc = log(1.d-4)
-                      end if
-
-                      slope = (x(nl)-x(nl+1))/ (pa-pc)
-                      xout(np) = x(nl+1) + slope* (pb-pc)
-                  end if
-              end do
-          end do
-      end if
-
-
-! place results in the return array;
-      xxout = xout
-
- END SUBROUTINE int1D
-
- FUNCTION virtual (tmp,rmix)
-!      This function returns virtual temperature in K, given temperature
-!      in K and mixing ratio in kg/kg.
-
-     real                              :: tmp, rmix, virtual
-
-     virtual=tmp*(0.622+rmix)/(0.622*(1.+rmix))
-
- END FUNCTION virtual
-
-END MODULE module_ForInterpolate
Index: /trunk/tools/module_ForInterpolate.f90
===================================================================
--- /trunk/tools/module_ForInterpolate.f90	(revision 1606)
+++ /trunk/tools/module_ForInterpolate.f90	(revision 1606)
@@ -0,0 +1,1671 @@
+! Module to interpolate values from a giving projection and pressure interpolation
+! To be included in a python
+! 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
+
+  USE module_generic
+
+  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
+
+  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 CoarselonlatFindExact(dx, dy, ilon, ilat, nxlon, nxlat, fracx, fracy, fraclon, fraclat,    &
+  iv, lonv, latv, per, Nperx, Npery, mindiff, ilonlat, mindiffLl)
+! Function to search a given value from a coarser version of the data
+
+  IMPLICIT NONE
+
+!  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
+  INTEGER, INTENT(in)                                    :: dx, dy, iv
+  REAL(r_k), DIMENSION(dx,dy), INTENT(in)                :: ilon, ilat
+  INTEGER, INTENT(in)                                    :: fracx, fracy
+  REAL(r_k), DIMENSION(Nperx,Npery), INTENT(in)          :: fraclon, fraclat
+  REAL(r_k), INTENT(in)                                  :: lonv, latv, per, mindiff
+  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
+  INTEGER                                                :: i
+  REAL(r_k), DIMENSION(Nperx,Npery)                      :: difffraclonlat
+  REAL(r_k)                                              :: mindifffracLl
+  INTEGER, DIMENSION(2)                                  :: ilonlatfrac
+  INTEGER                                                :: ixbeg, ixend, iybeg, iyend
+  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
+! iv: point in the input data
+! 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
+! frac[x/y]: Number of grid points for each fraction
+! fraclon, fraclat: longitude and latitude fractional matricies to perform the first guess
+! mindiff: authorized minimal distance between input and interpolated point
+! ilonlat: grid point on the total lon,lat matrix
+! mindiffLl: distance between input and interpolated point
+
+  fname = 'CoarselonlatFindExact'
+
+  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
+
+! 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
+    PRINT *,'Lluis!',fraclonv, '>=', lonv,'&', fraclatv, '>=', latv
+    IF (ilonlatfrac(1) > 1) THEN
+      ixbeg = (ilonlatfrac(1)-1)*fracx
+      ixend = ilonlatfrac(1)*fracx+1
+    ELSE
+      PRINT *,'Lluis 2'
+      ixbeg = 1
+      ixend = fracx+1
+    END IF
+    IF (ilonlatfrac(2) > 1) THEN
+      iybeg = (ilonlatfrac(2)-1)*fracy
+      iyend = ilonlatfrac(2)*fracy+1
+    ELSE
+      iybeg = 1
+      iyend = fracy+1
+    END IF
+  ELSE IF (fraclonv < lonv .AND. fraclatv >= latv) THEN
+    PRINT *,'Lluis!',fraclonv, '<', lonv,'&', fraclatv, '>=', latv
+    IF (ilonlatfrac(1) < Nperx) THEN
+      PRINT *,'Lluis 2'
+      IF (ilonlatfrac(1) /= 1) THEN
+        ixbeg = (ilonlatfrac(1)-1)*fracx
+        ixend = ilonlatfrac(1)*fracx+1
+      ELSE
+        ixbeg = 1
+        ixend = fracx+1
+      END IF
+    ELSE
+      ixbeg = Nperx*fracx
+      ixend = dx+1
+    END IF
+    IF (ilonlatfrac(2) > 1) THEN
+      iybeg = (ilonlatfrac(2)-1)*fracy
+      iyend = ilonlatfrac(2)*fracy+1
+    ELSE
+      iybeg = 1
+      iyend = fracy+1
+    END IF    
+  ELSE IF (fraclonv < lonv .AND. fraclatv < latv) THEN
+    PRINT *,'Lluis!',fraclonv, '<', lonv,'&', fraclatv, '<', latv
+    IF (ilonlatfrac(1) < Nperx) THEN
+      IF (ilonlatfrac(1) /= 1) THEN
+        ixbeg = (ilonlatfrac(1)-1)*fracx
+        ixend = ilonlatfrac(1)*fracx+1
+      ELSE
+        ixbeg = 1
+        ixend = fracx+1
+      END IF
+    ELSE
+      ixbeg = Nperx*fracx
+      ixend = dx+1
+    ENDIF
+    IF (ilonlatfrac(2) < Npery) THEN
+      IF (ilonlatfrac(2) /= 1) THEN
+        iybeg = (ilonlatfrac(2)-1)*fracy
+        iyend = ilonlatfrac(2)*fracy+1
+      ELSE
+        iybeg = 1
+        iyend = fracy+1
+      END IF
+    ELSE
+      iybeg = Npery*fracy
+      iyend = dy+1
+    END IF
+  ELSE IF (fraclonv >= lonv .AND. fraclatv < latv) THEN
+    PRINT *,'Llui!',fraclonv, '>=', lonv,'&', fraclatv, '<', latv
+    IF (ilonlatfrac(1) > 1) THEN
+      ixbeg = (ilonlatfrac(1)-1)*fracx
+      ixend = ilonlatfrac(1)*fracx+1
+    ELSE
+      ixbeg = 1
+      ixend = fracx+1
+    END IF
+    IF (ilonlatfrac(2) < Npery) THEN
+      IF (ilonlatfrac(2) /= 1) THEN
+        iybeg = (ilonlatfrac(2)-1)*fracy
+        iyend = ilonlatfrac(2)*fracy+1
+      ELSE
+        iybeg = 1
+        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)
+
+  IF (mindiffLl > mindiff) THEN 
+    difflonlat = SQRT((lon-lonv)**2. + (lat-latv)**2.)
+    mindiffLl = MINVAL(difflonlat)
+  END IF
+
+  IF (mindiffLl > mindiff) THEN
+    PRINT *,TRIM(ErrWarnMsg('err'))
+    PRINT *,'  ' // TRIM(fname) // ': not equivalent point closer than:',mindiff,' found!!'
+    PRINT *,'    at input point iv:', iv,' lon/lat:', lonv,', ',latv,' distance:',mindiffLl
+    PRINT *,'    Fraction values _______ (',Nperx,', ',Npery ,')'
+    PRINT *,'      fraclon'
+    DO i=1, Nperx
+      PRINT *,'      ',fraclon(i,:)
+    END DO
+    PRINT *,'      fraclat'
+    DO i=1, Nperx
+      PRINT *,'      ',fraclat(i,:)
+    END DO
+    PRINT *,'      frac lon, lat:', fraclon(ilonlatfrac(1),ilonlatfrac(2)), ' ,',                     &
+     fraclat(ilonlatfrac(1),ilonlatfrac(2))
+    PRINT *,'      mindifffracLl:', mindifffracLl, ' ilonlatfrac:', ilonlatfrac
+    PRINT *,'    Coarse values _______'
+    PRINT *,'      indices. x:', ixbeg, ', ', ixend, ' y:', iybeg, ', ', iyend
+    PRINT *,'      lon range:', '(',ilon(ixbeg,iybeg),', ',ilon(ixend,iyend),')'
+    PRINT *,'      lat range:', '(',ilat(ixbeg,iybeg),', ',ilat(ixend,iyend),')'
+    PRINT *,'      lon', UBOUND(lon)
+    DO i=1, ixend-ixbeg+1
+      PRINT *,'      ',lon(i,:)
+    END DO
+    PRINT *,'      lat', UBOUND(lat)
+    DO i=1, ixend-ixbeg+1
+      PRINT *,'      ',lat(i,:)
+    END DO
+    STOP
+  END IF
+
+  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 CoarselonlatFindExact
+
+SUBROUTINE lonlatFind(dx, dy, ilon, ilat, nxlon, nxlat, lonv, latv, ilonlat, mindiffLl)
+! Function to search a given value from a coarser version of the data
+
+  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), INTENT(in)                                  :: lonv, latv
+  REAL(r_k), DIMENSION(2), INTENT(in)                    :: nxlon, nxlat
+  INTEGER, DIMENSION(2), INTENT(out)                     :: ilonlat
+  REAL(r_k), INTENT(out)                                 :: mindiffLl
+! Local
+  REAL(r_k), DIMENSION(dx,dy)                            :: difflonlat
+  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
+
+  fname = 'lonlatFind'
+
+  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
+
+! Find point
+  difflonlat = SQRT((ilon-lonv)**2. + (ilat-latv)**2.)
+  mindiffLl = MINVAL(difflonlat)
+  ilonlat = index2DArrayR(difflonlat, dx, dy, mindiffLl)
+
+!  PRINT *,'mindiffLl:', mindiffLl, ' ilatlon:', ilatlon
+!  PRINT *,'lon, lat:', lon(ilonlat(1),ilonlat(2)), ' ,', lat(ilonlat(1),ilonlat(2))
+
+  RETURN
+
+END SUBROUTINE lonlatFind
+
+SUBROUTINE CoarseInterpolate(projlon, projlat, lonvs, latvs, percen, mindiff, inpt, ilonlat,          &
+  mindiffLl, dimx, dimy, Ninpts)
+! Subroutine which finds the closest grid point within a projection throughout a first guest 
+!   approche from percentages of the whole domain
+
+  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)               :: inpt, lonvs, latvs
+  REAL(r_k), INTENT(in)                                  :: mindiff, percen
+  INTEGER, DIMENSION(Ninpts,2), INTENT(out)              :: ilonlat
+  REAL(r_k), DIMENSION(Ninpts), INTENT(out)              :: mindiffLl
+
+! Local
+  INTEGER                                                :: iv,i,j
+  INTEGER                                                :: ierr
+  INTEGER                                                :: Ninpts1
+  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
+  REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: fractionlon, fractionlat
+  INTEGER                                                :: dfracdx, dfracdy, 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
+! inpt: whether the point has already been localized (1) or not (0)
+! ilonlat: Longitude and Latitude of the input points
+! mindiffLl: minimum difference between target and source longitude/latitude (in degrees)
+
+  fname = 'CoarseInterpolate'
+  Ninpts1 = Ninpts/100
+
+  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
+  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
+
+  PRINT *,'  ' // TRIM(fname) //' total space:', dimx, ', ', dimy, ' %', percen
+
+  dfracdx = INT(1./percen+1)
+  dfracdy = INT(1./percen+1)
+  fracdx = INT(dimx*percen)
+  fracdy = INT(dimy*percen)
+  PRINT *,'  ' // TRIM(fname) //' fraction:', dfracdx, ', ', dfracdy, ' freq:', fracdx,', ',fracdy
+
+  IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon)
+  ALLOCATE(fractionlon(dfracdx, dfracdy), STAT=ierr)
+  IF (ierr /= 0) THEN
+    PRINT *,TRIM(ErrWarnMsg('err'))
+    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlon' !!"
+    STOP
+  END IF
+  IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat)
+  ALLOCATE(fractionlat(dfracdx, dfracdy), STAT=ierr)
+  IF (ierr /= 0) THEN
+    PRINT *,TRIM(ErrWarnMsg('err'))
+    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlat' !!"
+    STOP
+  END IF
+
+  DO i=1,dfracdx
+    DO j=1,dfracdy
+      fractionlon(i,j) = projlon(fracdx*(i-1)+1,fracdy*(j-1)+1)
+      fractionlat(i,j) = projlat(fracdx*(i-1)+1,fracdy*(j-1)+1)
+!      PRINT *,'i,j:',i,', ',j,' frac ij:',fracdx*(i-1),', ',fracdy*(j-1),' lonlat:', fractionlon(i,j),&
+!        ', ',fractionlat(i,j)
+    END DO
+  END DO
+
+!  PRINT *,'  ' // TRIM(fname) // ' fractions of:'
+!  PRINT *,' lon _______ (',dfracdx,', ',dfracdy,')'
+!  DO i=1,dfracdx
+!    PRINT *,fractionlon(i,:)
+!  END DO
+!  PRINT *,' lat_______'
+!  DO i=1,dfracdx
+!    PRINT *,fractionlat(i,:)
+!  END DO
+
+  DO iv=1,Ninpts
+    IF (inpt(iv) == 0) THEN
+      CALL CoarselonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, fractionlon,        &
+        fractionlat, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, ilonlat(iv,:), mindiffLl(iv))
+
+      IF ((mindiffLl(iv) <= mindiff) .AND. .NOT.(ilonlat(iv,1) >= 0 .AND. ilonlat(iv,1) >= 0)) THEN
+        PRINT *,TRIM(ErrWarnMsg('err'))
+        PRINT *,'  ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv),        &
+          ' not relocated !!'
+        PRINT *,'    mindiffl:', mindiffLl(iv), ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2)
+        STOP
+      END IF
+
+    END IF
+  END DO
+
+END SUBROUTINE CoarseInterpolate
+
+SUBROUTINE CoarseInterpolateExact(projlon, projlat, lonvs, latvs, percen, mindiff, ivar, newvar,      &
+  newvarin, newvarinpt, newvarindiff, dimx, dimy, Ninpts)
+! Subroutine which finds the closest grid point within a projection throughout a first guest 
+!   and then whole domain approche from percentages of the whole domain 
+
+  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), INTENT(in)                                  :: 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,i,j
+  INTEGER                                                :: ierr
+  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                                                :: dfracdx, dfracdy, 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 = 'CoarseInterpolateExact'
+  Ninpts1 = Ninpts/100
+
+  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
+  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
+
+  PRINT *,'  ' // TRIM(fname) //' total space:', dimx, ', ', dimy, ' %', percen
+
+  dfracdx = INT(1./percen+1)
+  dfracdy = INT(1./percen+1)
+  fracdx = INT(dimx*percen)
+  fracdy = INT(dimy*percen)
+  PRINT *,'  ' // TRIM(fname) //' fraction:', dfracdx, ', ', dfracdy, ' freq:', fracdx,', ',fracdy
+
+  IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon)
+  ALLOCATE(fractionlon(dfracdx, dfracdy), STAT=ierr)
+  IF (ierr /= 0) THEN
+    PRINT *,TRIM(ErrWarnMsg('err'))
+    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlon' !!"
+    STOP
+  END IF
+  IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat)
+  ALLOCATE(fractionlat(dfracdx, dfracdy), STAT=ierr)
+  IF (ierr /= 0) THEN
+    PRINT *,TRIM(ErrWarnMsg('err'))
+    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlat' !!"
+    STOP
+  END IF
+
+  DO i=1,dfracdx
+    DO j=1,dfracdy
+      fractionlon(i,j) = projlon(fracdx*(i-1)+1,fracdy*(j-1)+1)
+      fractionlat(i,j) = projlat(fracdx*(i-1)+1,fracdy*(j-1)+1)
+!      PRINT *,'i,j:',i,', ',j,' frac ij:',fracdx*(i-1),', ',fracdy*(j-1),' lonlat:', fractionlon(i,j),&
+!        ', ',fractionlat(i,j)
+    END DO
+  END DO
+
+!  PRINT *,'  ' // TRIM(fname) // ' fractions of:'
+!  PRINT *,' lon _______ (',dfracdx,', ',dfracdy,')'
+!  DO i=1,dfracdx
+!    PRINT *,fractionlon(i,:)
+!  END DO
+!  PRINT *,' lat_______'
+!  DO i=1,dfracdx
+!    PRINT *,fractionlat(i,:)
+!  END DO
+
+  DO iv=1,Ninpts
+    IF (newvarinpt(iv) == 0) THEN
+      CALL CoarselonlatFindExact(dimx, dimy, projlon, projlat, extremelon, extremelat, fracdx, fracdy,&
+        fractionlon, fractionlat, iv, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, mindiff,        &
+        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 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
+
+  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)                                              :: inclon, inclat, maxdiffprojlonlat,        &
+    maxdiffinlonlat
+  REAL(r_k), DIMENSION(idimx,idimy)                      :: difflonlat
+  REAL(r_k), DIMENSION(idimx,idimy)                      :: idifflon, idifflat
+  REAL(r_k), DIMENSION(pdimx,pdimy)                      :: difflon, difflat
+  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat, ipos
+  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) /)
+ 
+  ! Maximum distance between grid points in input projection
+  idifflon = 0.
+  idifflat = 0.
+  idifflon(1:idimx-1,:) = inlonv(2:idimx,:)-inlonv(1:idimx-1,:)
+  idifflat(:,1:idimy-1) = inlatv(:,2:idimy)-inlatv(:,1:idimy-1)
+  maxdiffinlonlat = MAXVAL(SQRT(idifflon**2. + idifflat**2.))
+  ! Maximum distance between grid points in target projection
+  difflon = 0.
+  difflat = 0.
+  difflon(1:pdimx-1,:) = projlon(2:pdimx,:)-projlon(1:pdimx-1,:)
+  difflat(:,1:pdimy-1) = projlat(:,2:pdimy)-projlat(:,1:pdimy-1)
+  maxdiffprojlonlat = MAXVAL(SQRT(difflon**2. + difflat**2.))
+
+  IF (maxdiffinlonlat > maxdiffprojlonlat) THEN
+    PRINT *,TRIM(warnmsg)
+    PRINT *,'  ' //TRIM(fname)// '; input resolution: ', maxdiffinlonlat, ' is coarser than target:', &
+      maxdiffprojlonlat, ' !!'
+  END IF
+
+  ! Using case outside loop to be more efficient
+  SELECT CASE(TRIM(intkind))
+    CASE ('dis')
+      inclon = inlonv(2,1) - inlonv(1,1)
+      inclat = inlatv(1,2) - inlatv(1,1)
+
+      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)
+          IF ( (mindiffLl > maxdiffprojlonlat) .AND. (mindiffLl > maxdiffinlonlat)) THEN
+            outLlw(3,:,i,j) = 0.
+            outLlw(3,:,i,j) = -1.
+          ELSE
+            ! Getting the four surrounding values
+            iLl = index2DArrayR(difflonlat, idimx, idimy, mindiffLl)
+            IF ( (projlon(i,j) < inlonv(iLl(1),iLl(2)) .AND. inclon > 0.) .OR.                          &
+              (projlon(i,j) > inlonv(iLl(1),iLl(2)) .AND. inclon < 0.) ) 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) = MAX(iLl(1),1)
+              outLlw(1,2,i,j) = MAX(iLl(1),1)
+              outLlw(1,3,i,j) = MIN(iLl(1)+1,idimx)
+              outLlw(1,4,i,j) = MIN(iLl(1)+1,idimx)
+            END IF
+            IF ( (projlat(i,j) < inlatv(iLl(2),iLl(2)) .AND. inclat > 0.) .OR.                          &
+              (projlat(i,j) > inlatv(iLl(2),iLl(2)) .AND. inclat < 0.) )  THEN
+              outLlw(2,1,i,j) = MAX(iLl(2)-1,1)
+              outLlw(2,2,i,j) = MIN(iLl(2),idimy)
+              outLlw(2,3,i,j) = MIN(iLl(2),idimy)
+              outLlw(2,4,i,j) = MAX(iLl(2)-1,1)
+            ELSE
+              outLlw(2,1,i,j) = MAX(iLl(2),1)
+              outLlw(2,2,i,j) = MIN(iLl(2)+1,idimy)
+              outLlw(2,3,i,j) = MIN(iLl(2)+1,idimy)
+              outLlw(2,4,i,j) = MAX(iLl(2),1)
+            END IF
+            ! Computing distances
+!Keep the print for future checks?
+!            IF (MOD(i+j,200) == 0) THEN
+!              PRINT *,'center point:',i,j,'=', iLl,':',projlon(i,j),projlat(i,j),' incs',inclon,' ,',inclat
+!              PRINT *,'SW:', outLlw(1,1,i,j), outLlw(2,1,i,j),':',inlonv(outLlw(1,1,i,j), outLlw(2,1,i,j)),&
+!                inlatv(outLlw(1,1,i,j), outLlw(2,1,i,j))
+!              PRINT *,'NW:', outLlw(1,2,i,j), outLlw(2,2,i,j),':',inlonv(outLlw(1,2,i,j), outLlw(2,2,i,j)),&
+!                inlatv(outLlw(1,2,i,j), outLlw(2,2,i,j))
+!              PRINT *,'NE:', outLlw(1,3,i,j), outLlw(2,3,i,j),':',inlonv(outLlw(1,3,i,j), outLlw(2,3,i,j)),&
+!                inlatv(outLlw(1,3,i,j), outLlw(2,3,i,j))
+!              PRINT *,'SE:', outLlw(1,4,i,j), outLlw(2,4,i,j),':',inlonv(outLlw(1,4,i,j), outLlw(2,4,i,j)),&
+!                inlatv(outLlw(1,4,i,j), outLlw(2,4,i,j))
+!            END IF
+            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
+!              IF (i+j == 2) PRINT *,'iv:',iv,'dist:',dist,'weight:',outLlw(3,iv,i,j)
+            END DO
+            ! Normalizing
+            outLlw(3,:,i,j) = outLlw(3,:,i,j)/(SUM(outLlw(3,:,i,j)))
+!             IF (i+j == 2) PRINT *,'Normalized weights:',outLlw(3,:,i,j),':',SUM(outLlw(3,:,i,j))
+          END IF
+        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)
+          ipos = index2DArrayR(difflonlat, idimx, idimy, mindiffLl)
+          outLlw(1,1,i,j) = REAL(ipos(1))
+          outLlw(2,1,i,j) = REAL(ipos(2))
+          ! We do not want that values larger that the maximum distance between target grid points
+!          PRINT *,i,j,':',mindiffLl,'maxdiffLl:',maxdiffprojlonlat
+          IF ((mindiffLl .gt. maxdiffprojlonlat) .AND. (mindiffLl > maxdiffinlonlat)) THEN
+!            PRINT *,'  ' // TRIM(fname) // ': reprojected minimum distance to nearest grid point:',   &
+!              mindiffLl, ' larger than the maximum distance between grid points in target projection!!'
+            outLlw(3,1,i,j) = -1.
+          ELSE
+            outLlw(3,1,i,j) = 1.
+          END IF
+          ix = INT(outLlw(1,1,i,j))
+          iy = INT(outLlw(2,1,i,j))
+        END DO
+      END DO
+  END SELECT
+
+END SUBROUTINE LlInterpolateProjection
+
+SUBROUTINE var2D_IntProj(var2Din, inlonv, inlatv, projlon, projlat, intkind, mask, varout, idimx,     &
+  idimy, pdimx, pdimy)
+! Subroutine to interpolate a 2D variable
+
+  IMPLICIT NONE
+
+  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(idimx,idimy), INTENT(in)          :: var2Din
+  INTEGER, DIMENSION(idimx,idimy), INTENT(in)            :: mask
+  REAL(r_k), DIMENSION(pdimx,pdimy), 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.)
+! var2Din: 2D variable to interpolate
+! mask: mask of the intpu values (1: good, 0: none)
+! 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, pdimy
+          IF (outLlw(3,1,i,j) == -1.) THEN
+            varout(i,j) = fillVal64
+          ELSE 
+            varout(i,j) = 0.
+            DO iv=1, 4
+              ix = INT(outLlw(1,iv,i,j))
+              iy = INT(outLlw(2,iv,i,j))
+              IF (mask(ix,iy) == 1) THEN
+                w =  outLlw(3,iv,i,j)
+                varout(i,j) = varout(i,j) + w*var2Din(ix,iy)
+              END IF
+            END DO
+          END IF
+        END DO
+      END DO
+    CASE('npp')
+      DO i=1, pdimx
+        DO j=1, pdimy
+          ix = INT(outLlw(1,1,i,j))
+          iy = INT(outLlw(2,1,i,j))
+          IF ( (outLlw(3,1,i,j) == -1.) .OR. (mask(ix,iy) == 0) ) THEN
+            varout(i,j) = fillVal64
+          ELSE 
+            varout(i,j) = var2Din(ix,iy)*outLlw(3,1,i,j)
+          END IF
+        END DO
+      END DO
+  END SELECT
+
+END SUBROUTINE var2D_IntProj
+
+
+SUBROUTINE var3D_IntProj(var3Din, inlonv, inlatv, projlon, projlat, intkind, mask, varout, idimx,     &
+  idimy, pdimx, pdimy, d3)
+! Subroutine to interpolate a 3D variable
+
+  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(idimx,idimy,d3), INTENT(in)       :: var3Din
+  INTEGER, DIMENSION(idimx,idimy,d3), INTENT(in)         :: mask
+  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
+! mask: mask of the intpu values (1: good, 0: none)
+! 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, pdimy
+          IF (ALL(outLlw(3,:,i,j) == -1.)) THEN
+            varout(i,j,:) = fillVal64
+          ELSE 
+            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))
+                IF (mask(ix,iy,k) == 1) THEN
+                  w =  outLlw(3,iv,i,j)
+                  varout(i,j,k) = varout(i,j,k) + w*var3Din(ix,iy,k)
+                END IF
+              END DO
+            END DO
+          END IF
+        END DO
+      END DO
+    CASE('npp')
+      DO i=1, pdimx
+        DO j=1, pdimy
+          ix = INT(outLlw(1,1,i,j))
+          iy = INT(outLlw(2,1,i,j))
+          IF ( (outLlw(3,1,i,j) == -1.) .OR. (mask(ix,iy,1) == 0) ) THEN
+            varout(i,j,:) = fillVal64
+          ELSE 
+            DO k=1, d3
+              varout(i,j,k) = var3Din(ix,iy,k)*outLlw(3,1,i,j)
+            END DO
+          END IF
+        END DO
+      END DO
+  END SELECT
+
+END SUBROUTINE var3D_IntProj
+
+SUBROUTINE var4D_IntProj(var4Din, inlonv, inlatv, projlon, projlat, intkind, mask, varout, idimx,     &
+  idimy, pdimx, pdimy, d3, d4)
+! Subroutine to interpolate a 4D variable
+
+  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(idimx,idimy,d3,d4), INTENT(in)    :: var4Din
+  INTEGER, DIMENSION(idimx,idimy,d3,d4), INTENT(in)      :: mask
+  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
+! mask: mask of the intpu values (1: good, 0: none)
+! 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, pdimy
+          IF (ALL(outLlw(3,:,i,j) == -1.)) THEN
+            varout(i,j,:,:) = fillVal64
+          ELSE 
+            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))
+                  IF (mask(ix,iy,k,l) == 1) THEN
+                    w =  outLlw(3,iv,i,j)
+                    varout(i,j,k,l) = varout(i,j,k,l) + w*var4Din(ix,iy,k,l)
+                  END IF
+                END DO
+              END DO
+            END DO
+          END IF
+        END DO
+      END DO
+    CASE('npp')
+      DO i=1, pdimx
+        DO j=1, pdimy
+          ix = INT(outLlw(1,1,i,j))
+          iy = INT(outLlw(2,1,i,j))
+          IF ( (outLlw(3,1,i,j) == -1.) .OR. (mask(ix,iy,1,1) == 0) ) THEN
+            varout(i,j,:,:) = fillVal64
+          ELSE 
+            DO k=1, d3
+              DO l=1, d4
+                varout(i,j,k,l) = var4Din(ix,iy,k,l)*outLlw(3,1,i,j)
+              END DO
+            END DO
+          END IF
+        END DO
+      END DO
+  END SELECT
+
+END SUBROUTINE var4D_IntProj
+
+SUBROUTINE var5D_IntProj(var5Din, inlonv, inlatv, projlon, projlat, intkind, mask, varout, idimx,     &
+  idimy, pdimx, pdimy, d3, d4, d5)
+! Subroutine to interpolate a 5D variable
+
+  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(idimx,idimy,d3,d4,d5), INTENT(in) :: var5Din
+  INTEGER, DIMENSION(idimx,idimy,d3,d4,d5), INTENT(in)   :: mask
+  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
+! mask: mask of the intpu values (1: good, 0: none)
+! 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, pdimy
+          IF (ALL(outLlw(3,:,i,j) == -1.)) THEN
+            varout(i,j,:,:,:) = fillVal64
+          ELSE 
+            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))
+                    IF (mask(ix,iy,k,l,m) == 1) THEN
+                      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 IF
+                  END DO
+                END DO
+              END DO
+            END DO
+          END IF
+        END DO
+      END DO
+    CASE('npp')
+      DO i=1, pdimx
+        DO j=1, pdimy
+          ix = INT(outLlw(1,1,i,j))
+          iy = INT(outLlw(2,1,i,j))
+          IF ( (outLlw(3,1,i,j) == -1.) .OR. (mask(ix,iy,1,1,1) == 0) ) THEN
+            varout(i,j,:,:,:) = fillVal64
+          ELSE 
+            DO k=1, d3
+              DO l=1, d4
+                DO m=1, d5
+                  varout(i,j,k,l,m) = var5Din(ix,iy,k,l,m)*outLlw(3,1,i,j)
+                END DO
+              END DO
+            END DO
+          END IF
+        END DO
+      END DO
+  END SELECT
+
+END SUBROUTINE var5D_IntProj
+
+SUBROUTINE Interpolate(projlon, projlat, lonvs, latvs, mindiff, inpt, diffs, ilonlat, dimx, dimy,     &
+  Ninpts)
+! Subroutine which finds the closest grid point within a projection
+
+  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)               :: lonvs, latvs
+  REAL(r_k), INTENT(in)                                  :: mindiff
+  INTEGER, DIMENSION(Ninpts), INTENT(inout)              :: inpt
+  REAL(r_k), DIMENSION(Ninpts), INTENT(out)              :: diffs
+  INTEGER, DIMENSION(Ninpts,2), INTENT(out)              :: ilonlat
+
+! Local
+  INTEGER                                                :: iv
+  REAL(r_k)                                              :: mindiffLl
+  INTEGER                                                :: Ninpts1
+  REAL(r_k), DIMENSION(dimx,dimy)                        :: difflonlat
+  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
+  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
+! inpt: whether the point has already been localized
+! diffs: distance of point from the input data to the closest target point
+! ilonlat: longitude and latitude of the point
+! ncid: netCDF output file id
+
+  fname = 'Interpolate'
+  Ninpts1 = Ninpts/100
+
+  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
+  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
+
+  DO iv=1,Ninpts
+    IF (inpt(iv) <= 0) THEN
+! Not using the subroutine, not efficient!
+!      CALL lonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, lonvs(iv), latvs(iv),     &
+!        ilonlat, mindiffLl)
+
+      IF (lonvs(iv) < extremelon(1) .OR. lonvs(iv) > extremelon(2)) THEN
+        PRINT *, TRIM(ErrWarnMsg('err'))
+        PRINT *,'  ' // TRIM(fname) // ': longitude outside data range!!'
+        PRINT *,'    given value:', lonvs(iv),' outside (',extremelon(1),' ,',extremelon(2),' )'
+        STOP
+      END IF
+      IF (latvs(iv) < extremelat(1) .OR. latvs(iv) > extremelat(2)) THEN
+        PRINT *, TRIM(ErrWarnMsg('err'))
+        PRINT *,'  ' // TRIM(fname) // ': latitude outside data range!!'
+        PRINT *,'    given value:', latvs(iv),' outside (',extremelat(1),' ,',extremelat(2),' )'
+        STOP
+      END IF
+
+! Find point
+      difflonlat = SQRT((projlon-lonvs(iv))**2. + (projlat-latvs(iv))**2.)
+      mindiffLl = MINVAL(difflonlat)
+      ilonlat(iv,:) = index2DArrayR(difflonlat, dimx, dimy, mindiffLl)
+
+      IF (mindiffLl <= mindiff) THEN 
+!        percendone(iv,Ninpts,0.5,'done:')
+
+        IF (ilonlat(iv,1) >= 0 .AND. ilonlat(iv,2) >= 0) THEN
+          diffs(iv) = mindiffLl
+          inpt(iv) = 1
+!          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(iv,1), ' ilat:', ilonlat(iv,2)
+          STOP
+        END IF
+
+!        IF (MOD(iv,Ninpts1) == 0) newnc.sync()
+      ELSE
+! Because doing boxes and Goode is not conitnuos, we should jump this error message
+        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 Interpolate
+
+SUBROUTINE Interpolate1DLl(projlon, projlat, lonvs, latvs, mindiff, inpt, diffs, ilonlat, dimx, dimy, &
+  Ninpts)
+! Subroutine which finds the closest grid point within a projection with 1D longitudes and latitudes
+
+  IMPLICIT NONE
+
+!  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
+  INTEGER, INTENT(in)                                    :: dimx, dimy
+  REAL(r_k), DIMENSION(dimx), INTENT(in)                 :: projlon
+  REAL(r_k), DIMENSION(dimy), INTENT(in)                 :: projlat
+  INTEGER, INTENT(in)                                    :: Ninpts
+  REAL(r_k), DIMENSION(Ninpts), INTENT(in)               :: lonvs, latvs
+  REAL(r_k), INTENT(in)                                  :: mindiff
+  INTEGER, DIMENSION(Ninpts), INTENT(inout)              :: inpt
+  REAL(r_k), DIMENSION(Ninpts), INTENT(out)              :: diffs
+  INTEGER, DIMENSION(Ninpts,2), INTENT(out)              :: ilonlat
+
+! Local
+  INTEGER                                                :: iv
+  REAL(r_k)                                              :: mindifflo, mindiffLa, mindiffLl
+  INTEGER                                                :: Ninpts1
+  REAL(r_k), DIMENSION(dimx)                             :: difflon
+  REAL(r_k), DIMENSION(dimy)                             :: difflat
+  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
+  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
+! inpt: whether the point has already been localized
+! diffs: distance of point from the input data to the closest target point
+! ilonlat: longitude and latitude of the point
+! ncid: netCDF output file id
+
+  fname = 'Interpolate1DLl'
+  Ninpts1 = Ninpts/100
+
+  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
+  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
+
+  DO iv=1,Ninpts
+    IF (inpt(iv) <= 0) THEN
+! Not using the subroutine, not efficient!
+!      CALL lonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, lonvs(iv), latvs(iv),     &
+!        ilonlat, mindiffLl)
+
+      IF (lonvs(iv) < extremelon(1) .OR. lonvs(iv) > extremelon(2)) THEN
+        PRINT *, TRIM(ErrWarnMsg('err'))
+        PRINT *,'  ' // TRIM(fname) // ': longitude outside data range!!'
+        PRINT *,'    given value:', lonvs(iv),' outside (',extremelon(1),' ,',extremelon(2),' )'
+        STOP
+      END IF
+      IF (latvs(iv) < extremelat(1) .OR. latvs(iv) > extremelat(2)) THEN
+        PRINT *, TRIM(ErrWarnMsg('err'))
+        PRINT *,'  ' // TRIM(fname) // ': latitude outside data range!!'
+        PRINT *,'    given value:', latvs(iv),' outside (',extremelat(1),' ,',extremelat(2),' )'
+        STOP
+      END IF
+
+! Find point
+      difflon = SQRT((projlon-lonvs(iv))**2.)
+      difflat = SQRT((projlat-latvs(iv))**2.)
+      mindifflo = MINVAL(difflon)
+      mindiffLa = MINVAL(difflat)
+      mindifflL = SQRT(mindifflo*mindifflo + mindiffLa*mindiffLa)
+      ilonlat(iv,1) = index1DArrayR(difflon, dimx, mindifflo)
+      ilonlat(iv,2) = index1DArrayR(difflat, dimy, mindiffLa)
+!      PRINT *,'  Lluis: iv',iv,' lonvs:', lonvs(iv),' latvs:',latvs(iv)
+!      PRINT *,'  Lluis: mindifflo:', mindifflo,' ilonlat(1):',ilonlat(iv,1)
+!      PRINT *,'  Lluis: mindiffLa:', mindiffLa,' ilonlat(2):',ilonlat(iv,2)
+
+
+      IF (mindiffLl <= mindiff) THEN 
+!        percendone(iv,Ninpts,0.5,'done:')
+
+        IF (ilonlat(iv,1) >= 1 .AND. ilonlat(iv,2) >= 1) THEN
+          diffs(iv) = mindiffLl
+          inpt(iv) = 1
+!          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(iv,1), ' ilat:', ilonlat(iv,2)
+          STOP
+        END IF
+
+!        IF (MOD(iv,Ninpts1) == 0) newnc.sync()
+      ELSE
+! Because doing boxes and Goode is not conitnuos, we should jump this error message
+        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 Interpolate1DLl
+
+
+ SUBROUTINE interp (data_in, pres_field, interp_levels, psfc, ter, tk, qv, LINLOG, extrapolate,       &
+                     GEOPT, MISSING, data_out, ix, iy, iz, it, num_metgrid_levels)
+! Interpolation subroutine from the p_interp.F90 NCAR program
+!  Program to read wrfout data and interpolate to pressure levels
+!  The program reads namelist.pinterp
+!  November 2007 - Cindy Bruyere
+!
+     INTEGER, INTENT(IN)                                 :: ix, iy, iz, it
+     INTEGER, INTENT(IN)                                 :: num_metgrid_levels, LINLOG
+     REAL, DIMENSION(ix,iy,iz,it), INTENT(IN)            :: data_in, pres_field, tk, qv
+     REAL, DIMENSION(ix,iy,it), INTENT(IN)               :: psfc
+     REAL, DIMENSION(ix,iy), INTENT(IN)                  :: ter
+     REAL, DIMENSION(num_metgrid_levels), INTENT(IN)     :: interp_levels
+     INTEGER, INTENT(IN)                                 :: extrapolate
+     REAL, INTENT(IN)                                    :: MISSING
+     LOGICAL, INTENT(IN)                                 :: GEOPT
+     REAL, DIMENSION(ix,iy,num_metgrid_levels,it),                                                    &
+       INTENT(OUT)                                       :: data_out
+
+! Local
+     INTEGER                                             :: i, j, itt, k, kk, kin
+     REAL, DIMENSION(num_metgrid_levels)                 :: data_out1D
+     REAL, DIMENSION(iz)                                 :: data_in1D, pres_field1D
+     REAL, DIMENSION(ix, iy, num_metgrid_levels, it)     :: N
+     REAL                                                :: sumA, sumN, AVE_geopt
+
+!!!!!!! Variables
+! data_out: interpolated field
+! data_in: field to interpolate
+! pres_field: pressure field [Pa]
+! interp_levels: pressure levels to interpolate [hPa]
+! psfc: surface pressure [Pa]
+! ter: terrein height [m]
+! tk: temperature [K]
+! qv: mositure mizing ratio [kg/kg]
+! i[x/y/z/t]: size of the matrices
+! num_metgrid_levels: number of pressure values to interpolate
+! LINLOG: if abs(linlog)=1 use linear interp in pressure
+!         if abs(linlog)=2 linear interp in ln(pressure)
+! extrapolate: whether to set to missing value below/above model ground and top (0), or extrapolate (1)
+! GEOPT: Wether the file is the geopotential file or not
+! MISSING: Missing value
+
+     N = 1.0
+
+     expon=287.04*.0065/9.81
+
+     do itt = 1, it
+        do j = 1, iy
+        do i = 1, ix
+           data_in1D(:)    = data_in(i,j,:,itt)
+           pres_field1D(:) = pres_field(i,j,:,itt)
+           CALL int1D (data_out1D, data_in1D, pres_field1D, interp_levels, iz, num_metgrid_levels, LINLOG, MISSING)
+           data_out(i,j,:,itt) = data_out1D(:)
+        end do
+        end do
+     end do
+
+
+     ! Fill in missing values
+     IF ( extrapolate == 0 ) RETURN       !! no extrapolation - we are out of here
+
+     ! First find where about 400 hPa is located
+     kk = 0
+     find_kk : do k = 1, num_metgrid_levels
+        kk = k
+        if ( interp_levels(k) <= 40000. ) exit find_kk
+     end do find_kk
+
+     
+     IF ( GEOPT ) THEN     !! geopt is treated different below ground
+
+        do itt = 1, it
+           do k = 1, kk
+              do j = 1, iy
+              do i = 1, ix
+                 IF ( data_out(i,j,k,itt) == MISSING .AND. interp_levels(k) < psfc(i,j,itt) ) THEN
+
+!                We are below the first model level, but above the ground 
+
+                    data_out(i,j,k,itt) = ((interp_levels(k) - pres_field(i,j,1,itt))*ter(i,j)*9.81 +  &
+                                           (psfc(i,j,itt) - interp_levels(k))*data_in(i,j,1,itt) ) /   &
+                                          (psfc(i,j,itt) - pres_field(i,j,1,itt))
+
+                 ELSEIF ( data_out(i,j,k,itt) == MISSING ) THEN
+
+!                We are below both the ground and the lowest data level.
+
+!                First, find the model level that is closest to a "target" pressure
+!                level, where the "target" pressure is delta-p less that the local
+!                value of a horizontally smoothed surface pressure field.  We use
+!                delta-p = 150 hPa here. A standard lapse rate temperature profile
+!                passing through the temperature at this model level will be used
+!                to define the temperature profile below ground.  This is similar
+!                to the Benjamin and Miller (1990) method, except that for
+!                simplicity, they used 700 hPa everywhere for the "target" pressure.
+!                Code similar to what is implemented in RIP4
+
+                    ptarget = (psfc(i,j,itt)*.01) - 150.
+                    dpmin=1.e4
+                    kupper = 0
+                    loop_kIN : do kin=iz,1,-1
+                       kupper = kin
+                       dp=abs( (pres_field(i,j,kin,itt)*.01) - ptarget )
+                       if (dp.gt.dpmin) exit loop_kIN
+                       dpmin=min(dpmin,dp)
+                    enddo loop_kIN
+
+                    pbot=max(pres_field(i,j,1,itt),psfc(i,j,itt))
+                    zbot=min(data_in(i,j,1,itt)/9.81,ter(i,j))
+
+                    tbotextrap=tk(i,j,kupper,itt)*(pbot/pres_field(i,j,kupper,itt))**expon
+                    tvbotextrap=virtual(tbotextrap,qv(i,j,1,itt))
+
+                    data_out(i,j,k,itt) = (zbot+tvbotextrap/.0065*(1.-(interp_levels(k)/pbot)**expon))*9.81
+               
+                 ENDIF
+              enddo
+              enddo
+           enddo
+        enddo
+
+
+        !!! Code for filling missing data with an average - we don't want to do this
+        !!do itt = 1, it
+           !!loop_levels : do k = 1, num_metgrid_levels
+              !!sumA = SUM(data_out(:,:,k,itt), MASK = data_out(:,:,k,itt) /= MISSING)
+              !!sumN = SUM(N(:,:,k,itt), MASK = data_out(:,:,k,itt) /= MISSING)
+              !!IF ( sumN == 0. ) CYCLE loop_levels
+              !!AVE_geopt = sumA/sumN
+              !!WHERE ( data_out(:,:,k,itt) == MISSING )
+                 !!data_out(:,:,k,itt) = AVE_geopt
+              !!END WHERE
+           !!end do loop_levels
+        !!end do
+
+     END IF
+     
+     !!! All other fields and geopt at higher levels come here
+     do itt = 1, it
+        do j = 1, iy
+        do i = 1, ix
+          do k = 1, kk
+             if ( data_out(i,j,k,itt) == MISSING ) data_out(i,j,k,itt) = data_in(i,j,1,itt)
+          end do
+          do k = kk+1, num_metgrid_levels
+             if ( data_out(i,j,k,itt) == MISSING ) data_out(i,j,k,itt) = data_in(i,j,iz,itt)
+          end do
+        end do
+        end do
+     end do
+
+ END SUBROUTINE interp 
+
+ SUBROUTINE int1D(xxout, xxin, ppin, ppout, npin, npout, LINLOG, MISSING)
+
+! Modified from int2p - NCL code
+! routine to interpolate from one set of pressure levels
+! .   to another set  using linear or ln(p) interpolation
+!
+! NCL: xout = int2p (pin,xin,pout,linlog)
+! This code was originally written for a specific purpose.
+! .   Several features were added for incorporation into NCL's
+! .   function suite including linear extrapolation.
+!
+! nomenclature:
+!
+! .   ppin   - input pressure levels. The pin can be
+! .            be in ascending or descending order
+! .   xxin   - data at corresponding input pressure levels
+! .   npin   - number of input pressure levels >= 2
+! .   ppout  - output pressure levels (input by user)
+! .            same (ascending or descending) order as pin
+! .   xxout  - data at corresponding output pressure levels
+! .   npout  - number of output pressure levels
+! .   linlog - if abs(linlog)=1 use linear interp in pressure
+! .            if abs(linlog)=2 linear interp in ln(pressure)
+! .   missing- missing data code. 
+
+!                                                ! input types
+      INTEGER   :: npin,npout,linlog,ier
+      real      :: ppin(npin),xxin(npin),ppout(npout)
+      real      :: MISSING       
+     logical                                          :: AVERAGE
+!                                                ! output
+      real      :: xxout(npout)
+      INTEGER   :: j1,np,nl,nin,nlmax,nplvl
+      INTEGER   :: nlsave,np1,no1,n1,n2,nlstrt
+      real      :: slope,pa,pb,pc
+
+! automatic arrays
+      real      :: pin(npin),xin(npin),p(npin),x(npin)
+      real      :: pout(npout),xout(npout)
+
+
+      xxout = MISSING
+      pout  = ppout
+      p     = ppin
+      x     = xxin
+      nlmax = npin
+
+! exact p-level matches
+      nlstrt = 1
+      nlsave = 1
+      do np = 1,npout
+          xout(np) = MISSING
+          do nl = nlstrt,nlmax
+              if (pout(np).eq.p(nl)) then
+                  xout(np) = x(nl)
+                  nlsave = nl + 1
+                  go to 10
+              end if
+          end do
+   10     nlstrt = nlsave
+      end do
+
+      if (LINLOG.eq.1) then
+          do np = 1,npout
+              do nl = 1,nlmax - 1
+                  if (pout(np).lt.p(nl) .and. pout(np).gt.p(nl+1)) then
+                      slope = (x(nl)-x(nl+1))/ (p(nl)-p(nl+1))
+                      xout(np) = x(nl+1) + slope* (pout(np)-p(nl+1))
+                  end if
+              end do
+          end do
+      elseif (LINLOG.eq.2) then
+          do np = 1,npout
+              do nl = 1,nlmax - 1
+                  if (pout(np).lt.p(nl) .and. pout(np).gt.p(nl+1)) then
+                      pa = log(p(nl))
+                      pb = log(pout(np))
+! special case: in case someone inadvertently enter p=0.
+                      if (p(nl+1).gt.0.d0) then
+                          pc = log(p(nl+1))
+                      else
+                          pc = log(1.d-4)
+                      end if
+
+                      slope = (x(nl)-x(nl+1))/ (pa-pc)
+                      xout(np) = x(nl+1) + slope* (pb-pc)
+                  end if
+              end do
+          end do
+      end if
+
+
+! place results in the return array;
+      xxout = xout
+
+ END SUBROUTINE int1D
+
+ FUNCTION virtual (tmp,rmix)
+!      This function returns virtual temperature in K, given temperature
+!      in K and mixing ratio in kg/kg.
+
+     real                              :: tmp, rmix, virtual
+
+     virtual=tmp*(0.622+rmix)/(0.622*(1.+rmix))
+
+ END FUNCTION virtual
+
+END MODULE module_ForInterpolate
Index: unk/tools/module_generic.F90
===================================================================
--- /trunk/tools/module_generic.F90	(revision 1605)
+++ 	(revision )
@@ -1,95 +1,0 @@
-MODULE module_generic
-! Module with generic functions
-! Index1DArrayR: Function to provide the first index of a given value inside a 1D real array
-! Index2DArrayR: Function to provide the first index of a given value inside a 2D real array
-! ErrWarnMsg: Function to print error/warning message
-
-  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
-  ! Fill value at 64 bits
-  REAL(r_k)                                              :: fillval64 = 1.e20
-  CHARACTER(len=50), PARAMETER                           :: errormsg = 'ERROR -- error -- ERROR -- error'
-  CHARACTER(len=50), PARAMETER                           :: warnmsg = 'WARNING -- warning -- WARNING -- warning'
-
-  CONTAINS
-
-  INTEGER FUNCTION Index1DArrayR(array1D, d1, val)
-! Function to provide the first index of a given value inside a 1D real array
-
-    IMPLICIT NONE
-
-    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
-    INTEGER, INTENT(in)                                  :: d1
-    REAL(r_k), INTENT(in)                                :: val
-    REAL(r_k), DIMENSION(d1), INTENT(in)                 :: array1D
-
-! Local
-    INTEGER                                              :: i
-    CHARACTER(LEN=50)                                    :: fname
-
-    fname = 'Index1DArrayR'
-
-    Index1DArrayR = -1
-
-    DO i=1,d1
-      IF (array1d(i) == val) THEN
-        Index1DArrayR = i
-        EXIT
-      END IF
-    END DO
-
-  END FUNCTION Index1DArrayR
-
-  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
Index: /trunk/tools/module_generic.f90
===================================================================
--- /trunk/tools/module_generic.f90	(revision 1606)
+++ /trunk/tools/module_generic.f90	(revision 1606)
@@ -0,0 +1,95 @@
+MODULE module_generic
+! Module with generic functions
+! Index1DArrayR: Function to provide the first index of a given value inside a 1D real array
+! Index2DArrayR: Function to provide the first index of a given value inside a 2D real array
+! ErrWarnMsg: Function to print error/warning message
+
+  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
+  ! Fill value at 64 bits
+  REAL(r_k)                                              :: fillval64 = 1.e20
+  CHARACTER(len=50), PARAMETER                           :: errormsg = 'ERROR -- error -- ERROR -- error'
+  CHARACTER(len=50), PARAMETER                           :: warnmsg = 'WARNING -- warning -- WARNING -- warning'
+
+  CONTAINS
+
+  INTEGER FUNCTION Index1DArrayR(array1D, d1, val)
+! Function to provide the first index of a given value inside a 1D real array
+
+    IMPLICIT NONE
+
+    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
+    INTEGER, INTENT(in)                                  :: d1
+    REAL(r_k), INTENT(in)                                :: val
+    REAL(r_k), DIMENSION(d1), INTENT(in)                 :: array1D
+
+! Local
+    INTEGER                                              :: i
+    CHARACTER(LEN=50)                                    :: fname
+
+    fname = 'Index1DArrayR'
+
+    Index1DArrayR = -1
+
+    DO i=1,d1
+      IF (array1d(i) == val) THEN
+        Index1DArrayR = i
+        EXIT
+      END IF
+    END DO
+
+  END FUNCTION Index1DArrayR
+
+  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
