diff --git a/tools/NESTING/src/agrif_create_bathy.f90 b/tools/NESTING/src/agrif_create_bathy.f90
index fab9be5b7974675f2b55c05467c04d1966c440d5..cd92999a5d16fc1aa2e038f164df0efe9b709228 100644
--- a/tools/NESTING/src/agrif_create_bathy.f90
+++ b/tools/NESTING/src/agrif_create_bathy.f90
@@ -123,18 +123,20 @@ PROGRAM create_bathy
 
   ! check grids: if identical then do not interpolate
   identical_grids = .FALSE.
-  
-  IF(  SIZE(G0%nav_lat,1) == SIZE(G1%nav_lat,1) .AND. SIZE(G0%nav_lat,2) == SIZE(G1%nav_lat,2) .AND. &
-     & SIZE(G0%nav_lon,1) == SIZE(G1%nav_lon,1) .AND. SIZE(G0%nav_lon,2) == SIZE(G1%nav_lon,2) ) THEN
-     IF(  MAXVAL( ABS(G0%nav_lat(:,:)- G1%nav_lat(:,:)) ) < 0.0001 .AND. &
-        & MAXVAL( ABS(G0%nav_lon(:,:)- G1%nav_lon(:,:)) ) < 0.0001 ) THEN
-        WRITE(*,*) ''
-        WRITE(*,*) 'same grid between parent and child domains => NO INTERPOLATION'
-        WRITE(*,*) ''
-        G1%bathy_meter = G0%bathy_meter
-        identical_grids = .TRUE.
-     ENDIF
-  ENDIF
+ 
+!clem: comment these lines if one wants to create a zoom 1:1 and to change the bathymetry in the zoom
+!
+!  IF(  SIZE(G0%nav_lat,1) == SIZE(G1%nav_lat,1) .AND. SIZE(G0%nav_lat,2) == SIZE(G1%nav_lat,2) .AND. &
+!     & SIZE(G0%nav_lon,1) == SIZE(G1%nav_lon,1) .AND. SIZE(G0%nav_lon,2) == SIZE(G1%nav_lon,2) ) THEN
+!     IF(  MAXVAL( ABS(G0%nav_lat(:,:)- G1%nav_lat(:,:)) ) < 0.0001 .AND. &
+!        & MAXVAL( ABS(G0%nav_lon(:,:)- G1%nav_lon(:,:)) ) < 0.0001 ) THEN
+!        WRITE(*,*) ''
+!        WRITE(*,*) 'same grid between parent and child domains => NO INTERPOLATION'
+!        WRITE(*,*) ''
+!        G1%bathy_meter = G0%bathy_meter
+!        identical_grids = .TRUE.
+!     ENDIF
+!  ENDIF
   
   IF( .NOT.new_topo )   type_bathy_interp = 2   ! only one which works
   !
@@ -270,8 +272,8 @@ PROGRAM create_bathy
         ALLOCATE(masksrc(SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2)))
         ALLOCATE(bathy_test(nxfin,nyfin))
         !
-        WHERE(G0%bathy_meter.LE.0)   ;   masksrc = .FALSE.   ;
-        ELSEWHERE                    ;   masksrc = .TRUE.    ;
+        WHERE(G0%bathy_meter.LE.0.)   ;   masksrc = .FALSE.   ;
+        ELSEWHERE                     ;   masksrc = .TRUE.    ;
         END WHERE
         !            
         ! compute remapping matrix thanks to SCRIP package            
diff --git a/tools/NESTING/src/agrif_create_restart.f90 b/tools/NESTING/src/agrif_create_restart.f90
index 9197147a7c52458a43cb1c06ea7cd92205e35c88..7e53f3db702ee2184bcdad35cb3bdceff9b58a91 100644
--- a/tools/NESTING/src/agrif_create_restart.f90
+++ b/tools/NESTING/src/agrif_create_restart.f90
@@ -99,7 +99,7 @@ PROGRAM create_rstrt
   ! create this file
   !
   CALL set_child_name(restart_file,Child_file)
-  status = nf90_create(Child_file,NF90_WRITE,ncid)
+  status = nf90_create(Child_file,NF90_NETCDF4,ncid)
   status = nf90_close(ncid)
   WRITE(*,*) 'Child grid restart file name = ',TRIM(Child_file)      
   WRITE(*,*) ''
diff --git a/tools/NESTING/src/agrif_interpolation.f90 b/tools/NESTING/src/agrif_interpolation.f90
index e28a0ef55a9730f1be0afa92de716918ab6d1fa9..bb507ffc9d01e0bae76a80f4997b5dbb4ac9c168 100644
--- a/tools/NESTING/src/agrif_interpolation.f90
+++ b/tools/NESTING/src/agrif_interpolation.f90
@@ -210,21 +210,21 @@ CONTAINS
     SELECT CASE(typevar)
     CASE('T')
        IF(MOD(irafx,2)==1) THEN ! odd
-          zx = 1 ; zy = 1 ; jdecx = FLOOR(irafx/2.) ; jdecy = FLOOR(irafy/2.)
+          zx = 1 ; zy = 1 ; jdecx = FLOOR(REAL(irafx)/2.) ; jdecy = FLOOR(REAL(irafy)/2.)
        ELSE                     ! even
-          zx = 2 ; zy = 2 ; jdecx = FLOOR(irafx/2.) ; jdecy = FLOOR(irafy/2.)
+          zx = 2 ; zy = 2 ; jdecx = FLOOR(REAL(irafx)/2.) ; jdecy = FLOOR(REAL(irafy)/2.)
        ENDIF
     CASE('U')
        IF(MOD(irafx,2)==1) THEN ! odd
