Index: utils/tools/NESTING/src/agrif_create_bathy.f90
===================================================================
--- utils/tools/NESTING/src/agrif_create_bathy.f90	(revision 10381)
+++ utils/tools/NESTING/src/agrif_create_bathy.f90	(revision 12253)
@@ -87,4 +87,5 @@
   IF( .NOT.new_topo .AND. .NOT.partial_steps ) THEN   ! **** First option : no new topo file & no partial steps
      !                                                !                 ==> interpolate levels directly from parent file
+     !                                                !                     using bilinear interpolation
      !                                                ! ------------------------------------------------------------------
      WRITE(*,*) '*** First option: no new topo file & no partial steps'
@@ -169,5 +170,5 @@
      DO jj=1,nyfin
         DO ji=1,nxfin
-           IF (G1%bathy_level(ji,jj).LT.0.) THEN
+           IF (G1%bathy_level(ji,jj).LT.0) THEN
               PRINT *,'error in ',ji,jj,G1%bathy_level(ji,jj)
               STOP
@@ -175,5 +176,5 @@
         ENDDO
      ENDDO
-     WHERE ((G1%bathy_level.LT.3.).AND.(G1%bathy_level.GT.0.))   G1%bathy_level=3
+     WHERE ((G1%bathy_level.LT.3).AND.(G1%bathy_level.GT.0))   G1%bathy_level=3
      !
      ! remove closed seas
@@ -181,9 +182,9 @@
         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
+        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
@@ -192,8 +193,8 @@
            DO jj=2,nyfin-1
               DO ji=2,nxfin-1
-                 IF (G1%bathy_level(ji,jj).GT.0.) THEN
+                 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.
+                       IF (bathy_test(ji,jj).NE.1) nbadd = nbadd + 1
+                       bathy_test(ji,jj)=1
                     ENDIF
                  ENDIF
@@ -202,5 +203,5 @@
         ENDDO
 
-        WHERE (bathy_test.EQ.0.)   G1%bathy_level = 0.
+        WHERE (bathy_test.EQ.0.)   G1%bathy_level = 0
 
         DEALLOCATE(bathy_test)
@@ -214,5 +215,5 @@
      !!status = write_domcfg(TRIM(parent_domcfg_out),G0) ! do not activate it
 
-     WRITE(*,*) '****** Bathymetry successfully created if no new topo ******'
+     WRITE(*,*) '****** Bathymetry successfully created if no new topo and no partial steps ******'
      STOP
      !
@@ -277,120 +278,122 @@
      ! ==> It gives G1%bathy_meter from G0%bathy_meter
      ! ---------------------------------------------------------------------------------
