Index: /utils/tools_MERGE_2019/NESTING/src/agrif_connect_topo.f90
===================================================================
--- /utils/tools_MERGE_2019/NESTING/src/agrif_connect_topo.f90	(revision 12411)
+++ /utils/tools_MERGE_2019/NESTING/src/agrif_connect_topo.f90	(revision 12412)
@@ -251,18 +251,18 @@
     !
     !
-    DO j = 1,jpj
-       DO i = 1,jpi
+    DO jj = 1,jpj
+       DO ji = 1,jpi
           !
-          IF (Grid%bathy_meter(i,j) .EQ. 0.0 ) THEN
-             Grid%bathy_level(i,j)=0.
+          IF (Grid%bathy_meter(ji,jj) .EQ. 0.0 ) THEN
+             Grid%bathy_level(ji,jj)=0
           ELSE
              !	
-             k1=4
+             k1=2  ! clem: minimum levels = 2 ???
              DO WHILE (k1 .LT. (N-1))
-                IF ((Grid%bathy_meter(i,j).GE.gdepw(k1)) &
-                     .AND.(Grid%bathy_meter(i,j).LE.gdepw(k1+1))) EXIT
+                IF ((Grid%bathy_meter(ji,jj).GE.gdepw(k1)) &
+                     .AND.(Grid%bathy_meter(ji,jj).LE.gdepw(k1+1))) EXIT
                 k1=k1+1
              END DO
-             Grid%bathy_level(i,j)=k1
+             Grid%bathy_level(ji,jj)=k1
              !
           ENDIF
Index: /utils/tools_MERGE_2019/NESTING/src/agrif_create_bathy.f90
===================================================================
--- /utils/tools_MERGE_2019/NESTING/src/agrif_create_bathy.f90	(revision 12411)
+++ /utils/tools_MERGE_2019/NESTING/src/agrif_create_bathy.f90	(revision 12412)
@@ -84,197 +84,82 @@
      STOP
   ENDIF
-  !                                                   ! ------------------------------------------------------------------
-  IF( .NOT.new_topo .AND. .NOT.partial_steps ) THEN   ! **** First option : no new topo file & no partial steps
-     !                                                !                 ==> interpolate levels directly from parent file
-     !                                                ! ------------------------------------------------------------------
-     WRITE(*,*) '*** First option: no new topo file & no partial steps'
-     !      
-     ! read coarse grid coordinates
-     status = Read_Coordinates(TRIM(parent_coordinate_file),G0)    
-
-     ! read coarse grid bathymetry
-     IF( TRIM(parent_bathy_level) /= '' ) THEN
+  !
+  ! read fine and coarse grids coordinates file
+  status = Read_Coordinates(TRIM(parent_coordinate_file),G0)
+  status = Read_Coordinates(TRIM(child_coordinates),G1,Pacifique)
+  !
+  ! check error in size
+  IF( imax > SIZE(G0%nav_lon,1) .OR. jmax > SIZE(G0%nav_lon,2) .OR. imax <= imin .OR. jmax <= jmin ) THEN
+     WRITE(*,*) 'ERROR ***** bad child grid definition ...'
+     WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'       
+     STOP
+  ENDIF
+  IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN
+     WRITE(*,*) 'ERROR ***** bad child coordinates file ...'
+     WRITE(*,*) 'please check that child coordinates file has been created with the same namelist'
+     STOP
+  ENDIF
+  !
+  ! read bathymetry data set => G0%bathy_meter
+  IF( new_topo ) THEN                                       ! read G0%bathy_meter (on a reduced grid) and G1 coordinates
+     DEALLOCATE( G0%nav_lon, G0%nav_lat )
+     status = read_bathy_coord(TRIM(elevation_database),G0,G1,Pacifique)
+  ELSE                                                      ! read G0%bathy_meter (on the global grid)
+     IF( TRIM(parent_bathy_meter) /= '') THEN
+        status = read_bathy_meter(TRIM(parent_bathy_meter),G0)
+     ELSE
         status = Read_bathy_level(TRIM(parent_bathy_level),G0)
         CALL levels_to_meter(G0)
-     ELSEIF( TRIM(parent_bathy_level) == '' .AND. TRIM(parent_bathy_meter) /= '') THEN
-        status = read_bathy_meter(TRIM(parent_bathy_meter),G0)
-        CALL meter_to_levels(G0)
-     ENDIF
-     !           
-     ! read fine grid coordinates
-     status = Read_Coordinates(TRIM(child_coordinates),G1,pacifique)
-     !
-     ! stop if error in size
-     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 ...'
-        WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'       
-        STOP
-     ENDIF
-     IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN
-        WRITE(*,*) 'ERROR ***** bad child coordinates file ...'
-        WRITE(*,*) 'please check that child coordinates file has been created with the same namelist'
-        STOP
-     ENDIF
-     !
-     jpj = SIZE(G0%nav_lat,2)
-     jpi = SIZE(G0%nav_lat,1)    
-     !           
-     ! create logical array masksrc
-     ALLOCATE(masksrc(jpi,jpj))
-     WHERE(G0%bathy_level.LE.0)   ;   masksrc = .FALSE.   ;
-     ELSEWHERE                    ;   masksrc = .TRUE.    ;
-     END WHERE
-
+     ENDIF     
      ! change longitudes (from -180:180 to 0:360)
-     IF ( Pacifique ) THEN
+     IF(Pacifique) THEN
         WHERE(G0%nav_lon < 0.001)   G0%nav_lon = G0%nav_lon + 360.
      ENDIF