-          zx = 1 ; zy = 1 ; jdecx = irafx - 1 ; jdecy = FLOOR(irafy/2.)
+          zx = 1 ; zy = 1 ; jdecx = irafx - 1 ; jdecy = FLOOR(REAL(irafy)/2.)
        ELSE                     ! even
-          zx = 1 ; zy = 2 ; jdecx = irafx - 1 ; jdecy = FLOOR(irafy/2.)
+          zx = 1 ; zy = 2 ; jdecx = irafx - 1 ; jdecy = FLOOR(REAL(irafy)/2.)
        ENDIF
     CASE('V')
        IF(MOD(irafx,2)==1) THEN ! odd
-          zx = 1 ; zy = 1 ; jdecx = FLOOR(irafx/2.) ; jdecy = irafy - 1
+          zx = 1 ; zy = 1 ; jdecx = FLOOR(REAL(irafx)/2.) ; jdecy = irafy - 1
        ELSE                     ! even
-          zx = 2 ; zy = 1 ; jdecx = FLOOR(irafx/2.) ; jdecy = irafy - 1
+          zx = 2 ; zy = 1 ; jdecx = FLOOR(REAL(irafx)/2.) ; jdecy = irafy - 1
        ENDIF
     CASE('F')
        IF(MOD(irafx,2)==1) THEN ! odd
@@ -238,21 +238,24 @@ CONTAINS
     DO jj = 1, nyf
 
        jjf = jj - jdecy
-       jjc = j_min + FLOOR((jjf-1.) / irafy) 
-       
+       jjc = j_min + FLOOR(REAL(jjf-1) / REAL(irafy)) 
+
        DO ji = 1, nxf
           
           jif = ji - jdecx 
-          jic = i_min + FLOOR((jif-1.) / irafx) 
-          
-          Bx = MOD( zx*jif-1, zx*irafx ) / REAL(zx*irafx)
-          By = MOD( zy*jjf-1, zy*irafy ) / REAL(zy*irafy)
+          jic = i_min + FLOOR(REAL(jif-1) / REAL(irafx)) 
+
+          Bx = REAL(MODULO( zx*jif-1, zx*irafx )) / REAL(zx*irafx)
+          By = REAL(MODULO( zy*jjf-1, zy*irafy )) / REAL(zy*irafy)
           Ax = 1. - Bx
           Ay = 1. - By
 
           jic1 = MIN( nxc, jic+1 ) ! avoid out of bounds for tabin below
           jjc1 = MIN( nyc, jjc+1 ) !             --
 
+          jic = MAX( 1, jic ) ! avoid out of bounds
+          jjc = MAX( 1, jjc ) ! avoid out of bounds
+
           tabout(ji,jj) = ( Bx * tabin(jic1,jjc ) + Ax * tabin(jic,jjc ) ) * Ay + &
              &            ( Bx * tabin(jic1,jjc1) + Ax * tabin(jic,jjc1) ) * By
        END DO
@@ -503,7 +506,7 @@ CONTAINS
     REAL*8, POINTER, DIMENSION(:,:) :: lonParent,latParent => NULL()
     REAL*8, POINTER, DIMENSION(:,:) :: lonChild,latChild,latlon_temp => NULL()
     REAL*8, POINTER, DIMENSION(:,:,:,:) :: tabinterp4d,tabvar1,tabvar2,tabvar3 => NULL()
-    REAL*8, POINTER, DIMENSION(:,:,:) :: tabinterp3d,tabvar3d => NULL()
+    REAL*8, POINTER, DIMENSION(:,:,:) :: tabinterp3d,tabvar3d,tabvar01,tabvar02,tabvar03 => NULL()
     REAL*8, POINTER, DIMENSION(:) :: timedepth_temp,depth => NULL()
     REAL*8,DIMENSION(:,:),POINTER :: matrix => NULL()
     INTEGER,DIMENSION(:),POINTER :: src_add,dst_add  => NULL()
@@ -515,11 +518,17 @@ CONTAINS
     LOGICAL,DIMENSION(:,:),POINTER :: masksrc => NULL()
     LOGICAL :: Interpolation,conservation,Pacifique,Extrapolation,land_level
     !      
-    INTEGER :: deptht,time,i,status,ncid,t,ii,j,nb,numlon,numlat
+    INTEGER :: deptht,time,i,status,ncid,t,ii,j,nb,numlon,numlat,unlimdimid
 
     !        
     TYPE(Coordinates) :: G0,G1      
     !
+    status = nf90_open(TRIM(filename),NF90_NOWRITE,ncid)
+    IF (status/=nf90_noerr) THEN
+       WRITE(*,*)"unable to open netcdf file : ",TRIM(filename)
+       STOP
+    ENDIF
+    !
     !*****************
     !If coarse grid is masked possibility to activate an extrapolation process
     !
@@ -541,7 +550,6 @@ CONTAINS
     !
     CALL Read_Ncdf_dim('x',filename,numlon)
     CALL Read_Ncdf_dim('y',filename,numlat)
-    CALL Read_Ncdf_dim('time_counter',filename,time)
     IF ( Dims_Existence( 'deptht' , filename ) ) THEN
        CALL Read_Ncdf_dim('deptht',filename,deptht)
     ELSE IF ( Dims_Existence( 'depthu' , filename ) ) THEN
@@ -555,7 +563,28 @@ CONTAINS
     ELSE
        deptht = N
     ENDIF