-     !                                                               ! ----------------------------- 
-     IF( type_bathy_interp == 0 .OR. type_bathy_interp == 1 ) THEN   ! arithmetic or median averages
-        !                                                            ! ----------------------------- 
-        ALLOCATE(trouble_points(nxfin,nyfin))
-        trouble_points(:,:) = 0
-        !                       
-        DO jj = 2, nyfin
-           DO ji = 2, nxfin
-              !	    
-              ! fine grid cell extension               
-              Cell_lonmin = MIN(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj))
-              Cell_lonmax = MAX(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj))
-              Cell_latmin = MIN(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj))
-              Cell_latmax = MAX(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj)) 
-              !               
-              ! look for points in G0 (bathy dataset) contained in the fine grid cells  
-              iimin = 1
-              DO WHILE( G0%nav_lon(iimin,1) < Cell_lonmin ) 
-                 iimin = iimin + 1
+     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. identical_grids ) THEN
+        !                                                               ! ----------------------------- 
+        IF( type_bathy_interp == 0 .OR. type_bathy_interp == 1 ) THEN   ! arithmetic or median averages
+           !                                                            ! ----------------------------- 
+           ALLOCATE(trouble_points(nxfin,nyfin))
+           trouble_points(:,:) = 0
+           !                       
+           DO jj = 2, nyfin
+              DO ji = 2, nxfin
+                 !	    
+                 ! fine grid cell extension               
+                 Cell_lonmin = MIN(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj))
+                 Cell_lonmax = MAX(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj))
+                 Cell_latmin = MIN(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj))
+                 Cell_latmax = MAX(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj)) 
+                 !               
+                 ! look for points in G0 (bathy dataset) contained in the fine grid cells  
+                 iimin = 1
+                 DO WHILE( G0%nav_lon(iimin,1) < Cell_lonmin ) 
+                    iimin = iimin + 1
+                 ENDDO
+                 !               
+                 jjmin = 1
+                 DO WHILE( G0%nav_lat(iimin,jjmin) < Cell_latmin ) 
+                    jjmin = jjmin + 1
+                 ENDDO
+                 !                
+                 iimax = iimin 
+                 DO WHILE( G0%nav_lon(iimax,1) <= Cell_lonmax ) 
+                    iimax = iimax + 1
+                    iimax = MIN( iimax,SIZE(G0%bathy_meter,1))
+                 ENDDO
+                 !                               
+                 jjmax = jjmin 
+                 DO WHILE( G0%nav_lat(iimax,jjmax) <= Cell_latmax ) 
+                    jjmax = jjmax + 1
+                    jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2))
+                 ENDDO
+                 !
+                 IF( ln_agrif_domain ) THEN
+                    iimax = iimax-1
+                    jjmax = jjmax-1
+                 ELSE
+                    iimax = MAX(iimin,iimax-1)
+                    jjmax = MAX(jjmin,jjmax-1)
+                 ENDIF
+                 !               
+                 iimin = MAX( iimin,1 )
+                 jjmin = MAX( jjmin,1 )
+                 iimax = MIN( iimax,SIZE(G0%bathy_meter,1))
+                 jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2))
+
+                 nxhr = iimax - iimin + 1
+                 nyhr = jjmax - jjmin + 1                    
+
+                 IF( nxhr == 0 .OR. nyhr == 0 ) THEN
+                    !
+                    trouble_points(ji,jj) = 1
+                    !
+                 ELSE
+                    !
+                    ALLOCATE( vardep(nxhr,nyhr), mask_oce(nxhr,nyhr) )
+                    vardep(:,:) = G0%bathy_meter(iimin:iimax,jjmin:jjmax)
+                    !
+                    WHERE( vardep(:,:) .GT. 0. )   ;   mask_oce = 1 ;
+                    ELSEWHERE                      ;   mask_oce = 0 ;
+                    ENDWHERE
+                    !
+                    nxyhr = nxhr*nyhr
+                    IF( SUM(mask_oce) < 0.5*(nxyhr) ) THEN ! if more than half of the points are on land then bathy fine = 0
+                       G1%bathy_meter(ji,jj) = 0.
+                    ELSE
+                       IF( type_bathy_interp == 0 ) THEN     ! Arithmetic average
+                          G1%bathy_meter(ji,jj) = SUM( vardep(:,:) * mask_oce(:,:) ) / SUM( mask_oce(:,:) )
+                       ELSE                                  ! Median average
+                          ALLOCATE(vardep1d(nxyhr))
+                          vardep1d = RESHAPE(vardep,(/ nxyhr /) )
+                          !!CALL ssort(vardep1d,nxyhr)
+                          CALL quicksort(vardep1d,1,nxyhr)
+                          !
+                          ! Calculate median
+                          IF (MOD(nxyhr,2) .NE. 0) THEN
+                             G1%bathy_meter(ji,jj) = vardep1d( nxyhr/2 + 1 )
+                          ELSE
+                             G1%bathy_meter(ji,jj) = 0.5 * ( vardep1d(nxyhr/2) + vardep1d(nxyhr/2+1) )
+                          END IF
+                          DEALLOCATE(vardep1d)   
+                       ENDIF
+                    ENDIF
+                    DEALLOCATE (mask_oce,vardep)
+                    !
+                 ENDIF
               ENDDO
