Index: /utils/tools/NESTING/src/agrif_create_bathy.f90
===================================================================
--- /utils/tools/NESTING/src/agrif_create_bathy.f90	(revision 12258)
+++ /utils/tools/NESTING/src/agrif_create_bathy.f90	(revision 12259)
@@ -84,4 +84,59 @@
      STOP
   ENDIF
+  !
+  ! 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)
+     CALL meter_to_levels(G0)
+  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)
+        CALL meter_to_levels(G0)
+     ENDIF
+     ! change longitudes (from -180:180 to 0:360)
+     IF(Pacifique) THEN
+        WHERE(G0%nav_lon < 0.001)   G0%nav_lon = G0%nav_lon + 360.
+     ENDIF
+  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 .AND. .NOT.partial_steps ) THEN   ! **** First option : no new topo file & no partial steps
@@ -90,53 +145,22 @@
      !                                                ! ------------------------------------------------------------------
      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
-        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
-
-     ! change longitudes (from -180:180 to 0:360)
-     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)
-     !      
+     !
+     IF( .NOT. identical_grids ) THEN
+        !           
+        WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Bilinear interpolation ...'
+        !
+        ! create logical array masksrc
+        ALLOCATE(masksrc(SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2)))
+        WHERE(G0%bathy_level.LE.0)   ;   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,G1%bathy_meter,nxfin,nyfin,matrix,src_add,dst_add)  
+        DEALLOCATE(masksrc)
+        !
+     ENDIF
+     
      ! compute constant bathymetry for Parent-Child bathymetry connection
      CALL init_constant_bathy(G0%bathy_meter,bathy_fin_constant)
@@ -210,7 +234,15 @@
      CALL levels_to_meter(G1) ! needed for domcfg
      !
+     ! 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
-     status = Write_Bathy_level(TRIM(child_level),G1)
-     status = write_domcfg(TRIM(child_domcfg),G1)
+     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)
      !!status = write_domcfg(TRIM(parent_domcfg_out),G0) ! do not activate it
 
@@ -222,39 +254,7 @@
      !                                                ! ----------------------------------------------------- 
      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
-     !               
+     !
+     ! === Here: G0 is the grid associated with the new topography (as gebco or etopo) ===
+     !
      ! what type of interpolation for bathymetry
      IF( type_bathy_interp == 0 ) THEN
@@ -269,7 +269,4 @@
      ENDIF
      !
-     ! 1st allocation of child grid bathy
-     ALLOCATE(G1%bathy_meter(nxfin,nyfin))
-     G1%bathy_meter(:,:)=0.                       
      !
      ! ---------------------------------------------------------------------------------
@@ -278,17 +275,4 @@
      ! ==> It gives G1%bathy_meter from G0%bathy_meter
      ! ---------------------------------------------------------------------------------
-     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
@@ -392,4 +376,6 @@
            ENDIF
 
+           DEALLOCATE(trouble_points)
+
         ENDIF
         !                                                       ! ----------------------------- 
@@ -400,9 +386,7 @@
            ALLOCATE(bathy_test(nxfin,nyfin))
            !
-           !Where(G0%bathy_meter.le.0.00001) 
-           !  masksrc = .false.
-           !ElseWhere
-           masksrc = .TRUE.
-           !End where                       
+           WHERE(G0%bathy_level.LE.0)   ;   masksrc = .FALSE.   ;
+           ELSEWHERE                    ;   masksrc = .TRUE.    ;
+           END WHERE
            !            
            ! compute remapping matrix thanks to SCRIP package            
@@ -446,5 +430,5 @@
      ! 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(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)))
      !                       
@@ -619,9 +603,10 @@
      !
      ! 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)
-     
+     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)