-
+    IF    ( Dims_Existence( 'time_counter' , filename ) ) THEN
+       CALL Read_Ncdf_dim('time_counter',filename,time)
+       ! check that time_counter is the unlimited dim
+       status = nf90_inquire(ncid, unlimiteddimid = unlimdimid)
+       IF ( unlimdimid == -1 ) THEN
+          WRITE(*,*)"time_counter should be the unlimited dimension in this file : ",filename
+          WRITE(*,*)" use nco command: ncks -O --mk_rec_dmn time_counter ",filename," ",filename
+          STOP
+       ENDIF
+    ELSEIF( Dims_Existence( 'time' , filename ) ) THEN
+       CALL Read_Ncdf_dim('time',filename,time)
+       ! check that time is the unlimited dim
+       status = nf90_inquire(ncid, unlimiteddimid = unlimdimid)
+       IF ( unlimdimid == -1 ) THEN
+          WRITE(*,*)"time should be the unlimited dimension in this file : ",filename
+          WRITE(*,*)" use nco command: ncks -O --mk_rec_dmn time ",filename," ",filename
+          STOP
+       ENDIF
+    ELSE
+       time = 0
+    ENDIF
+    
     !
     ! retrieve netcdf variable name
     !
@@ -579,7 +608,7 @@ CONTAINS
     status = Read_Coordinates(Childcoordinates,G1,Pacifique)
     !
     ! check consistency of informations read in namelist 
-    !      
+    !
     IF( imax > SIZE(G0%glamt,1) .OR. jmax > SIZE(G0%glamt,2) .OR. &
          imax <= imin .OR. jmax <= jmin ) THEN                    
        WRITE(*,*) 'ERROR ***** bad child grid definition ...'
@@ -597,7 +626,6 @@ CONTAINS
        STOP
     ENDIF
     !      
-
     !
     ! Initialization of T-mask thanks to bathymetry 
     !
@@ -610,7 +638,6 @@ CONTAINS
        CALL Init_mask(parent_bathy_level,G0,1,1)
        !
     ENDIF
-
     !      
     ! select coordinates to use according to variable position
     !
@@ -619,23 +646,67 @@ CONTAINS
 
     SELECT CASE(posvar)
     CASE('T')
-       lonParent = G0%glamt
-       latParent = G0%gphit
-       lonChild  = G1%glamt
-       latChild  = G1%gphit
-       mask      = G1%tmask
+       IF ( Vars_Existence( 'glamt' , parent_coordinate_file ) ) THEN
+          lonParent = G0%glamt
+          latParent = G0%gphit
+       ELSEIF( Vars_Existence( 'nav_lon' , parent_coordinate_file ) ) THEN
+          lonParent = G0%nav_lon
+          latParent = G0%nav_lat
+       ENDIF
+       IF ( Vars_Existence( 'glamt' , Childcoordinates ) ) THEN
+          lonChild  = G1%glamt
+          latChild  = G1%gphit
+       ELSEIF( Vars_Existence( 'nav_lon' , Childcoordinates ) ) THEN
+          lonChild  = G1%nav_lon
+          latChild  = G1%nav_lat
+       ENDIF
+       !!IF ( Vars_Existence( 'tmask' , Childcoordinates ) ) THEN
+          mask = G1%tmask
+       !!ELSEIF( Vars_Existence( 'top_level' , Childcoordinates ) ) THEN
+       !!   mask = MAX( 1., G1%top_level )
+       !!ELSE
+       !!   mask = 1.
+       !!ENDIF
     CASE('U')
-       lonParent = G0%glamu
-       latParent = G0%gphiu
-       lonChild  = G1%glamu
-       latChild  = G1%gphiu
-       mask      = G1%umask
+       IF ( Vars_Existence( 'glamu' , parent_coordinate_file ) ) THEN
+          lonParent = G0%glamu
+          latParent = G0%gphiu
+       ELSE
+          WRITE(*,*) 'ERROR ***** missing glamu/gphiu coordinates in ', TRIM(parent_coordinate_file)
+          STOP
+       ENDIF
+       IF ( Vars_Existence( 'glamu' , Childcoordinates ) ) THEN
+          lonChild  = G1%glamu
+          latChild  = G1%gphiu
+       ELSE
+          WRITE(*,*) 'ERROR ***** missing glamu/gphiu coordinates in ', TRIM(Childcoordinates)
+          STOP
+       ENDIF
+       !!IF ( Vars_Existence( 'umask' , Childcoordinates ) ) THEN
+          mask = G1%umask
+       !!ELSE
+       !!   mask = 1.
+       !!ENDIF
     CASE('V')
-       lonParent = G0%glamv
-       latParent = G0%gphiv
-       lonChild  = G1%glamv
-       latChild  = G1%gphiv
-       mask      = G1%vmask
+       IF ( Vars_Existence( 'glamv' , parent_coordinate_file ) ) THEN
+          lonParent = G0%glamv
+          latParent = G0%gphiv
+       ELSE
+          WRITE(*,*) 'ERROR ***** missing glamv/gphiv coordinates in ', TRIM(parent_coordinate_file)
+          STOP
+       ENDIF
+       IF ( Vars_Existence( 'glamv' , Childcoordinates ) ) THEN
+          lonChild  = G1%glamv
+          latChild  = G1%gphiv
+       ELSE
+          WRITE(*,*) 'ERROR ***** missing glamv/gphiv coordinates in ', TRIM(Childcoordinates)
+          STOP
+       ENDIF
+       !!IF ( Vars_Existence( 'vmask' , Childcoordinates ) ) THEN
+          mask = G1%vmask
+       !!ELSE
+       !!   mask = 1.
+       !!ENDIF
     END SELECT
 
     DEALLOCATE(G0%glamu,G0%glamv,G0%gphiu,G0%gphiv)