-     !             
-     ALLOCATE(G1%bathy_meter(nxfin,nyfin))
-     !
-     ! compute remapping matrix thanks to SCRIP package
-     CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add)
-     CALL make_remap(G0%bathy_meter,G1%bathy_meter,nxfin,nyfin,matrix,src_add,dst_add)  
-     DEALLOCATE(masksrc)
-     !      
-     ! compute constant bathymetry for Parent-Child bathymetry connection
-     CALL init_constant_bathy(G0%bathy_meter,bathy_fin_constant)
-     !
-     ! replace child bathymetry by parent bathymetry at the boundaries
-     IF( ln_agrif_domain ) THEN
-        boundary = npt_copy*irafx + nbghostcellsfine + 1 
-     ELSE
-        boundary = npt_copy*irafx
-     ENDIF
-     G1%bathy_meter(1:boundary,:) = bathy_fin_constant(1:boundary,:)
-     G1%bathy_meter(:,1:boundary) = bathy_fin_constant(:,1:boundary)
-     G1%bathy_meter(nxfin-boundary+1:nxfin,:) = bathy_fin_constant(nxfin-boundary+1:nxfin,:)
-     G1%bathy_meter(:,nyfin-boundary+1:nyfin) = bathy_fin_constant(:,nyfin-boundary+1:nyfin)
-     DEALLOCATE(bathy_fin_constant)
-     !                  
-     ! bathymetry smoothing (everywhere except at the boundaries) 
-     IF( smoothing ) THEN
-        IF( ln_agrif_domain ) THEN
-           CALL smooth_topo(G1%bathy_meter(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter)
-        ELSE
-           CALL smooth_topo(G1%bathy_meter(boundary+1:nxfin-boundary,boundary+1:nyfin-boundary),nbiter)
-        ENDIF
-     ENDIF
-     !
-     ! From meters to levels
-     CALL meter_to_levels(G1)
-     G1%bathy_level=NINT(G1%bathy_level)
-     !	    
-     ! Check errors in bathy
-     DO jj=1,nyfin
-        DO ji=1,nxfin
-           IF (G1%bathy_level(ji,jj).LT.0.) THEN
-              PRINT *,'error in ',ji,jj,G1%bathy_level(ji,jj)
-              STOP
-           ENDIF
-        ENDDO
-     ENDDO
-     WHERE ((G1%bathy_level.LT.3.).AND.(G1%bathy_level.GT.0.))   G1%bathy_level=3
-     !
-     ! remove closed seas
-     IF (removeclosedseas) THEN
-        ALLOCATE(bathy_test(nxfin,nyfin))
-
-        bathy_test=0.
-        WHERE (G1%bathy_level(1,:)     .GT.0.)   bathy_test(1,:)=1
-        WHERE (G1%bathy_level(nxfin,:) .GT.0.)   bathy_test(nxfin,:)=1
-        WHERE (G1%bathy_level(:,1)     .GT.0.)   bathy_test(:,1)=1
-        WHERE (G1%bathy_level(:,nyfin) .GT.0.)   bathy_test(:,nyfin)=1
-
-        nbadd = 1
-        DO WHILE (nbadd.NE.0)
-           nbadd = 0
-           DO jj=2,nyfin-1
-              DO ji=2,nxfin-1
-                 IF (G1%bathy_level(ji,jj).GT.0.) THEN
-                    IF (MAX(bathy_test(ji,jj+1),bathy_test(ji,jj-1),bathy_test(ji-1,jj),bathy_test(ji+1,jj)).EQ.1) THEN
-                       IF (bathy_test(ji,jj).NE.1.) nbadd = nbadd + 1
-                       bathy_test(ji,jj)=1.
-                    ENDIF
-                 ENDIF
-              ENDDO
-           ENDDO
-        ENDDO
-
-        WHERE (bathy_test.EQ.0.)   G1%bathy_level = 0.
-
-        DEALLOCATE(bathy_test)
-     ENDIF
-     !
-     CALL levels_to_meter(G1) ! needed for domcfg
-     !
-     ! store interpolation result in output file
-     status = Write_Bathy_level(TRIM(child_level),G1)
-     status = write_domcfg(TRIM(child_domcfg),G1)
-     !!status = write_domcfg(TRIM(parent_domcfg_out),G0) ! do not activate it
-
-     WRITE(*,*) '****** Bathymetry successfully created if no new topo ******'
-     STOP
-     !
-     !                                                ! -----------------------------------------------------
-  ELSE                                                ! **** Second option : new topo file or partial steps     
-     !                                                ! ----------------------------------------------------- 
-     WRITE(*,*) '*** Second option : new topo or partial steps'
-
-     ! read fine and coarse grids coordinates file
-     status = Read_Coordinates(TRIM(parent_coordinate_file),G0)
-     status = Read_Coordinates(TRIM(child_coordinates),G1,Pacifique)
-     !
-     ! check error in size
-     IF( imax > SIZE(G0%nav_lon,1) .OR. jmax > SIZE(G0%nav_lon,2) .OR. imax <= imin .OR. jmax <= jmin ) THEN                    
-        WRITE(*,*) 'ERROR ***** bad child grid definition ...'
-        WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'       
-        STOP
-     ENDIF
-     IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN
-        WRITE(*,*) 'ERROR ***** bad child coordinates file ...'
-        WRITE(*,*) 'please check that child coordinates file has been created with the same namelist'
-        STOP
-     ENDIF
-     !
-     ! === From here on: G0 is the grid associated with the new topography (as gebco or etopo) ===
-     !
-     ! read bathymetry data set => G0%bathy_meter
-     IF( new_topo ) THEN                                       ! read G0%bathy_meter (on a reduced grid) and G1 coordinates
-        DEALLOCATE( G0%nav_lon, G0%nav_lat )
-        status = read_bathy_coord(TRIM(elevation_database),G0,G1,Pacifique)
-     ELSE                                                      ! read G0%bathy_meter (on the global grid)
-        IF( TRIM(parent_bathy_level) /= '' ) THEN
-           status = Read_bathy_level(TRIM(parent_bathy_level),G0)
-           CALL levels_to_meter(G0)
-        ELSEIF( TRIM(parent_bathy_level) == '' .AND. TRIM(parent_bathy_meter) /= '') THEN
-           status = read_bathy_meter(TRIM(parent_bathy_meter),G0)
-        ENDIF
-        IF(Pacifique) THEN
-           WHERE(G0%nav_lon < 0.0001)   G0%nav_lon = G0%nav_lon + 360.
-        ENDIF
-     ENDIF
-     !               
-     ! what type of interpolation for bathymetry
-     IF( type_bathy_interp == 0 ) THEN
-        WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Arithmetic average ...'
-     ELSE IF( type_bathy_interp == 1 ) THEN
-        WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Median average ...'
-     ELSE IF( type_bathy_interp == 2 ) THEN     
-        WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Bilinear interpolation ...'
-     ELSE     
-        WRITE(*,*) 'bad value for type_bathy_interp variable ( must be 0, 1 or 2 )'
-        STOP 
-     ENDIF
-     !
-     ! 1st allocation of child grid bathy
-     ALLOCATE(G1%bathy_meter(nxfin,nyfin))
-     G1%bathy_meter(:,:)=0.                       
-     !
-     ! ---------------------------------------------------------------------------------
-     ! ===                 Bathymetry of the fine grid (step1)                       ===
-     ! ---------------------------------------------------------------------------------
-     ! ==> It gives G1%bathy_meter from G0%bathy_meter
-     ! ---------------------------------------------------------------------------------
+  ENDIF
+  !
+  ! 1st allocation of child grid bathy
+  ALLOCATE(G1%bathy_meter(nxfin,nyfin))
+  G1%bathy_meter(:,:)=0.
+
+  ! 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
+  
+  IF( .NOT.new_topo )   type_bathy_interp = 2   ! only one which works
+  !
+  !
+  ! what type of interpolation for bathymetry
+  IF( type_bathy_interp == 0 ) THEN
+     WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Arithmetic average ...'
+  ELSE IF( type_bathy_interp == 1 ) THEN
+     WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Median average ...'
+  ELSE IF( type_bathy_interp == 2 ) THEN     
+     WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Bilinear interpolation ...'
+  ELSE     
+     WRITE(*,*) 'bad value for type_bathy_interp variable ( must be 0, 1 or 2 )'
+     STOP 
+  ENDIF
+  !
+  !
+  ! ---------------------------------------------------------------------------------
+  ! ===                 Bathymetry of the fine grid (step1)                       ===
+  ! ---------------------------------------------------------------------------------
+  ! ==> It gives G1%bathy_meter from G0%bathy_meter
+  ! ---------------------------------------------------------------------------------
+  
+  ! === Here: G0 is the grid associated with the new topography (as gebco or etopo) ===
+  
+  IF( .NOT. identical_grids ) THEN
      !                                                               ! ----------------------------- 
      IF( type_bathy_interp == 0 .OR. type_bathy_interp == 1 ) THEN   ! arithmetic or median averages
@@ -375,258 +260,243 @@
            type_bathy_interp = 2
         ENDIF
-           
+
+        DEALLOCATE(trouble_points)
+
      ENDIF
      !                                                       ! ----------------------------- 
      IF( type_bathy_interp == 2) THEN                        ! Bilinear interpolation
         !                                                    ! ----------------------------- 
-        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
-              PRINT*,'same grid between ',elevation_database,' and child domain'    
-              G1%bathy_meter = G0%bathy_meter 
-              identical_grids = .TRUE.                          
-           ENDIF
-        ENDIF
-
-        IF( .NOT. identical_grids ) THEN 
-
-           ALLOCATE(masksrc(SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2)))
-           ALLOCATE(bathy_test(nxfin,nyfin))
-           !
-           !Where(G0%bathy_meter.le.0.00001) 
-           !  masksrc = .false.
-           !ElseWhere
-              masksrc = .TRUE.
-           !End where                       
-           !            
-           ! compute remapping matrix thanks to SCRIP package            
-           CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add)
-           CALL make_remap(G0%bathy_meter,bathy_test,nxfin,nyfin,matrix,src_add,dst_add)  
-           !                                  
-           G1%bathy_meter = bathy_test               
-           !            
-           DEALLOCATE(masksrc)
-           DEALLOCATE(bathy_test) 
-
-        ENDIF
+
+        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.    ;
+        END WHERE
         !            