-              !               
-              jjmin = 1
-              DO WHILE( G0%nav_lat(iimin,jjmin) < Cell_latmin ) 
-                 jjmin = jjmin + 1
-              ENDDO
-              !                
-              iimax = iimin 
-              DO WHILE( G0%nav_lon(iimax,1) <= Cell_lonmax ) 
-                 iimax = iimax + 1
-                 iimax = MIN( iimax,SIZE(G0%bathy_meter,1))
-              ENDDO
-              !                               
-              jjmax = jjmin 
-              DO WHILE( G0%nav_lat(iimax,jjmax) <= Cell_latmax ) 
-                 jjmax = jjmax + 1
-                 jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2))
-              ENDDO
-              !
-              IF( ln_agrif_domain ) THEN
-                 iimax = iimax-1
-                 jjmax = jjmax-1
-              ELSE
-                 iimax = MAX(iimin,iimax-1)
-                 jjmax = MAX(jjmin,jjmax-1)
-              ENDIF
-              !               
-              iimin = MAX( iimin,1 )
-              jjmin = MAX( jjmin,1 )
-              iimax = MIN( iimax,SIZE(G0%bathy_meter,1))
-              jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2))
-
-              nxhr = iimax - iimin + 1
-              nyhr = jjmax - jjmin + 1                    
-
-              IF( nxhr == 0 .OR. nyhr == 0 ) THEN
-                 !
-                 trouble_points(ji,jj) = 1
-                 !
-              ELSE
-                 !
-                 ALLOCATE( vardep(nxhr,nyhr), mask_oce(nxhr,nyhr) )
-                 vardep(:,:) = G0%bathy_meter(iimin:iimax,jjmin:jjmax)
-                 !
-                 WHERE( vardep(:,:) .GT. 0. )   ;   mask_oce = 1 ;
-                 ELSEWHERE                      ;   mask_oce = 0 ;
-                 ENDWHERE
-                 !
-                 nxyhr = nxhr*nyhr
-                 IF( SUM(mask_oce) < 0.5*(nxyhr) ) THEN ! if more than half of the points are on land then bathy fine = 0
-                    G1%bathy_meter(ji,jj) = 0.
-                 ELSE
-                    IF( type_bathy_interp == 0 ) THEN     ! Arithmetic average
-                       G1%bathy_meter(ji,jj) = SUM( vardep(:,:) * mask_oce(:,:) ) / SUM( mask_oce(:,:) )
-                    ELSE                                  ! Median average
-                       ALLOCATE(vardep1d(nxyhr))
-                       vardep1d = RESHAPE(vardep,(/ nxyhr /) )
-                       !!CALL ssort(vardep1d,nxyhr)
-                       CALL quicksort(vardep1d,1,nxyhr)
-                       !
-                       ! Calculate median
-                       IF (MOD(nxyhr,2) .NE. 0) THEN
-                          G1%bathy_meter(ji,jj) = vardep1d( nxyhr/2 + 1 )
-                       ELSE
-                          G1%bathy_meter(ji,jj) = 0.5 * ( vardep1d(nxyhr/2) + vardep1d(nxyhr/2+1) )
-                       END IF
-                       DEALLOCATE(vardep1d)   
-                    ENDIF
-                 ENDIF
-                 DEALLOCATE (mask_oce,vardep)
-                 !
-              ENDIF
            ENDDO
-        ENDDO
-
-        IF( SUM( trouble_points ) > 0 ) THEN
-           PRINT*,'too much empty cells, proceed to bilinear interpolation'
-           type_bathy_interp = 2
-        ENDIF
-           
-     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.                          
+
+           IF( SUM( trouble_points ) > 0 ) THEN
+              PRINT*,'too much empty cells, proceed to bilinear interpolation'
+              type_bathy_interp = 2
            ENDIF
-        ENDIF
-
-        IF( .NOT. identical_grids ) THEN 
+
+        ENDIF
+        !                                                       ! ----------------------------- 
+        IF( type_bathy_interp == 2) THEN                        ! Bilinear interpolation
+           !                                                    ! ----------------------------- 
 
            ALLOCATE(masksrc(SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2)))
@@ -400,5 +403,5 @@
            !  masksrc = .false.
            !ElseWhere
-              masksrc = .TRUE.
+           masksrc = .TRUE.
            !End where                       
            !            
@@ -411,8 +414,8 @@
            DEALLOCATE(masksrc)
            DEALLOCATE(bathy_test) 
-
+           
         ENDIF
         !            
-     ENDIF
+     ENDIF ! not identical grids
      ! ---
      ! At this stage bathymetry in meters has already been interpolated on fine grid