@@ -665,8 +736,9 @@ CONTAINS
     IF ( Dims_Existence( 'depthu' , filename ) ) CALL Write_Ncdf_dim('depthu',Child_file,deptht)
     IF ( Dims_Existence( 'depthv' , filename ) ) CALL Write_Ncdf_dim('depthv',Child_file,deptht)
     IF ( Dims_Existence( 'depthw' , filename ) ) CALL Write_Ncdf_dim('depthw',Child_file,deptht)
-    IF ( Dims_Existence( 'z' , filename ) ) CALL Write_Ncdf_dim('z',Child_file,deptht)
-    CALL Write_Ncdf_dim('time_counter',Child_file,0)
+    IF ( Dims_Existence( 'z'      , filename ) ) CALL Write_Ncdf_dim('z'     ,Child_file,deptht)
+    IF ( Dims_Existence( 'time_counter' , filename ) ) CALL Write_Ncdf_dim('time_counter',Child_file,time)
+    IF ( Dims_Existence( 'time'         , filename ) ) CALL Write_Ncdf_dim('time'        ,Child_file,time)
 
     IF( deptht .NE. 1 .AND. deptht .NE. N ) THEN
        WRITE(*,*) '***'
@@ -733,6 +805,17 @@ CONTAINS
           varname = TRIM(Ncdf_varname(i))
           Interpolation = .FALSE.
           !
+          !copy time from input file to output file
+       CASE('time')
+          ALLOCATE(timedepth_temp(time))
+          CALL Read_Ncdf_var('time',filename,timedepth_temp) 
+          CALL Write_Ncdf_var('time','time',  &
+               Child_file,timedepth_temp,'float')
+          CALL Copy_Ncdf_att('time',TRIM(filename),Child_file)
+          DEALLOCATE(timedepth_temp)
+          varname = TRIM(Ncdf_varname(i))
+          Interpolation = .FALSE.
+          !
           !copy deptht from input file to output file
        CASE('deptht')
           ALLOCATE(depth(deptht))
@@ -806,84 +889,305 @@ CONTAINS
           !	         
        END SELECT
 
-       ! //////////////// INTERPOLATION FOR 3D VARIABLES /////////////////////////////////////
+       ! //////////////// INTERPOLATION FOR 2D VARIABLES /////////////////////////////////////
        ! 
-       IF( Interpolation .AND. Get_NbDims(TRIM(varname),TRIM(filename)) == 3 ) THEN
+       IF( Interpolation .AND. Get_NbDims(TRIM(varname),TRIM(filename)) == 2 ) THEN
           !      
           ALLOCATE(detected_pts(numlon,numlat,N))
           ALLOCATE(masksrc(numlon,numlat))                                           
           !
-          ! ******************************LOOP ON TIME*******************************************
-          !loop on time
-          DO t = 1,time
-             !                   
-             IF(extrapolation) THEN
-                WRITE(*,*) 'interpolation/extrapolation ',TRIM(varname),' for time t = ',t
-             ELSE
-                WRITE(*,*) 'interpolation ',TRIM(varname),' for time t = ',t
-             ENDIF
-             !                            
-             ALLOCATE(tabvar3d(numlon,numlat,1))      
-             ALLOCATE(tabinterp3d(nxfin,nyfin,1))
-             !                    
-             CALL Read_Ncdf_var(varname,filename,tabvar3d,t)                                              
+          IF(extrapolation) THEN
+             WRITE(*,*) 'interpolation/extrapolation ',TRIM(varname)
+          ELSE
+             WRITE(*,*) 'interpolation ',TRIM(varname)
+          ENDIF
+          !                            
+          ALLOCATE(tabvar3d(numlon,numlat,1))      
+          ALLOCATE(tabinterp3d(nxfin,nyfin,1))
+          !                    
+          CALL Read_Ncdf_var(varname,filename,tabvar3d)                                              
+          !
+          ! search points where extrapolation is required
+          ! 
+          IF(Extrapolation) THEN
+             WHERE( tabvar3d .GE. 1.e+20 ) tabvar3d = 0.
+             CALL extrap_detect(G0,G1,detected_pts(:,:,1),1)
+             CALL correct_field_2d(detected_pts(:,:,1),tabvar3d,G0,masksrc,'posvar')
+          ELSE
+             masksrc = .TRUE.
+          ENDIF
+
+          SELECT CASE(TRIM(interp_type))
+          CASE('bilinear')
+             CALL get_remap_matrix(latParent,latChild,   &
+                lonParent,lonChild,                    &
+                masksrc,matrix,src_add,dst_add)
+
+          CASE('bicubic')
+             CALL get_remap_bicub(latParent,latChild,   &
+                lonParent,lonChild,   &
+                masksrc,matrix,src_add,dst_add)
+
+          END SELECT
+          !                                                        
+          !      
+          SELECT CASE(TRIM(interp_type))
+          CASE('bilinear')                                                       
+             CALL make_remap(tabvar3d(:,:,1),tabinterp3d(:,:,1),nxfin,nyfin, &
+                matrix,src_add,dst_add)     
+          CASE('bicubic')                                   
+             CALL make_bicubic_remap(tabvar3d(:,:,1),masksrc,tabinterp3d(:,:,1),nxfin,nyfin, &
+                matrix,src_add,dst_add)                        
+          END SELECT
+          !                     
+          IF( conservation ) CALL Correctforconservation(tabvar3d(:,:,1),tabinterp3d(:,:,1), &
+             G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin)  
+          !      
+          IF(Extrapolation) tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,1)     
+!!          tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,1)
+          !                     
+          dimnames(1)='x'
+          dimnames(2)='y'
+          !                                     
+          CALL Write_Ncdf_var(TRIM(varname),dimnames(1:2),&
+             Child_file,tabinterp3d(:,:,1),'float')
+          !      
+          DEALLOCATE(tabinterp3d)
+          DEALLOCATE(tabvar3d)                     
+          !                   
+          DEALLOCATE(detected_pts)
+          IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,dst_add,src_add)
+          DEALLOCATE( masksrc)
+          
+          CALL Copy_Ncdf_att(TRIM(varname),TRIM(filename),Child_file)                        
+          !
+          
+       ! //////////////// INTERPOLATION FOR 3D VARIABLES /////////////////////////////////////
+       ! 
+       ELSEIF( Interpolation .AND. Get_NbDims(TRIM(varname),TRIM(filename)) == 3 ) THEN
+          !      
+          IF( DimUnlimited_Var(TRIM(varname),TRIM(filename)) ) THEN
+
+             ALLOCATE(detected_pts(numlon,numlat,N))
+             ALLOCATE(masksrc(numlon,numlat))                                           
              !