-     ENDIF
-     ! ---
-     ! At this stage bathymetry in meters has already been interpolated on fine grid
-     !                    => G1%bathy_meter(nxfin,nyfin)
+        ! compute remapping matrix thanks to SCRIP package            
+        CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add)
+        CALL make_remap(G0%bathy_meter,bathy_test,nxfin,nyfin,matrix,src_add,dst_add)  
+        !                                  
+        G1%bathy_meter = bathy_test               
+        !            
+        DEALLOCATE(masksrc)
+        DEALLOCATE(bathy_test) 
+
+     ENDIF
+     !            
+  ENDIF ! not identical grids
+  ! ---
+  ! At this stage bathymetry in meters has already been interpolated on fine grid
+  !                    => G1%bathy_meter(nxfin,nyfin)
+  !
+  ! Also G0 was the grid from the new bathymetry data set (etopo, gebco...) and not the coarse grid
+  ! ---
+  !
+  ! ---------------------------------------------------------------------------------
+  ! ===                 Bathymetry of the fine grid (step2)                       ===
+  ! ---------------------------------------------------------------------------------
+  ! ==> It gives an update of G1%bathy_meter and G1%bathy_level
+  ! ---------------------------------------------------------------------------------
+  ! From here on: G0 is the coarse grid
+  !
+  ! Coarse grid bathymetry : G0%bathy_meter (on the global grid)
+  IF( TRIM(parent_bathy_meter) /= '') THEN
+     status = read_bathy_meter(TRIM(parent_bathy_meter),G0)
+  ELSE
+     status = Read_bathy_level(TRIM(parent_bathy_level),G0)
+     CALL levels_to_meter(G0)
+  ENDIF
+
+  ! Coarse grid coordinatees : G0 coordinates
+  DEALLOCATE(G0%nav_lat,G0%nav_lon)
+  status = Read_coordinates(TRIM(parent_coordinate_file),G0)
+
+  ! allocate temporary arrays                  
+  IF (.NOT.ASSOCIATED(G0%gdepw_ps))       ALLOCATE(G0%gdepw_ps    (SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2)))
+  IF (.NOT.ASSOCIATED(G1%gdepw_ps))       ALLOCATE(G1%gdepw_ps    (SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2)))
+  IF (.NOT.ASSOCIATED(gdepw_ps_interp))   ALLOCATE(gdepw_ps_interp(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2)))
+  !                       
+  IF( ln_agrif_domain ) THEN
+     boundary = npt_copy*irafx + nbghostcellsfine + 1
+  ELSE
+     boundary = npt_copy*irafx
+  ENDIF
+  !
+  ! compute G0%gdepw_ps and G1%gdepw_ps
+  CALL get_partial_steps(G0) 
+  CALL get_partial_steps(G1)
+  CALL bathymetry_control(G0%Bathy_level)
+  
+  ! ---------------------------------------
+  ! Bathymetry at the boundaries (npt_copy)                      
+  ! ---------------------------------------
+  ! 1st step: interpolate coarse bathy on the fine grid (using partial steps or not)
+  IF( ln_agrif_domain ) THEN                   
+     CALL Check_interp(G0,gdepw_ps_interp)
+  ELSE
+     gdepw_ps_interp = 0. * G1%gdepw_ps
+     !!CALL agrif_interp(G0%gdepw_ps,gdepw_ps_interp,'T')
+     CALL init_constant_bathy(G0%gdepw_ps,gdepw_ps_interp)
+  ENDIF
+
+  IF (.NOT.ASSOCIATED(G1%wgt))   ALLOCATE(G1%wgt(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2)))
+  G1%wgt(:,:) = 0.
+  IF ((.NOT.ASSOCIATED(G0%wgt)).AND.bathy_update) THEN 
+     ALLOCATE(G0%wgt(SIZE(G0%nav_lat,1),SIZE(G0%nav_lat,2)))
+     G0%wgt(:,:) = 0.
+  ENDIF
+  !
+!!$  IF( new_topo ) THEN ! clem: no, do it even when there is no new topo
+  ! 2nd step: copy parent bathymetry at the boundaries
+  DO jj=1,nyfin   ! West and East
+     IF ( gdepw_ps_interp(nbghostcellsfine+1,jj) > 0. ) THEN
+        G1%gdepw_ps(1:boundary,jj) = gdepw_ps_interp(1:boundary,jj) 
+        G1%wgt(1:boundary,jj) = 1.
+     ELSE
+        G1%gdepw_ps(1:nbghostcellsfine+1,jj)=0. 
+     ENDIF
      !