-             ! search points where extrapolation is required
-             ! 
-             IF(Extrapolation) THEN
-                 WHERE( tabvar3d .GE. 1.e+20 ) tabvar3d = 0.
-                IF (t .EQ. 1. ) CALL extrap_detect(G0,G1,detected_pts(:,:,1),1)
-                CALL correct_field_2d(detected_pts(:,:,1),tabvar3d,G0,masksrc,'posvar')
-             ELSE
-                masksrc = .TRUE.
-             ENDIF
+             ! ******************************LOOP ON TIME*******************************************
+             !loop on time
+             DO t = 1,time
+                !                   
+                IF(extrapolation) THEN
+                   WRITE(*,*) 'interpolation/extrapolation ',TRIM(varname),' for time t = ',t
+                ELSE
+                   WRITE(*,*) 'interpolation ',TRIM(varname),' for time t = ',t
+                ENDIF
+                !                            
+                ALLOCATE(tabvar3d(numlon,numlat,1))      
+                ALLOCATE(tabinterp3d(nxfin,nyfin,1))
+                !                    
+                CALL Read_Ncdf_var(varname,filename,tabvar3d,t)                                              
+                !
+                ! search points where extrapolation is required
+                ! 
+                IF(Extrapolation) THEN
+                   WHERE( tabvar3d .GE. 1.e+20 ) tabvar3d = 0.
+                   IF (t .EQ. 1. ) CALL extrap_detect(G0,G1,detected_pts(:,:,1),1)
+                   CALL correct_field_2d(detected_pts(:,:,1),tabvar3d,G0,masksrc,'posvar')
+                ELSE
+                   masksrc = .TRUE.
+                ENDIF
+
+                IF (t.EQ.1 ) THEN 
 
-             IF (t.EQ.1 ) THEN 
+                   SELECT CASE(TRIM(interp_type))
+                   CASE('bilinear')
+                      CALL get_remap_matrix(latParent,latChild,   &
+                         lonParent,lonChild,                    &
+                         masksrc,matrix,src_add,dst_add)
+                   CASE('bicubic')
+                      CALL get_remap_bicub(latParent,latChild,   &
+                         lonParent,lonChild,   &
+                         masksrc,matrix,src_add,dst_add)
 
+                   END SELECT
+                   !                                                        
+                ENDIF
+                !      
                 SELECT CASE(TRIM(interp_type))
-                CASE('bilinear')
-                   CALL get_remap_matrix(latParent,latChild,   &
-                        lonParent,lonChild,                    &
-                        masksrc,matrix,src_add,dst_add)
+                CASE('bilinear')                                                       
+                   CALL make_remap(tabvar3d(:,:,1),tabinterp3d(:,:,1),nxfin,nyfin, &
+                      matrix,src_add,dst_add)     
+                CASE('bicubic')                                   
+                   CALL make_bicubic_remap(tabvar3d(:,:,1),masksrc,tabinterp3d(:,:,1),nxfin,nyfin, &
+                      matrix,src_add,dst_add)                        
+                END SELECT
+                !                     
+                IF( conservation ) CALL Correctforconservation(tabvar3d(:,:,1),tabinterp3d(:,:,1), &
+                   G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin)  
+                !      
+                IF(Extrapolation) tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,1)     
+!!                tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,1)
+                !                     
+                dimnames(1)='x'
+                dimnames(2)='y'
+                dimnames(3)='time_counter'
+                !                                     
+                CALL Write_Ncdf_var(TRIM(varname),dimnames(1:3),&
+                   Child_file,tabinterp3d,t,'float')
+                !      
+                DEALLOCATE(tabinterp3d)
+                DEALLOCATE(tabvar3d)                     
+                !end loop on time
+             END DO
+             !                   
+             DEALLOCATE(detected_pts)
+             IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,dst_add,src_add)
+             DEALLOCATE( masksrc)
 
-                CASE('bicubic')
-                   CALL get_remap_bicub(latParent,latChild,   &
-                        lonParent,lonChild,   &
-                        masksrc,matrix,src_add,dst_add)
+          ELSE
 
-                END SELECT
-                !                                                        
-             ENDIF
-             !      
-             SELECT CASE(TRIM(interp_type))
-             CASE('bilinear')                                                       
-                CALL make_remap(tabvar3d(:,:,1),tabinterp3d(:,:,1),nxfin,nyfin, &
-                     matrix,src_add,dst_add)     
-             CASE('bicubic')                                   
-                CALL make_bicubic_remap(tabvar3d(:,:,1),masksrc,tabinterp3d(:,:,1),nxfin,nyfin, &
-                     matrix,src_add,dst_add)                        
-             END SELECT
-             !                     
-             IF( conservation ) CALL Correctforconservation(tabvar3d(:,:,1),tabinterp3d(:,:,1), &
-                  G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin)  
-             !      
-             IF(Extrapolation) tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,1)     
-             !                     
              dimnames(1)='x'
              dimnames(2)='y'
-             dimnames(3)='time_counter'
-             !                                     
-             CALL Write_Ncdf_var(TRIM(varname),dimnames(1:3),&
-                  Child_file,tabinterp3d,t,'float')
-             !      
-             DEALLOCATE(tabinterp3d)
-             DEALLOCATE(tabvar3d)                     
-             !end loop on time
-          END DO
-          !                   
-          DEALLOCATE(detected_pts)
-          IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,dst_add,src_add)
-          DEALLOCATE( masksrc)
+             IF ( Dims_Existence( 'deptht' , filename ) ) dimnames(3)='deptht'  
+             IF ( Dims_Existence( 'depthu' , filename ) ) dimnames(3)='depthu'
+             IF ( Dims_Existence( 'depthv' , filename ) ) dimnames(3)='depthv'
+             IF ( Dims_Existence( 'depthw' , filename ) ) dimnames(3)='depthw'
+             IF ( Dims_Existence( 'z' , filename ) ) dimnames(3)='z'
+             !
+             ! loop on vertical levels
+             ! 
+             DO nb = 1,deptht
+                !        
+                ALLOCATE(masksrc(numlon,numlat))
+                ALLOCATE(detected_pts(numlon,numlat,N))          
+                !
+                ! point detection et level n                                     
+                !
+                land_level = .FALSE.
+                IF( Extrapolation ) THEN
+                   IF(MAXVAL(mask(:,:,nb))==0.) land_level = .TRUE.
+                ENDIF
+
+
+                IF ( land_level ) THEN
+                   !
+                   WRITE(*,*) 'only land points on level ',nb
+                   ALLOCATE(tabinterp3d(nxfin,nyfin,1)) 
+                   tabinterp3d = 0.e0
+                   !                              
+                   CALL Write_Ncdf_var(TRIM(varname),dimnames, &
+                      Child_file,tabinterp3d,nb,'float')
+                   DEALLOCATE(tabinterp3d)
+                   !
+                ELSE
+                   !                       
+                   ALLOCATE(tabvar01(numlon,numlat,1))    ! level k
+                   IF(Extrapolation) ALLOCATE(tabvar02(numlon,numlat,1))    ! level k-1
+                   IF(Extrapolation) ALLOCATE(tabvar03(numlon,numlat,1))    ! level k-2                       
+                   ALLOCATE(tabinterp3d(nxfin,nyfin,1)) 
+                   !
+                   IF(Extrapolation) THEN                                  
+                      IF(nb==1) THEN
+                         CALL Read_Ncdf_var(varname,filename,tabvar01,nb)    
+                         WHERE( tabvar01 .GE. 1.e+20 ) tabvar01 = 0.
+                      ELSE IF (nb==2) THEN
+                         CALL Read_Ncdf_var(varname,filename,tabvar02,nb-1)
+                         CALL Read_Ncdf_var(varname,filename,tabvar01,nb)            
+                         WHERE( tabvar01 .GE. 1.e+20 ) tabvar01 = 0.
+                         WHERE( tabvar02 .GE. 1.e+20 ) tabvar02 = 0.
+                      ELSE 
+                         CALL Read_Ncdf_var(varname,filename,tabvar03,nb-2)
+                         CALL Read_Ncdf_var(varname,filename,tabvar02,nb-1)
+                         CALL Read_Ncdf_var(varname,filename,tabvar01,nb)
+                         WHERE( tabvar01 .GE. 1.e+20 ) tabvar01 = 0.
+                         WHERE( tabvar02 .GE. 1.e+20 ) tabvar02 = 0.
+                         WHERE( tabvar03 .GE. 1.e+20 ) tabvar03 = 0.
+                      ENDIF
+                      !                              
+                      CALL extrap_detect(G0,G1,detected_pts(:,:,nb),nb)
+
+                      ALLOCATE(tabvar1(numlon,numlat,1,1))    ! level k
+                      ALLOCATE(tabvar2(numlon,numlat,1,1))    ! level k-1
+                      ALLOCATE(tabvar3(numlon,numlat,1,1))    ! level k-2
+                      tabvar1(:,:,1,1) = tabvar01(:,:,1)
+                      tabvar2(:,:,1,1) = tabvar02(:,:,1)
+                      tabvar3(:,:,1,1) = tabvar03(:,:,1)
+                      CALL correct_field(detected_pts(:,:,nb),tabvar1,tabvar2,&
+                         tabvar3,G0,depth,masksrc,nb,'posvar')
+                      tabvar01(:,:,1) = tabvar1(:,:,1,1)
+                      tabvar02(:,:,1) = tabvar2(:,:,1,1)
+                      tabvar03(:,:,1) = tabvar3(:,:,1,1)                      
+                      DEALLOCATE(tabvar1,tabvar2,tabvar3)
+                      DEALLOCATE(tabvar02,tabvar03)
+
+                   ELSE 
+                      CALL Read_Ncdf_var(varname,filename,tabvar01,nb)
+                      IF(MAXVAL(tabvar01(:,:,1))==0.) land_level = .TRUE.
+                      masksrc = .TRUE. 
+                   ENDIF
+
+                   IF( Extrapolation ) THEN
+                      WRITE(*,*) 'interpolation/extrapolation ',TRIM(varname),' for vertical level = ',nb   
+                   ELSE
+                      WRITE(*,*) 'interpolation ',TRIM(varname),' for vertical level = ',nb
+                   ENDIF
+                   !
+                   SELECT CASE(TRIM(interp_type))
+                   CASE('bilinear')
+                      CALL get_remap_matrix(latParent,latChild,   &
+                         lonParent,lonChild,   &
+                         masksrc,matrix,src_add,dst_add)
+
+                   CASE('bicubic')
+                      CALL get_remap_bicub(latParent,latChild,   &
+                         lonParent,lonChild,   &
+                         masksrc,matrix,src_add,dst_add)
+                      !                             
+                   END SELECT
+                   !                                                        
+                   !      
+                   SELECT CASE(TRIM(interp_type))
+                      !                             
+                   CASE('bilinear')                                                       
+                      CALL make_remap(tabvar01(:,:,1),tabinterp3d(:,:,1),nxfin,nyfin, &
+                         matrix,src_add,dst_add)     
+                   CASE('bicubic')                                   
+                      CALL make_bicubic_remap(tabvar01(:,:,1),masksrc,tabinterp3d(:,:,1),nxfin,nyfin, &
+                         matrix,src_add,dst_add)                        
+                   END SELECT
+                   !                     
+                   IF( conservation ) CALL Correctforconservation(tabvar01(:,:,1),tabinterp3d(:,:,1), &
+                      G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin)  
+
+                   !
+                   ALLOCATE(tabinterp4d(nxfin,nyfin,1,1))
+                   tabinterp4d(:,:,1,1) = tabinterp3d(:,:,1)
+                   IF(Extrapolation) CALL check_extrap(G1,tabinterp4d,nb)
+                   tabinterp3d(:,:,1) = tabinterp4d(:,:,1,1)
+                   DEALLOCATE(tabinterp4d)
+                   !     
+                   IF(Extrapolation) tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,nb)     
+!!                   tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,nb)
+                   !                     
+                   CALL Write_Ncdf_var(TRIM(varname),dimnames, &
+                      Child_file,tabinterp3d,nb,'float')
+                   !      
+                   DEALLOCATE(tabinterp3d)
+                   DEALLOCATE(tabvar01)                     
+                   !
+                   !
+                ENDIF
+
+                !
+                IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,dst_add,src_add)
+                DEALLOCATE( masksrc )
+                DEALLOCATE(detected_pts)       
+                !
+                ! end loop on vertical levels
+                !      
+             END DO
+
+             
+          ENDIF
 
           CALL Copy_Ncdf_att(TRIM(varname),TRIM(filename),Child_file)                        
           !