-     ! Also G0 was the grid from the new bathymetry data set (etopo, gebco...) and not the coarse grid
-     ! ---
+     IF ( gdepw_ps_interp(nxfin-nbghostcellsfine,jj) > 0.) THEN
+        G1%gdepw_ps(nxfin-boundary+1:nxfin,jj)=gdepw_ps_interp(nxfin-boundary+1:nxfin,jj)
+        G1%wgt(nxfin-boundary+1:nxfin,jj) = 1.
+     ELSE
+        G1%gdepw_ps(nxfin-nbghostcellsfine:nxfin,jj) = 0.
+     ENDIF
+  END DO
+  !
+  DO ji=1,nxfin    ! South and North
+     IF (gdepw_ps_interp(ji,nbghostcellsfine+1)>0.) THEN
+        G1%gdepw_ps(ji,1:boundary) = gdepw_ps_interp(ji,1:boundary)
+        G1%wgt(ji,1:boundary) = 1.
+     ELSE
+        G1%gdepw_ps(ji,1:nbghostcellsfine+1)=0. 
+     ENDIF
      !
-     ! ---------------------------------------------------------------------------------
-     ! ===                 Bathymetry of the fine grid (step2)                       ===
-     ! ---------------------------------------------------------------------------------
-     ! ==> It gives an update of G1%bathy_meter and G1%bathy_level
-     ! ---------------------------------------------------------------------------------
-     ! From here on: G0 is the coarse grid
-     !
-     ! Coarse grid bathymetry : G0%bathy_meter (on the global grid)
-     IF( TRIM(parent_bathy_meter) /= '') THEN
-        status = read_bathy_meter(TRIM(parent_bathy_meter),G0)
+     IF (gdepw_ps_interp(ji,nyfin-nbghostcellsfine)>0.) THEN
+        G1%gdepw_ps(ji,nyfin-boundary+1:nyfin)=gdepw_ps_interp(ji,nyfin-boundary+1:nyfin)
+        G1%wgt(ji,nyfin-boundary+1:nyfin) = 1.
      ELSE
-        status = Read_bathy_level(TRIM(parent_bathy_level),G0)
-        CALL levels_to_meter(G0)
-     ENDIF
-     
-     ! Coarse grid coordinatees : G0 coordinates
-     DEALLOCATE(G0%nav_lat,G0%nav_lon)
-     status = Read_coordinates(TRIM(parent_coordinate_file),G0)
-
-     ! allocate temporary arrays                  
-     IF (.NOT.ASSOCIATED(G0%gdepw_ps))       ALLOCATE(G0%gdepw_ps    (SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2)))
-     IF (.NOT.ASSOCIATED(G1%gdepw_ps))       ALLOCATE(G1%gdepw_ps    (SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2)))                  
-     IF (.NOT.ASSOCIATED(gdepw_ps_interp))   ALLOCATE(gdepw_ps_interp(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2)))
-     !                       
-     IF( ln_agrif_domain ) THEN
-        boundary = npt_copy*irafx + nbghostcellsfine + 1
-     ELSE
-        boundary = npt_copy*irafx
-     ENDIF
-     !
-     ! compute G0%gdepw_ps and G1%gdepw_ps
-     CALL get_partial_steps(G0) 
-     CALL get_partial_steps(G1)
-     CALL bathymetry_control(G0%Bathy_level)
-
-     ! ---------------------------------------
-     ! Bathymetry at the boundaries (npt_copy)                      
-     ! ---------------------------------------
-     ! 1st step: interpolate coarse bathy on the fine grid (using partial steps or not)
-     IF( partial_steps .AND. ln_agrif_domain ) THEN                   
-        CALL Check_interp(G0,gdepw_ps_interp)
-     ELSE
-        gdepw_ps_interp = 0. * G1%gdepw_ps
-        !!CALL agrif_interp(G0%gdepw_ps,gdepw_ps_interp,'T')
-        CALL init_constant_bathy(G0%gdepw_ps,gdepw_ps_interp)
-     ENDIF
-
-     IF (.NOT.ASSOCIATED(G1%wgt))   ALLOCATE(G1%wgt(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2)))
-     G1%wgt(:,:) = 0.
-     IF ((.NOT.ASSOCIATED(G0%wgt)).AND.bathy_update) THEN 
-        ALLOCATE(G0%wgt(SIZE(G0%nav_lat,1),SIZE(G0%nav_lat,2)))
-        G0%wgt(:,:) = 0.
-     ENDIF
-     !
-     ! 2nd step: copy parent bathymetry at the boundaries
-     DO jj=1,nyfin   ! West and East
-        IF ( gdepw_ps_interp(nbghostcellsfine+1,jj) > 0. ) THEN
-           G1%gdepw_ps(1:boundary,jj) = gdepw_ps_interp(1:boundary,jj) 
-           G1%wgt(1:boundary,jj) = 1.
-        ELSE
-           G1%gdepw_ps(1:nbghostcellsfine+1,jj)=0. 
-        ENDIF
-        !
-        IF ( gdepw_ps_interp(nxfin-nbghostcellsfine,jj) > 0.) THEN
-           G1%gdepw_ps(nxfin-boundary+1:nxfin,jj)=gdepw_ps_interp(nxfin-boundary+1:nxfin,jj)
-           G1%wgt(nxfin-boundary+1:nxfin,jj) = 1.
-        ELSE
-           G1%gdepw_ps(nxfin-nbghostcellsfine:nxfin,jj) = 0.
-        ENDIF
+        G1%gdepw_ps(ji,nyfin-nbghostcellsfine:nyfin) = 0.
+     ENDIF
+  END DO
+  !
+  !clem: recalculate interpolation everywhere before linear connection (useless to me??)
+  IF( ln_agrif_domain ) THEN                
+     gdepw_ps_interp = 0.
+     CALL Check_interp(G0,gdepw_ps_interp)
+  ENDIF
+  !
+  ! -------------------------------------------------------
+  ! Bathymetry between boundaries and interior (npt_connect)                 
+  ! --------------------------------------------------------
+  ! Make linear connection (on npt_connect*irafx points) between the boundaries and the interior
+  IF( ln_agrif_domain ) THEN
+     boundary = (npt_copy + npt_connect)*irafx + nbghostcellsfine + 1
+  ELSE
+     boundary = (npt_copy + npt_connect)*irafx
+  ENDIF
+
+  IF( npt_connect > 0 ) THEN
+     WRITE(*,*) ' linear connection on ',npt_connect,'coarse grid points'
+
+     wghts = 1.
+     DO ji = boundary - npt_connect*irafx + 1 , boundary
+        wghts = wghts - (1. / (npt_connect*irafx + 1. ) )
+        DO jj=1,nyfin
+           IF (G1%gdepw_ps(nbghostcellsfine+1,jj) > 0.)       G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))  
+        END DO
      END DO