diff --git a/tools/NESTING/src/agrif_readwrite.f90 b/tools/NESTING/src/agrif_readwrite.f90
index 87e50abdd926ee0a4bdfc79b725645f22b07dfb5..f86aeea5eb1c4bac189fad4f1615d0464279a8db 100644
--- a/tools/NESTING/src/agrif_readwrite.f90
+++ b/tools/NESTING/src/agrif_readwrite.f90
@@ -1128,7 +1128,7 @@ CONTAINS
     USE io_netcdf
     !      
     CHARACTER(*) name
-    INTEGER :: nx,ny,k,i,j,jpiglo,jpjglo
+    INTEGER :: nx,ny,k,i,j,jpiglo,jpjglo,ji,jj
     TYPE(Coordinates) :: Grid
     REAL*8, POINTER, DIMENSION(:,:) ::zwf => NULL()
     !
@@ -1145,7 +1145,6 @@ CONTAINS
     nx = SIZE(Grid%Bathy_level,1)
     ny = SIZE(Grid%Bathy_level,2)
     !
-    !      
     ALLOCATE(Grid%tmask(nx,ny,N), &
          Grid%umask(nx,ny,N), &
          Grid%vmask(nx,ny,N), &
@@ -1154,21 +1153,34 @@ CONTAINS
     DO k = 1,N
        !      
        WHERE(Grid%Bathy_level(:,:) <= k-1 )    
-          Grid%tmask(:,:,k) = 0
+          Grid%tmask(:,:,k) = 0.
        ELSEWHERE
-          Grid%tmask(:,:,k) = 1
+          Grid%tmask(:,:,k) = 1.
        END WHERE
        !
     END DO
     !
-    Grid%umask(1:nx-1,:,:) = Grid%tmask(1:nx-1,:,:)*Grid%tmask(2:nx,:,:)
-    Grid%vmask(:,1:ny-1,:) = Grid%tmask(:,1:ny-1,:)*Grid%tmask(:,2:ny,:)
+    DO jj=1,ny
+       DO ji=1,nx-1
+          Grid%umask(ji,jj,:) = Grid%tmask(ji,jj,:) * Grid%tmask(ji+1,jj,:)
+       ENDDO
+    ENDDO
+    DO jj=1,ny-1
+       DO ji=1,nx
+          Grid%vmask(ji,jj,:) = Grid%tmask(ji,jj,:)*Grid%tmask(ji,jj+1,:)
+       ENDDO
+    ENDDO
     !
     Grid%umask(nx,:,:) = Grid%tmask(nx,:,:)
     Grid%vmask(:,ny,:) = Grid%tmask(:,ny,:)
-    !      
-    Grid%fmask(1:nx-1,1:ny-1,:) = Grid%tmask(1:nx-1,1:ny-1,:)*Grid%tmask(2:nx,1:ny-1,:)* &
-         Grid%tmask(1:nx-1,2:ny,:)*Grid%tmask(2:nx,2:ny,:) 
+
+    !    
+    DO jj=1,ny-1
+       DO ji=1,nx-1  
+          Grid%fmask(ji,jj,:) = Grid%tmask(ji,jj,:)*Grid%tmask(ji+1,jj,:)* &
+                                Grid%tmask(ji,jj+1,:)*Grid%tmask(ji+1,jj+1,:) 
+       ENDDO
+    ENDDO
     !
     Grid%fmask(nx,:,:) = Grid%tmask(nx,:,:)
     Grid%fmask(:,ny,:) = Grid%tmask(:,ny,:)
@@ -1233,6 +1245,8 @@ CONTAINS
     nx = SIZE(Grid%Bathy_level,1)
     ny = SIZE(Grid%Bathy_level,2) 
     !
+    WRITE(*,*) 'nx=',nx,' ny=',ny
+    !
     WRITE(*,*) 'Init masks in T points'    
     !     
     ALLOCATE( Grid%tmask(nx,ny,N) )
diff --git a/tools/NESTING/src/bilinear_interp.f90 b/tools/NESTING/src/bilinear_interp.f90
index 695a700bc009d3990ec8e0a778b5b9eec05ac584..444eca5154004e8a32800ae413d1dc1654136a0b 100644
--- a/tools/NESTING/src/bilinear_interp.f90
+++ b/tools/NESTING/src/bilinear_interp.f90
@@ -205,7 +205,6 @@ CONTAINS
     ALLOCATE(grid1_mask(SIZE(grid1_lat,1)*SIZE(grid1_lat,2)))
     CALL logtab2Dto1D(mask,grid1_mask)
     !      
-    !      Write(*,*) ,'grid1_mask = ',grid1_mask                 
     !
     ! degrees to radian
     !
@@ -529,6 +528,7 @@ CONTAINS
     ! bilinear weights for four corners
     !      
     REAL*8, DIMENSION(4) :: wgts            
+    INTEGER, DIMENSION(1) :: wmax !!clem nearest
     !
     REAL*8 :: &
          plat, plon,       &  ! lat/lon coords of destination point
@@ -653,6 +653,11 @@ CONTAINS
              wgts(3) = iguess*jguess
              wgts(4) = (one-iguess)*jguess
              !	    	        
+!!clem nearest
+!wmax=MAX( 1, MAXLOC( wgts(1:4), MASK = (grid1_mask(src_add(1:4)) == .TRUE.) ) )
+!wgts(1:4)=0.
+!wgts(wmax)=1.
+!!clem
              !
              CALL store_link_bilin(dst_add, src_add, wgts, nmap)
 
@@ -669,7 +674,7 @@ CONTAINS
           !*** search for bilinear failed - use a distance-weighted
           !*** average instead (this is typically near the pole)
           !***
-       ELSE IF (src_add(1) < 0) THEN
+       ELSE IF (src_add(1) <= 0) THEN
 
           src_add = ABS(src_add)
           icount = 0
@@ -693,6 +698,11 @@ CONTAINS
              wgts(3) = src_lats(3)/sum_wgts
              wgts(4) = src_lats(4)/sum_wgts
              !
+!!clem nearest
+!wmax=MAX( 1, MAXLOC( wgts(1:4), MASK = (grid1_mask(src_add(1:4)) == .TRUE.) ) )
+!wgts(1:4)=0.
+!wgts(wmax)=1.
+!!clem
              grid2_frac(dst_add) = one
              CALL store_link_bilin(dst_add, src_add, wgts, nmap)
           ENDIF
diff --git a/tools/NESTING/src/io_netcdf.f90 b/tools/NESTING/src/io_netcdf.f90
index 98a407653f3b6ffefe2f3094a54ed6413908e2a8..b870cf5dc87fbf88daba1f972f83c43441831b3c 100644
--- a/tools/NESTING/src/io_netcdf.f90
+++ b/tools/NESTING/src/io_netcdf.f90
@@ -1603,4 +1603,53 @@ CONTAINS
     !
   END FUNCTION Dims_Existence
   !
+  LOGICAL FUNCTION Vars_Existence( varname , filename )
+    !
+    CHARACTER(*),INTENT(in) :: varname,filename
+    INTEGER :: status,ncid,varid    
+    !      
+    status = nf90_open(TRIM(filename),NF90_NOWRITE,ncid)
+    IF (status/=nf90_noerr) THEN    
+       WRITE(*,*)"unable to open netcdf file : ",TRIM(filename)
+       STOP
+    ENDIF
+    status = nf90_inq_varid(ncid,TRIM(varname),varid)
+    !      
+    IF (status/=nf90_noerr) THEN
+       Vars_Existence = .FALSE.
+    ELSE
+       Vars_Existence = .TRUE.
+    ENDIF
+    !
+    RETURN
+    !
+  END FUNCTION Vars_Existence
+  !
+  LOGICAL FUNCTION DimUnlimited_Var( varname , filename )
+    !
+    CHARACTER(*),INTENT(in) :: varname,filename
+    INTEGER :: ji, ncid, varid, status, numDims, unlimDimID
+    INTEGER, DIMENSION(nf90_max_var_dims) :: VarDimIds
+    !      
+    status = nf90_open(TRIM(filename),NF90_NOWRITE,ncid)
+    IF (status/=nf90_noerr) THEN    
+       WRITE(*,*)"unable to open netcdf file : ",TRIM(filename)
+       STOP
+    ENDIF
+    status = nf90_inquire(ncid, unlimiteddimid = unlimdimid)
+    !
+    status = nf90_inq_varid(ncid,TRIM(varname),varid)
+    !
+    status = nf90_inquire_variable(ncid, varid , ndims = numdims)
+    status = nf90_inquire_variable(ncid, varid , dimids = vardimids(:numdims))
+
+    DimUnlimited_Var = .FALSE.
+    DO ji = 1, numdims
+       IF( vardimids(ji) == unlimdimid ) DimUnlimited_Var = .TRUE.
+    ENDDO
+
+    RETURN
+    !
+  END FUNCTION DimUnlimited_Var
+  
 END MODULE io_netcdf