-     !
-     DO ji=1,nxfin    ! South and North
-        IF (gdepw_ps_interp(ji,nbghostcellsfine+1)>0.) THEN
-           G1%gdepw_ps(ji,1:boundary) = gdepw_ps_interp(ji,1:boundary)
-           G1%wgt(ji,1:boundary) = 1.
-        ELSE
-           G1%gdepw_ps(ji,1:nbghostcellsfine+1)=0. 
-        ENDIF
-        !
-        IF (gdepw_ps_interp(ji,nyfin-nbghostcellsfine)>0.) THEN
-           G1%gdepw_ps(ji,nyfin-boundary+1:nyfin)=gdepw_ps_interp(ji,nyfin-boundary+1:nyfin)
-           G1%wgt(ji,nyfin-boundary+1:nyfin) = 1.
-        ELSE
-           G1%gdepw_ps(ji,nyfin-nbghostcellsfine:nyfin) = 0.
-        ENDIF
+
+     wghts = 1.
+     DO ji = nxfin - (boundary - npt_connect*irafx), nxfin - boundary +1 , -1
+        wghts = wghts - (1. / (npt_connect*irafx + 1. ) )
+        DO jj=1,nyfin
+           IF (G1%gdepw_ps(nxfin-nbghostcellsfine,jj) > 0.)   G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))
+        END DO
      END DO
-     !
-     !clem: recalculate interpolation everywhere before linear connection (useless to me)
-     IF( partial_steps .AND. ln_agrif_domain ) THEN                
-        gdepw_ps_interp = 0.
-        CALL Check_interp(G0,gdepw_ps_interp)
-     ENDIF
-     !
-     ! -------------------------------------------------------
-     ! Bathymetry between boundaries and interior (npt_connect)                 
-     ! --------------------------------------------------------
-     ! Make linear connection (on npt_connect*irafx points) between the boundaries and the interior
-     IF( ln_agrif_domain ) THEN
-        boundary = (npt_copy + npt_connect)*irafx + nbghostcellsfine + 1
-     ELSE
-        boundary = (npt_copy + npt_connect)*irafx
-     ENDIF
-
-     IF( npt_connect > 0 ) THEN
-        WRITE(*,*) ' linear connection on ',npt_connect,'coarse grid points'
-
-        wghts = 1.
-        DO ji = boundary - npt_connect*irafx + 1 , boundary
-           wghts = wghts - (1. / (npt_connect*irafx + 1. ) )
-           DO jj=1,nyfin
-              IF (G1%gdepw_ps(nbghostcellsfine+1,jj) > 0.)       G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))  
-           END DO
+
+     wghts = 1.
+     DO jj = boundary - npt_connect*irafy + 1 , boundary
+        wghts = wghts - (1. / (npt_connect*irafy + 1. ) )
+        DO ji=1,nxfin
+           IF (G1%gdepw_ps(ji,nbghostcellsfine+1) > 0.)       G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))
         END DO
-
-        wghts = 1.
-        DO ji = nxfin - (boundary - npt_connect*irafx), nxfin - boundary +1 , -1
-           wghts = wghts - (1. / (npt_connect*irafx + 1. ) )
-           DO jj=1,nyfin
-              IF (G1%gdepw_ps(nxfin-nbghostcellsfine,jj) > 0.)   G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))
-           END DO
+     END DO
+
+     wghts = 1.
+     DO jj = nyfin - (boundary - npt_connect*irafy) , nyfin - boundary +1, -1
+        wghts = wghts - (1. / (npt_connect*irafy + 1. ) )
+        DO ji=1,nxfin
+           IF (G1%gdepw_ps(ji,nyfin-nbghostcellsfine) > 0.)   G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))
         END DO
-
-        wghts = 1.
-        DO jj = boundary - npt_connect*irafy + 1 , boundary
-           wghts = wghts - (1. / (npt_connect*irafy + 1. ) )
-           DO ji=1,nxfin
-              IF (G1%gdepw_ps(ji,nbghostcellsfine+1) > 0.)       G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))
-           END DO
-        END DO
-
-        wghts = 1.
-        DO jj = nyfin - (boundary - npt_connect*irafy) , nyfin - boundary +1, -1
-           wghts = wghts - (1. / (npt_connect*irafy + 1. ) )
-           DO ji=1,nxfin
-              IF (G1%gdepw_ps(ji,nyfin-nbghostcellsfine) > 0.)   G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))
-           END DO
-        END DO
-        IF (.NOT.identical_grids) THEN
-           G1%gdepw_ps(:,:) = (1.-G1%wgt(:,:)) * G1%gdepw_ps(:,:) + gdepw_ps_interp(:,:)*G1%wgt(:,:)
-        ENDIF
-
-     ENDIF
-
-     ! replace G1%bathy_meter by G1%gdepw_ps
-     G1%bathy_meter = G1%gdepw_ps
-     !                     
-     ! --------------------
-     ! Bathymetry smoothing                 
-     ! --------------------
-     IF( smoothing .AND. (.NOT.identical_grids) ) THEN
-        ! Chanut: smoothing everywhere then discard result in connection zone
-        CALL smooth_topo(G1%gdepw_ps(:,:),nbiter)
-        WHERE (G1%wgt(:,:)==0.) G1%bathy_meter(:,:) = G1%gdepw_ps(:,:)
-     ELSE
-        WRITE(*,*) 'No smoothing process only connection is carried out'
-     ENDIF
-     !
-     ! ------------------
-     ! Remove closed seas
-     ! ------------------
-     IF (removeclosedseas) THEN
-        ALLOCATE(bathy_test(nxfin,nyfin))
-        bathy_test=0.
-        WHERE (G1%bathy_meter(1,:)    .GT.0.)   bathy_test(1,:)=1
-        WHERE (G1%bathy_meter(nxfin,:).GT.0.)   bathy_test(nxfin,:)=1
-        WHERE (G1%bathy_meter(:,1)    .GT.0.)   bathy_test(:,1)=1
-        WHERE (G1%bathy_meter(:,nyfin).GT.0.)   bathy_test(:,nyfin)=1
-
-        nbadd = 1
-        DO WHILE (nbadd.NE.0)
-           nbadd = 0
-           DO jj=2,nyfin-1
-              DO ji=2,nxfin-1
-                 IF (G1%bathy_meter(ji,jj).GT.0.) THEN
-                    IF (MAX(bathy_test(ji,jj+1),bathy_test(ji,jj-1),bathy_test(ji-1,jj),bathy_test(ji+1,jj)).EQ.1) THEN
-                       IF (bathy_test(ji,jj).NE.1.) nbadd = nbadd + 1
-                       bathy_test(ji,jj)=1.
-                    ENDIF
-
+     END DO
+     IF (.NOT.identical_grids) THEN
+        G1%gdepw_ps(:,:) = (1.-G1%wgt(:,:)) * G1%gdepw_ps(:,:) + gdepw_ps_interp(:,:)*G1%wgt(:,:)
+     ENDIF
+
+  ENDIF
+!!$  ENDIF
+
+  ! replace G1%bathy_meter by G1%gdepw_ps
+  G1%bathy_meter = G1%gdepw_ps
+  !                     
+  ! --------------------
+  ! Bathymetry smoothing                 
+  ! --------------------
+  IF( smoothing .AND. (.NOT.identical_grids) ) THEN
+     ! Chanut: smoothing everywhere then discard result in connection zone
+     CALL smooth_topo(G1%gdepw_ps(:,:),nbiter)
+     WHERE (G1%wgt(:,:)==0.) G1%bathy_meter(:,:) = G1%gdepw_ps(:,:)
+  ELSE
+     WRITE(*,*) 'No smoothing process only connection is carried out'
+  ENDIF
+  !
+  ! ------------------
+  ! Remove closed seas
+  ! ------------------
+  IF (removeclosedseas) THEN
+     ALLOCATE(bathy_test(nxfin,nyfin))
+     bathy_test=0.
+     WHERE (G1%bathy_meter(1,:)    .GT.0.)   bathy_test(1,:)=1
+     WHERE (G1%bathy_meter(nxfin,:).GT.0.)   bathy_test(nxfin,:)=1
+     WHERE (G1%bathy_meter(:,1)    .GT.0.)   bathy_test(:,1)=1
+     WHERE (G1%bathy_meter(:,nyfin).GT.0.)   bathy_test(:,nyfin)=1
+
+     nbadd = 1
+     DO WHILE (nbadd.NE.0)
+        nbadd = 0
+        DO jj=2,nyfin-1
+           DO ji=2,nxfin-1
+              IF (G1%bathy_meter(ji,jj).GT.0.) THEN
+                 IF (MAX(bathy_test(ji,jj+1),bathy_test(ji,jj-1),bathy_test(ji-1,jj),bathy_test(ji+1,jj)).EQ.1) THEN
+                    IF (bathy_test(ji,jj).NE.1.) nbadd = nbadd + 1
+                    bathy_test(ji,jj)=1.
                  ENDIF
-              ENDDO
+
+              ENDIF
            ENDDO
         ENDDO
-        WHERE (bathy_test.EQ.0.)   G1%bathy_meter = 0.
-        DEALLOCATE(bathy_test)
-     ENDIF
-     !
-     IF( partial_steps ) THEN
-        CALL get_partial_steps(G1)  ! recompute bathy_level and gdepw_ps for G1 (and correct bathy_meter)
-     ELSE
-        CALL meter_to_levels(G1)    ! convert bathymetry from meters to levels
-     ENDIF
-     !
-     ! update parent grid
-     IF(bathy_update)   CALL Update_Parent_Bathy( G0,G1 ) 
-     
-     IF(bathy_update)   status = Write_Bathy_meter(TRIM(parent_bathy_meter_updated),G0)
-     IF(bathy_update)   status = write_domcfg(TRIM(parent_domcfg_updated),G0)
-     
-     ! store interpolation result in output file
-     IF( TRIM(parent_bathy_level) /= '' )   status = Write_Bathy_level(TRIM(child_level),G1)
-     IF( TRIM(parent_bathy_meter) /= '' )   status = Write_Bathy_meter(TRIM(child_meter),G1)
-     IF( TRIM(parent_domcfg_out)  /= '' )   status = write_domcfg(TRIM(child_domcfg),G1)
-     !
-     WRITE(*,*) '****** Bathymetry successfully created ******'
-     STOP
-  ENDIF
+     ENDDO
+     WHERE (bathy_test.EQ.0.)   G1%bathy_meter = 0.
+     DEALLOCATE(bathy_test)
+  ENDIF
+  !
+  CALL get_partial_steps(G1)  ! recompute bathy_level and gdepw_ps for G1 (and correct bathy_meter)
+  !
+  ! update parent grid
+  IF(bathy_update) THEN
+     CALL Update_Parent_Bathy( G0,G1 ) 
+     status = Write_Bathy_meter(TRIM(parent_bathy_meter_updated),G0)
+     status = write_domcfg(TRIM(parent_domcfg_updated),G0)
+  ENDIF
+  !
+  ! store interpolation result in output file
+  IF( TRIM(parent_bathy_level) /= '' )   status = Write_Bathy_level(TRIM(child_level),G1)
+  IF( TRIM(parent_bathy_meter) /= '' )   status = Write_Bathy_meter(TRIM(child_meter),G1)
+  IF( TRIM(parent_domcfg_out)  /= '' )   status = write_domcfg(TRIM(child_domcfg),G1)
+  !
+  WRITE(*,*) '****** Bathymetry successfully created ******'
+  STOP
   !
 END PROGRAM create_bathy
Index: /utils/tools_MERGE_2019/NESTING/src/agrif_interpolation.f90
===================================================================
--- /utils/tools_MERGE_2019/NESTING/src/agrif_interpolation.f90	(revision 12411)
+++ /utils/tools_MERGE_2019/NESTING/src/agrif_interpolation.f90	(revision 12412)
@@ -198,10 +198,13 @@
     CHARACTER(*) :: typevar   
 
-    INTEGER :: nxf,nyf,zx,zy
-    INTEGER :: ji,jj,jif,jjf,jic,jjc,jdecx,jdecy
+    INTEGER :: nxf,nyf,nxc,nyc,zx,zy
+    INTEGER :: ji,jj,jif,jjf,jic,jjc,jic1,jjc1,jdecx,jdecy
     REAL*8 :: Ax, Bx, Ay, By
 
     nxf = SIZE(tabout,1)
     nyf = SIZE(tabout,2)
+
+    nxc = SIZE(tabin,1)
+    nyc = SIZE(tabin,2)
 
     SELECT CASE(typevar)
@@ -248,6 +251,9 @@
           Ay = 1. - By
 
-          tabout(ji,jj) = ( Bx * tabin(jic+1,jjc  ) + Ax * tabin(jic,jjc  ) ) * Ay + &
-             &            ( Bx * tabin(jic+1,jjc+1) + Ax * tabin(jic,jjc+1) ) * By
+          jic1 = MIN( nxc, jic+1 ) ! avoid out of bounds for tabin below
+          jjc1 = MIN( nyc, jjc+1 ) !             --
+
+          tabout(ji,jj) = ( Bx * tabin(jic1,jjc ) + Ax * tabin(jic,jjc ) ) * Ay + &
+             &            ( Bx * tabin(jic1,jjc1) + Ax * tabin(jic,jjc1) ) * By
        END DO
     END DO
Index: /utils/tools_MERGE_2019/NESTING/src/agrif_partial_steps.f90
===================================================================
--- /utils/tools_MERGE_2019/NESTING/src/agrif_partial_steps.f90	(revision 12411)
+++ /utils/tools_MERGE_2019/NESTING/src/agrif_partial_steps.f90	(revision 12412)
@@ -174,9 +174,9 @@
     END DO
     !
+!!$    IF( partial_steps ) THEN
     ! Compute bathy_level for ocean points (i.e. the number of ocean levels)
     ! find the number of ocean levels such that the last level thickness
     ! is larger than the minimum of e3zps_min and e3zps_rat * e3t (where
     ! e3t is the reference level thickness   
-    !     
     DO jk = N-1, 1, -1
        zdepth = gdepw(jk) + MIN( e3zps_min, e3t(jk)*e3zps_rat )
@@ -184,9 +184,30 @@
           DO ji = 1, jpi
              IF( 0. < Grid%bathy_meter(ji,jj) .AND. Grid%bathy_meter(ji,jj) <= zdepth ) &
-                  Grid%bathy_level(ji,jj) = jk-1
+                Grid%bathy_level(ji,jj) = jk-1
           END DO
        END DO
     END DO
-
+!!$    ELSE
+!!$       DO jj = 1,jpj
+!!$          DO ji = 1,jpi
+!!$             !
+!!$             IF (Grid%bathy_meter(ji,jj) .EQ. 0.0 ) THEN
+!!$                Grid%bathy_level(ji,jj)=0
+!!$             ELSE
+!!$                !	
+!!$                k1=2  ! clem: minimum levels = 2 ???
+!!$                DO WHILE (k1 .LT. (N-1))
+!!$                   IF ((Grid%bathy_meter(ji,jj).GE.gdepw(k1)) &
+!!$                      .AND.(Grid%bathy_meter(ji,jj).LE.gdepw(k1+1))) EXIT
+!!$                   k1=k1+1
+!!$                END DO
+!!$                Grid%bathy_level(ji,jj)=k1
+!!$                !
+!!$             ENDIF
+!!$             !
+!!$          END DO
+!!$       END DO
+!!$
+!!$    ENDIF
 
     CALL bathymetry_control(grid%bathy_level)
@@ -675,5 +696,5 @@
   SUBROUTINE bathymetry_control(mbathy)
 
-    INTEGER ::   i, j, jl           
+    INTEGER ::   ji, jj, jl           
     INTEGER ::   icompt, ibtest, ikmax          
     REAL*8, DIMENSION(:,:) :: mbathy     
@@ -692,13 +713,13 @@
     DO jl = 1, 2
        !
-       DO j = 2, SIZE(mbathy,2)-1
-          DO i = 2, SIZE(mbathy,1)-1
-
-             ibtest = MAX( mbathy(i-1,j), mbathy(i+1,j),mbathy(i,j-1),mbathy(i,j+1) )
+       DO jj = 2, SIZE(mbathy,2)-1
+          DO ji = 2, SIZE(mbathy,1)-1
+
+             ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj),mbathy(ji,jj-1),mbathy(ji,jj+1) )
              !               
-             IF( ibtest < mbathy(i,j) ) THEN
+             IF( ibtest < mbathy(ji,jj) ) THEN
                 !                  
-                WRITE(*,*) 'grid-point(i,j)= ',i,j,'is changed from',mbathy(i,j),' to ', ibtest
-                mbathy(i,j) = ibtest
+                WRITE(*,*) 'grid-point(i,j)= ',ji,jj,'is changed from',mbathy(ji,jj),' to ', ibtest
+                mbathy(ji,jj) = ibtest
                 icompt = icompt + 1
                 !
@@ -720,7 +741,7 @@
 
     ikmax = 0
-    DO j = 1, SIZE(mbathy,2)
+    DO jj = 1, SIZE(mbathy,2)
        DO ji = 1, SIZE(mbathy,1)
-          ikmax = MAX( ikmax, NINT(mbathy(i,j)) )
+          ikmax = MAX( ikmax, NINT(mbathy(ji,jj)) )
        END DO
     END DO
Index: /utils/tools_MERGE_2019/NESTING/src/agrif_readwrite.f90
===================================================================
--- /utils/tools_MERGE_2019/NESTING/src/agrif_readwrite.f90	(revision 12411)
+++ /utils/tools_MERGE_2019/NESTING/src/agrif_readwrite.f90	(revision 12412)
@@ -241,5 +241,5 @@
     CALL write_ncdf_var('nav_lon'        ,dimnames,name,Grid%nav_lon    ,'float')
     CALL write_ncdf_var('nav_lat'        ,dimnames,name,Grid%nav_lat    ,'float')
-    CALL write_ncdf_var(parent_level_name,dimnames,name,Grid%bathy_level,'float')
+    CALL write_ncdf_var(parent_level_name,dimnames,name,NINT(Grid%bathy_level),'integer')
     !
     CALL copy_ncdf_att('nav_lon'        ,TRIM(parent_bathy_level),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon))
@@ -545,13 +545,14 @@
      ln_sco = 0
      ln_isfcav = 0
-     IF( partial_steps ) THEN
+!!$     IF( partial_steps ) THEN
         ln_zps = 1
         ln_zco = 0
-        CALL zgr_zps( nx, ny, gdept_1d, gdepw_1d, e3t_1d, e3w_1d, zbathy, &        ! input
-        &             mbathy, e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 )  ! output
-     ELSE
-        ln_zps = 0
-        ln_zco = 1
-     ENDIF
+!!$     ELSE
+!!$        ln_zps = 0
+!!$        ln_zco = 1
+!!$     ENDIF
+
+     CALL zgr_zps( nx, ny, gdept_1d, gdepw_1d, e3t_1d, e3w_1d, zbathy, &        ! input
+        &          mbathy, e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 )  ! output
 
      ! top/bottom levels
@@ -865,17 +866,23 @@
       zmax = gdepw_1d(N) + e3t_1d(N)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(N-1) )
       zbathy(:,:) = MIN( zmax ,  zbathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat)
-      WHERE( zbathy(:,:) == 0. )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0
+      WHERE( zbathy(:,:) == 0. )   ;   mbathy(:,:) = 0     ! land  : set mbathy to 0
       ELSEWHERE                    ;   mbathy(:,:) = N-1   ! ocean : initialize mbathy to the max ocean level
       END WHERE
 
-      ! Compute mbathy for ocean points (i.e. the number of ocean levels)
-      ! find the number of ocean levels such that the last level thickness
-      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
-      ! e3t_1d is the reference level thickness
-      DO jk = N-1, 1, -1
-         zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
-         WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
-      END DO
-
+!!$      IF( partial_steps ) THEN
+         ! Compute mbathy for ocean points (i.e. the number of ocean levels)
+         ! find the number of ocean levels such that the last level thickness
+         ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
+         ! e3t_1d is the reference level thickness
+         DO jk = N-1, 1, -1
+            zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
+            WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
+         END DO
+!!$      ELSE
+!!$         DO jk = 1, N
+!!$            WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) >= gdepw_1d(jk) )   mbathy(:,:) = jk-1
+!!$         END DO        
+!!$      ENDIF
+      
       ! Scale factors and depth at T- and W-points
       DO jk = 1, N                        ! intitialization to the reference z-coordinate
@@ -891,54 +898,54 @@
       !
       ! Scale factors and depth at T- and W-points
-!      IF ( ln_isfcav == 0 ) THEN
-         DO jj = 1, ny
-            DO ji = 1, nx
-               ik = mbathy(ji,jj)
-               IF( ik > 0 ) THEN               ! ocean point only
-                  ! max ocean level case
-                  IF( ik == N-1 ) THEN
-                     zdepwp = zbathy(ji,jj)
-                     ze3tp  = zbathy(ji,jj) - gdepw_1d(ik)
-                     ze3wp = 0.5 * e3w_1d(ik) * ( 1. + ( ze3tp/e3t_1d(ik) ) )
-                     e3t_0(ji,jj,ik  ) = ze3tp
-                     e3t_0(ji,jj,ik+1) = ze3tp
-                     e3w_0(ji,jj,ik  ) = ze3wp
-                     e3w_0(ji,jj,ik+1) = ze3tp
-                     gdepw_0(ji,jj,ik+1) = zdepwp
-                     gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
-                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
-                     !
-                  ELSE                         ! standard case
-                     IF( zbathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = zbathy(ji,jj)
-                     ELSE                                        ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
-                     ENDIF
-   !gm Bug?  check the gdepw_1d
-                     !       ... on ik
-                     gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   &
-                        &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
-                        &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
-                     e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
-                        &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
-                     e3w_0(ji,jj,ik) = 0.5 * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2. * gdepw_1d(ik) )   &
-                        &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
-                     !       ... on ik+1
-                     e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
-                     e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
-                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
+      ! IF ( ln_isfcav == 0 ) THEN
+      DO jj = 1, ny
+         DO ji = 1, nx
+            ik = mbathy(ji,jj)
+            IF( ik > 0 ) THEN               ! ocean point only
+               ! max ocean level case
+               IF( ik == N-1 ) THEN
+                  zdepwp = zbathy(ji,jj)
+                  ze3tp  = zbathy(ji,jj) - gdepw_1d(ik)
+                  ze3wp = 0.5 * e3w_1d(ik) * ( 1. + ( ze3tp/e3t_1d(ik) ) )
+                  e3t_0(ji,jj,ik  ) = ze3tp
+                  e3t_0(ji,jj,ik+1) = ze3tp
+                  e3w_0(ji,jj,ik  ) = ze3wp
+                  e3w_0(ji,jj,ik+1) = ze3tp
+                  gdepw_0(ji,jj,ik+1) = zdepwp
+                  gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
+                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
+                  !
+               ELSE                         ! standard case
+                  IF( zbathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = zbathy(ji,jj)
+                  ELSE                                        ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
                   ENDIF
+                  !gm Bug?  check the gdepw_1d
+                  !       ... on ik
+                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   &
+                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
+                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
+                  e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
+                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
+                  e3w_0(ji,jj,ik) = 0.5 * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2. * gdepw_1d(ik) )   &
+                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
+                  !       ... on ik+1
+                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
+                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
+                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
                ENDIF
-            END DO
+            ENDIF
          END DO
-         !
-!         DO jj = 1, ny
-!            DO ji = 1, nx
-!               ik = mbathy(ji,jj)
-!               IF( ik > 0 ) THEN               ! ocean point only
-!                  e3tp (ji,jj) = e3t_0(ji,jj,ik)
-!                  e3wp (ji,jj) = e3w_0(ji,jj,ik)
-!               ENDIF
-!            END DO
-!         END DO
-!      END IF
+      END DO
+      !
+      ! DO jj = 1, ny
+      !    DO ji = 1, nx
+      !       ik = mbathy(ji,jj)
+      !       IF( ik > 0 ) THEN               ! ocean point only
+      !          e3tp (ji,jj) = e3t_0(ji,jj,ik)
+      !          e3wp (ji,jj) = e3w_0(ji,jj,ik)
+      !       ENDIF
+      !    END DO
+      ! END DO
+      ! END IF
       !
       ! Scale factors and depth at U-, V-, UW and VW-points
Index: /utils/tools_MERGE_2019/NESTING/src/agrif_types.f90
===================================================================
--- /utils/tools_MERGE_2019/NESTING/src/agrif_types.f90	(revision 12411)
+++ /utils/tools_MERGE_2019/NESTING/src/agrif_types.f90	(revision 12412)
@@ -192,4 +192,9 @@
        ENDIF
        !
+       IF( .NOT.partial_steps ) THEN
+          WRITE(*,*) 'E R R O R: partial_steps must be set to true otherwise very thin bottom layers can be created'
+          STOP
+       ENDIF
+       !
     ELSE
        !
Index: /utils/tools_MERGE_2019/SIREN/src/create_boundary.F90
===================================================================
--- /utils/tools_MERGE_2019/SIREN/src/create_boundary.F90	(revision 12411)
+++ /utils/tools_MERGE_2019/SIREN/src/create_boundary.F90	(revision 12412)
@@ -937,5 +937,5 @@
    ELSE
 
-      WRITE(cl_errormsg,*) " ERROR : can't find "//TRIM(cl_namelist)
+      WRITE(cl_errormsg,*) " ERROR : can't find "//TRIM(cd_namelist)
       CALL fct_help(cp_myname,cl_errormsg) 
       CALL EXIT(1)
Index: /utils/tools_MERGE_2019/SIREN/src/docsrc/5_changeLog.md
===================================================================
--- /utils/tools_MERGE_2019/SIREN/src/docsrc/5_changeLog.md	(revision 12411)
+++ /utils/tools_MERGE_2019/SIREN/src/docsrc/5_changeLog.md	(revision 12412)
@@ -5,10 +5,10 @@
 # Release 2019-12-03 {#rev_2019-12-03}
 ## New features
-- M Siren/src/iom_cdf.f90 :
+- M src/iom_cdf.f90 :
 	- write netcdf file as netcdf4
 
 # Release 2019-11-05 {#rev_2019-11-05}
 ## New features
-- M Siren/src/function.f90
+- M src/function.f90
 - M src/create_bathy.f90 :
 	- add help and version optional arguments
