Index: /utils/tools/NESTING/namelist_nordic1
===================================================================
--- /utils/tools/NESTING/namelist_nordic1	(revision 10024)
+++ /utils/tools/NESTING/namelist_nordic1	(revision 10025)
@@ -1,4 +1,4 @@
 &input_output
-    iom_activated = true 
+    iom_activated = .true. 
 /
 
@@ -6,15 +6,17 @@
     parent_coordinate_file = 'coordinates.nc'
     parent_meshmask_file   = 'meshmask.nc'  
+    parent_domcfg_output   = 'domain_cfg.nc'
 /    
     
 &bathymetry    
-    new_topo = true
+    new_topo = .true.
     elevation_database = 'GEBCO_2014_2D.nc'
     elevation_name     = 'elevation'
-    smoothing          = true
+    smoothing          = .true.
     smoothing_factor   = 0.6
+    ln_agrif_domain    = .true.
     npt_connect        = 2      ! default = 3
     npt_copy           = 2      ! default = 2
-    removeclosedseas   = true
+    removeclosedseas   = .true.
     type_bathy_interp  = 0 
     rn_hmin            = -3 
@@ -28,6 +30,8 @@
     rho  = 4
     rhot = 4
-    bathy_update        = true
+    nbghostcellsfine = 3
+    bathy_update = .true.
     updated_parent_file = 'bathy_updated.nc'				
+    updated_parent_domcfg = 'domain_cfg_updated.nc'
 /
 
@@ -42,5 +46,5 @@
     pa1       = 245.58132232490
     N         = 31
-    ldbletanh = .FALSE.
+    ldbletanh = .false.
     pa2       = 0
     ppkth2    = 0
@@ -49,5 +53,5 @@
     
 &partial_cells    	
-    partial_steps      = true
+    partial_steps      = .true.
     parent_bathy_meter = 'bathy_meter.nc'
     parent_batmet_name = 'Bathymetry'
Index: /utils/tools/NESTING/namelist_nordic1_update
===================================================================
--- /utils/tools/NESTING/namelist_nordic1_update	(revision 10024)
+++ /utils/tools/NESTING/namelist_nordic1_update	(revision 10025)
@@ -1,4 +1,4 @@
 &input_output
-    iom_activated = true 
+    iom_activated = .true. 
 /
 
@@ -6,4 +6,5 @@
     parent_coordinate_file = 'coordinates.nc'
     parent_meshmask_file   = 'meshmask.nc'  
+    parent_domcfg_output   = 'domain_cfg.nc'
 /    
     
@@ -14,7 +15,8 @@
     smoothing          = true
     smoothing_factor   = 0.6
+    ln_agrif_domain    = .true.
     npt_connect        = 2    
     npt_copy           = 2 
-    removeclosedseas   = true
+    removeclosedseas   = .true.
     type_bathy_interp  = 2 
     rn_hmin            = -3 
@@ -28,6 +30,8 @@
     rho  = 4
     rhot = 4
+    nbghostcellsfine = 3
     bathy_update = true
     updated_parent_file = 'bathy_updated_parent.nc'
+    updated_parent_domcfg = 'domain_cfg_updated_parent.nc'
 /
 
@@ -42,5 +46,5 @@
     pa1       = 245.58132232490
     N         = 31
-    ldbletanh = .FALSE.
+    ldbletanh = .false.
     pa2       = 0
     ppkth2    = 0
@@ -49,5 +53,5 @@
     
 &partial_cells    	
-    partial_steps      = true
+    partial_steps      = .true.
     parent_bathy_meter = 'bathy_meter.nc'
     parent_batmet_name = 'Bathymetry'
Index: /utils/tools/NESTING/namelist_nordic2
===================================================================
--- /utils/tools/NESTING/namelist_nordic2	(revision 10024)
+++ /utils/tools/NESTING/namelist_nordic2	(revision 10025)
@@ -1,4 +1,4 @@
 &input_output
-    iom_activated = true
+    iom_activated = .true.
 /
 
@@ -6,15 +6,17 @@
     parent_coordinate_file = '1_coordinates.nc'
     parent_meshmask_file   = '1_meshmask.nc'
+    parent_domcfg_output   = '1_domain_cfg.nc'
 /    
     
 &bathymetry    
-    new_topo           = true
+    new_topo           = .true.
     elevation_database = 'GEBCO_2014_2D.nc'
     elevation_name     = 'elevation'
-    smoothing          = true
+    smoothing          = .true.
     smoothing_factor   = 0.6
+    ln_agrif_domain    = .true.
     npt_connect        = 2 
     npt_copy           = 2 
-    removeclosedseas   = true
+    removeclosedseas   = .true.
     type_bathy_interp  = 0 
     rn_hmin            = -3
@@ -28,6 +30,8 @@
     rho  = 3
     rhot = 3
-    bathy_update = true
+    nbghostcellsfine = 3
+    bathy_update = .true.
     updated_parent_file = '1_bathy_updated.nc'
+    updated_parent_domcfg = '1_domain_cfg_updated.nc'
 /
 
@@ -42,5 +46,5 @@
     pa1       = 245.58132232490
     N         = 31
-    ldbletanh = .FALSE.
+    ldbletanh = .false.
     pa2       = 0
     ppkth2    = 0
@@ -49,5 +53,5 @@
     
 &partial_cells    	
-    partial_steps      = true
+    partial_steps      = .true.
     parent_bathy_meter = '1_bathy_meter.nc'
     parent_batmet_name = 'Bathymetry'
Index: /utils/tools/NESTING/src/agrif_connect_topo.f90
===================================================================
--- /utils/tools/NESTING/src/agrif_connect_topo.f90	(revision 10024)
+++ /utils/tools/NESTING/src/agrif_connect_topo.f90	(revision 10025)
@@ -51,38 +51,56 @@
     TYPE(Coordinates) :: Grid     
     !
-    !
-    diff = 0
-    IF(MOD(rho,2) .EQ. 0) diff = 1
-    !
-    indx = 1 + nbghostcellsfine + CEILING(irafx/2.0) + diff      
-    indy = 1 + nbghostcellsfine + CEILING(irafy/2.0) + diff
-    bornex = 1+nbghostcellsfine + CEILING(irafx/2.0) + diff - irafx
-    borney = 1+nbghostcellsfine + CEILING(irafy/2.0) + diff - irafy
-    bornex2 = nxfin - (nbghostcellsfine) + irafx - CEILING(irafx/2.0) 
-    borney2 = nyfin - (nbghostcellsfine) + irafy - CEILING(irafy/2.0) 
-    !    
-    ALLOCATE(bathy_fin_constant(bornex-FLOOR(irafx/2.0):bornex2+FLOOR(irafx/2.0), &
-         borney-FLOOR(irafy/2.0):borney2+FLOOR(irafy/2.0)))
-    !
-    DO j = borney,borney2,irafy
-
-       jpt = jmin + 1 + nbghostcellsfine + (j-indy)/irafy
-       IF(j<=1) jpt = jmin + 1
-
-       DO i = bornex,bornex2,irafx
-
-          ipt = imin + 1 + nbghostcellsfine + (i-indx)/irafx
-          IF(i<=1) ipt = imin + 1
-          !       
-          DO jj = j-FLOOR(irafy/2.0),j+FLOOR(irafy/2.0)-diff
-             DO ii = i-FLOOR(irafx/2.0),i+FLOOR(irafx/2.0)-diff
-
-                bathy_fin_constant(ii,jj) = coarse_bathy(ipt,jpt)
-
+    IF( ln_agrif_domain ) THEN
+    !
+       diff = 0
+       IF(MOD(rho,2) .EQ. 0) diff = 1
+       !
+       indx = 1 + nbghostcellsfine + CEILING(irafx/2.0) + diff      
+       indy = 1 + nbghostcellsfine + CEILING(irafy/2.0) + diff
+       bornex = 1+nbghostcellsfine + CEILING(irafx/2.0) + diff - irafx
+       borney = 1+nbghostcellsfine + CEILING(irafy/2.0) + diff - irafy
+       bornex2 = nxfin - (nbghostcellsfine) + irafx - CEILING(irafx/2.0) 
+       borney2 = nyfin - (nbghostcellsfine) + irafy - CEILING(irafy/2.0) 
+       !    
+       ALLOCATE(bathy_fin_constant(bornex-FLOOR(irafx/2.0):bornex2+FLOOR(irafx/2.0), &
+          borney-FLOOR(irafy/2.0):borney2+FLOOR(irafy/2.0)))
+       !
+       DO j = borney,borney2,irafy
+          
+          jpt = jmin + 1 + nbghostcellsfine + (j-indy)/irafy
+          IF(j<=1) jpt = jmin + 1
+          
+          DO i = bornex,bornex2,irafx
+             
+             ipt = imin + 1 + nbghostcellsfine + (i-indx)/irafx
+             IF(i<=1) ipt = imin + 1
+             !       
+             DO jj = j-FLOOR(irafy/2.0),j+FLOOR(irafy/2.0)-diff
+                DO ii = i-FLOOR(irafx/2.0),i+FLOOR(irafx/2.0)-diff
+                   
+                   bathy_fin_constant(ii,jj) = coarse_bathy(ipt,jpt)
+                   
+                END DO
              END DO
+             
           END DO
-
-       END DO
-    END DO
+       END DO
+       
+    ELSE
+
+       ALLOCATE(bathy_fin_constant(1:nxfin,1:nyfin))
+
+       DO j = 1,nyfin-irafy+1,irafy
+          jpt = jmin + FLOOR( (j - 1.) / irafy )
+          !
+          DO i = 1,nxfin-irafx+1,irafx
+             ipt = imin + FLOOR( (i - 1.) / irafx )
+             !
+             bathy_fin_constant(i:i+irafx-1,j:j+irafy-1) = coarse_bathy(ipt,jpt)
+             !
+          END DO
+       END DO
+ 
+    ENDIF
     !
     !
Index: /utils/tools/NESTING/src/agrif_create_bathy.f90
===================================================================
--- /utils/tools/NESTING/src/agrif_create_bathy.f90	(revision 10024)
+++ /utils/tools/NESTING/src/agrif_create_bathy.f90	(revision 10025)
@@ -13,4 +13,5 @@
   USE agrif_partial_steps
   USE agrif_connect_topo
+  USE agrif_interpolation
   !
   IMPLICIT NONE
@@ -40,21 +41,21 @@
   !      
   CHARACTER(len=80) :: namelistname
-  CHARACTER*100 :: Childmeter_file,Childlevel_file,Child_coordinates,child_ps     
-  LOGICAL,DIMENSION(:,:),POINTER :: masksrc => NULL()  
+  CHARACTER*100     :: Childmeter_file,Childlevel_file,Child_coordinates,child_ps, Child_domcfg    
   LOGICAL :: identical_grids     
-  INTEGER,DIMENSION(:,:),ALLOCATABLE ::mask_oce,trouble_points
   INTEGER :: i,j,num_links,nb,nbadd,status,narg,iargc     
-  INTEGER,DIMENSION(:),POINTER :: src_add,dst_add => NULL() 
-  INTEGER :: numlatfine,numlonfine,numlat,numlon,pos,pos2
-  REAL*8,DIMENSION(:,:),POINTER :: matrix,interpdata => NULL()     
-  REAL*8, DIMENSION(:,:),POINTER :: bathy_fin_constant => NULL()  
-  REAL*8,DIMENSION(:,:),ALLOCATABLE :: bathy_test,vardep,glamhr,gphihr
-  REAL*8,DIMENSION(:),ALLOCATABLE :: vardep1d
-  REAL*8, DIMENSION(:,:),POINTER :: gdepw_ps_interp => NULL() 
-  REAL*8, DIMENSION(:,:),POINTER :: save_gdepw,rx,ry,maskedtopo
+  INTEGER :: jpj,jpi,pos,pos2
+  LOGICAL,DIMENSION(:,:),POINTER     :: masksrc => NULL()  
+  INTEGER,DIMENSION(:,:),ALLOCATABLE :: mask_oce,trouble_points
+  INTEGER,DIMENSION(:)  ,POINTER     :: src_add,dst_add => NULL() 
+  REAL*8, DIMENSION(:,:),POINTER     :: matrix,interpdata => NULL()     
+  REAL*8, DIMENSION(:,:),POINTER     :: bathy_fin_constant => NULL()  
+  REAL*8, DIMENSION(:,:),ALLOCATABLE :: bathy_test,vardep,glamhr,gphihr
+  REAL*8, DIMENSION(:)  ,ALLOCATABLE :: vardep1d
+  REAL*8, DIMENSION(:,:),POINTER     :: gdepw_ps_interp => NULL() 
+  REAL*8, DIMENSION(:,:),POINTER     :: save_gdepw,rx,ry,maskedtopo
   REAL*8  :: Cell_lonmin,Cell_lonmax,Cell_latmin,Cell_latmax,wghts
-  LOGICAL :: Pacifique=.FALSE.
+  LOGICAL :: Pacifique = .FALSE.
   INTEGER :: boundary,xpos,ypos,iimin,iimax,jjmax,jjmin
-  INTEGER :: nbloops,nxhr,nyhr,ji,jj,nbiter,nbloopmax
+  INTEGER :: nxhr,nyhr,ji,jj,nbiter
   INTEGER :: ipt,jpt,iloc,jloc
   INTEGER, DIMENSION(2) :: i_min,i_max,j_min,j_max
@@ -70,216 +71,201 @@
   !
   ! read input file (namelist.input)
-  !
   CALL read_namelist(namelistname)
   !      
   ! define names of child grid files
+  CALL set_child_name(parent_coordinate_file,Child_coordinates) 
+  IF( TRIM(parent_meshmask_file) .NE. '/NULL' )    CALL set_child_name(parent_meshmask_file,Childlevel_file)            
+  IF( TRIM(parent_meshmask_file) .NE. '/NULL' )    CALL set_child_name(parent_domcfg_output,Child_domcfg) 
   !
-  CALL set_child_name(parent_coordinate_file,Child_coordinates) 
-  IF( TRIM(parent_meshmask_file) .NE. '/NULL' ) &
-       CALL set_child_name(parent_meshmask_file,Childlevel_file)            
-  !
-  !
-  !
-  !
-  !
-  !------------------------------------------------------------------
-  ! ****First option : no new topo file / no partial steps
-  ! interpolate levels directly from parent file
-  !------------------------------------------------------------------
-  !
-  !
-  !
-  !
-  !
-  !
-  IF(.NOT.new_topo .AND. .NOT.partial_steps) THEN      
+  !                                                   ! ------------------------------------------------------------------
+  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'
      !      
-     WRITE(*,*) 'First option'
-     !read coarse grid bathymetry and coordinates file
-     !           
-     WRITE(*,*) 'No new topo file ...'
+     ! read coarse grid bathymetry and coordinates
      status = Read_Coordinates(TRIM(parent_coordinate_file),G0)    
      status = Read_bathy_level(TRIM(parent_meshmask_file),G0)
      !           
-     IF( imax > SIZE(G0%glamt,1) .OR. jmax > SIZE(G0%glamt,2) .OR. &
-          imax <= imin .OR. jmax <= jmin ) THEN                   
+     ! 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'       
-        WRITE(*,*) ' '
         STOP
      ENDIF
-     !
-     !read fine grid coordinates file
-     !      
-     status = Read_Coordinates(TRIM(Child_coordinates),G1,pacifique)
-     !  
      IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN
-        !
         WRITE(*,*) 'ERROR ***** bad child coordinates file ...'
-        WRITE(*,*) ' '
-        WRITE(*,*) 'please check that child coordinates file '
-        WRITE(*,*) 'has been created with the same namelist '
-        WRITE(*,*) ' '
+        WRITE(*,*) 'please check that child coordinates file has been created with the same namelist'
         STOP
-        !
-     ENDIF
-     !
-     numlat =  SIZE(G0%nav_lat,2)
-     numlon =  SIZE(G0%nav_lat,1)    
-     numlatfine =  SIZE(G1%nav_lat,2)
-     numlonfine =  SIZE(G1%nav_lat,1)  
+     ENDIF
+     !
+     jpj = SIZE(G0%nav_lat,2)
+     jpi = SIZE(G0%nav_lat,1)    
      !           
-     ALLOCATE(masksrc(numlon,numlat))
-     !
      ! create logical array masksrc
-     !
-     WHERE(G0%bathy_level.LE.0) 
-        masksrc = .FALSE.
-     ELSEWHERE
-        masksrc = .TRUE.
+     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.
-        END WHERE
-     ENDIF
-
-
-     !-----------------          
-     ! compute remapping matrix thanks to SCRIP package
-     !
-     ! remapping process
-     !              
+        WHERE(G0%nav_lon < 0.001)   G0%nav_lon = G0%nav_lon + 360.
+     ENDIF
+
+     ! From levels to meters
      ALLOCATE(G1%bathy_meter(nxfin,nyfin))
      CALL levels_to_meter(G0)
      !             
-     !             Call levels_to_meter(G1)
-     !             
-     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)  
-     !                                  
-     !            
+     ! 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)
      !
-     boundary = npt_copy*irafx + nbghostcellsfine + 1 
-     !
-     ! connection carried out by copying parent grid values for the fine points
-     ! corresponding to 3 coarse grid cells at the boundaries
-     !                 
+     ! 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)
      !                  
-     CALL smooth_topo(G1%bathy_meter(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter)
-     !             
+     ! 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)
-     !             
-     DEALLOCATE(bathy_fin_constant)
-     !
-     !
-     !------------------------------------------------------------------
-     ! ****Second option : new topo file or/and partial steps     
-     !------------------------------------------------------------------ 
-     !
-     !
-     !
-     !
-     !
-     !
-     !
-     !
-  ELSE
-     !
-     WRITE(*,*) 'Second option : partial steps'
+     G1%bathy_level=NINT(G1%bathy_level)
+     !	    
+     ! Check errors in bathy
+     DO j=1,nyfin
+        DO i=1,nxfin
+           IF (G1%bathy_level(i,j).LT.0.) THEN
+              PRINT *,'error in ',i,j,G1%bathy_level(i,j)
+              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 j=2,nyfin-1
+              DO i=2,nxfin-1
+                 IF (G1%bathy_level(i,j).GT.0.) THEN
+                    IF (MAX(bathy_test(i,j+1),bathy_test(i,j-1),bathy_test(i-1,j),bathy_test(i+1,j)).EQ.1) THEN
+                       IF (bathy_test(i,j).NE.1.) nbadd = nbadd + 1
+                       bathy_test(i,j)=1.
+                    ENDIF
+                 ENDIF
+              ENDDO
+           ENDDO
+        ENDDO
+
+        WHERE (bathy_test.EQ.0.)   G1%bathy_level = 0.
+
+        DEALLOCATE(bathy_test)           
+     ENDIF
+     !
+     ! store interpolation result in output file
+     CALL levels_to_meter(G1)   ! From levels to meters
+     status = Write_Bathy_level(TRIM(Childlevel_file),G1)
+     status = write_domcfg(TRIM(Child_domcfg),G1)
+
+     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)
-     !                        
-     IF( imax > SIZE(G0%nav_lon,1) .OR. jmax > SIZE(G0%nav_lon,2) .OR. &
-          imax <= imin .OR. jmax <= jmin ) THEN                    
+     !
+     ! 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'       
-        WRITE(*,*) ' '
         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(*,*) ' '
-        WRITE(*,*) 'please check that child coordinates file '
-        WRITE(*,*) 'has been created with the same namelist '
-        WRITE(*,*) ' '
+        WRITE(*,*) 'please check that child coordinates file has been created with the same namelist'
         STOP
-        !
-     ENDIF
-     !      
-     ! read coarse grid bathymetry file
-     !---
-     IF( new_topo ) THEN
-        WRITE(*,*) 'use new topo file : ',TRIM(elevation_database)
+     ENDIF
+     !
+     ! === From here on: G0 is the grid associated with the new topography (like 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_meter(TRIM(elevation_database),G0,G1,Pacifique)
-     ELSE
-        WRITE(*,*) 'no new topo file'
-        status = Read_Bathymeter(TRIM(parent_bathy_meter),G0)
+     ELSE                                                      ! read G0%bathy_meter (on the global grid)
+        status = Read_Bathymeter(TRIM(parent_bathy_meter),G0)   
         IF(Pacifique) THEN
-           WHERE(G0%nav_lon < 0.0001) 
-              G0%nav_lon = G0%nav_lon + 360.
-           END WHERE
-        ENDIF
-     ENDIF
-     !---            
-     numlatfine =  SIZE(G1%nav_lat,2)
-     numlonfine =  SIZE(G1%nav_lat,1) 
-     !  
+           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.                       
-
-     WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid'
-
-     IF( type_bathy_interp == 0 ) THEN
-        WRITE(*,*) 'Arithmetic average ...'
-     ELSE IF( type_bathy_interp == 1 ) THEN
-        WRITE(*,*) 'Median average ...'
-     ELSE IF( type_bathy_interp == 2 ) THEN     
-        WRITE(*,*) 'Bilinear interpolation ...'
-     ELSE     
-        WRITE(*,*) 'bad value for type_bathy_interp variable ( must be 0,1 or 2 )'
-        STOP 
-     ENDIF
-     !
-     !************************************
-     !MEDIAN AVERAGE or ARITHMETIC AVERAGE
-     !************************************
-     !
-     IF( type_bathy_interp == 0 .OR. type_bathy_interp == 1 ) THEN 
-        !
+     !
+     ! ---------------------------------------------------------------------------------
+     ! ===                 Bathymetry of the fine grid (step1)                       ===
+     ! ---------------------------------------------------------------------------------
+     ! ==> 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
-        !
-        !  POINT DETECTION
-        !
+        trouble_points(:,:) = 0
         !                       
-        DO jj = 2,numlatfine
-           DO ji = 2,numlonfine
+        DO jj = 2,nyfin
+           DO ji = 2,nxfin
               !	    
-              ! FINE GRID CELLS DEFINITION               
-              !
+              ! 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))
@@ -287,6 +273,5 @@
               Cell_latmax = MAX(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj)) 
               !               
-              ! SEARCH FOR ALL POINTS CONTAINED IN THIS CELL
-              !
+              ! look for points in G0 (bathy dataset) contained in the fine grid cells  
               iimin = 1
               DO WHILE( G0%nav_lon(iimin,1) < Cell_lonmin ) 
@@ -302,5 +287,5 @@
               DO WHILE( G0%nav_lon(iimax,1) <= Cell_lonmax ) 
                  iimax = iimax + 1
-              iimax = MIN( iimax,SIZE(G0%bathy_meter,1))
+                 iimax = MIN( iimax,SIZE(G0%bathy_meter,1))
               ENDDO
               !                               
@@ -308,10 +293,14 @@
               DO WHILE( G0%nav_lat(iimax,jjmax) <= Cell_latmax ) 
                  jjmax = jjmax + 1
-              jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2))
-
+                 jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2))
               ENDDO
               !
-              iimax = iimax-1
-              jjmax = jjmax-1
+              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 )
@@ -324,26 +313,23 @@
 
               IF( nxhr == 0 .OR. nyhr == 0 ) THEN
+                 !
                  trouble_points(ji,jj) = 1
+                 !
               ELSE
-
-                 ALLOCATE( vardep(nxhr,nyhr) )
-                 ALLOCATE( mask_oce(nxhr,nyhr) )
-                 mask_oce = 0       
-
+                 !
+                 ALLOCATE( vardep(nxhr,nyhr), mask_oce(nxhr,nyhr) )
                  vardep(:,:) = G0%bathy_meter(iimin:iimax,jjmin:jjmax)
-
-                 WHERE( vardep(:,:) .GT. 0. )  mask_oce = 1
-
-!                 IF( SUM(mask_oce) == 0 ) THEN
-                 IF( SUM(mask_oce) < 0.5*(nxhr*nyhr) ) THEN
+                 !
+                 WHERE( vardep(:,:) .GT. 0. )   ;   mask_oce = 1 ;
+                 ELSEWHERE                      ;   mask_oce = 0 ;
+                 ENDWHERE
+                 !
+                 IF( SUM(mask_oce) < 0.5*(nxhr*nyhr) ) 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		    
-                       !
-                       vardep(:,:) = vardep(:,:)*mask_oce(:,:) - 100000*(1-mask_oce(:,:))
+                    IF( type_bathy_interp == 0 ) THEN     ! Arithmetic average
+                       G1%bathy_meter(ji,jj) = SUM( vardep(:,:) * mask_oce(:,:) ) / SUM( mask_oce(:,:) )
+                    ELSE                                  ! Median average
+                       vardep(:,:) = vardep(:,:) * mask_oce(:,:) - 100000 * ( 1 - mask_oce(:,:) )
                        ALLOCATE(vardep1d(nxhr*nyhr))
                        vardep1d = RESHAPE(vardep,(/ nxhr*nyhr /) )
@@ -351,9 +337,8 @@
                        !
                        ! Calculate median
-                       !
                        IF (MOD(SUM(mask_oce),2) .NE. 0) THEN
-                          G1%bathy_meter(ji,jj) = vardep1d( SUM(mask_oce)/2 + 1)
+                          G1%bathy_meter(ji,jj) = vardep1d( SUM(mask_oce)/2 + 1 )
                        ELSE
-                          G1%bathy_meter(ji,jj) = ( vardep1d(SUM(mask_oce)/2) + vardep1d(SUM(mask_oce)/2+1) )/2.0
+                          G1%bathy_meter(ji,jj) = 0.5 * ( vardep1d(SUM(mask_oce)/2) + vardep1d(SUM(mask_oce)/2+1) )
                        END IF
                        DEALLOCATE(vardep1d)   
@@ -361,5 +346,5 @@
                  ENDIF
                  DEALLOCATE (mask_oce,vardep)
-
+                 !
               ENDIF
            ENDDO
@@ -367,31 +352,23 @@
 
         IF( SUM( trouble_points ) > 0 ) THEN
-           PRINT*,'too much empty cells, proceed to bilinear interpolation !!'
+           PRINT*,'too much empty cells, proceed to bilinear interpolation'
            type_bathy_interp = 2
         ENDIF
-
-     ENDIF
-
-     !
-     ! create logical array masksrc
-     !
-     IF( type_bathy_interp == 2) THEN 
-        !
-
-        !            
+           
+     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( 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
+               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 
-              PRINT*,'same grid between ',elevation_database,' and child domain'    
               identical_grids = .TRUE.                          
            ENDIF
         ENDIF
-
 
         IF( .NOT. identical_grids ) THEN 
@@ -400,19 +377,13 @@
            ALLOCATE(bathy_test(nxfin,nyfin))
            !
-           !                    Where(G0%bathy_meter.le.0.00001) 
-           !	                masksrc = .false.
-           !	            ElseWhere
-           !
-           masksrc = .TRUE.
-           !
-           !	            End where                       
+           !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)  
+           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               
@@ -424,142 +395,146 @@
         !            
      ENDIF
-     !
-     !
-     !
-     !------------------------------------------------------------------------------------------
-     ! ! ****Third  option : Partial Steps
-     ! The code includes the 
-     ! option to include partial cells at the bottom 
-     ! in order to better resolve topographic variations
-     !------------------------------------------------------------------------------------------
-     !
-     ! At this step bathymetry in meters has already been interpolated on fine grid
-     !
-     ! 	                
-     IF( partial_steps ) THEN                
-        !                  
+     ! ---
+     ! 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( partial_steps ) THEN                     
         status = Read_Bathymeter(TRIM(parent_bathy_meter),G0)
-        DEALLOCATE(G0%nav_lat,G0%nav_lon)
-        status = Read_coordinates(TRIM(parent_coordinate_file),G0)
-        !------------------------                  
-        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)))
-        !
-        !                       
-        WRITE(*,*) 'Coarse grid : '
-        CALL get_partial_steps(G0) 
-        WRITE(*,*) ' '
-        WRITE(*,*) 'Fine grid : '
-        CALL get_partial_steps(G1)                 ! compute gdepw_ps for G1
-        CALL bathymetry_control(G0%Bathy_level)    !    
-        CALL Check_interp(G0,gdepw_ps_interp)      ! interpolation in connection zone (3 coarse grid cells)
-        !
-        boundary = npt_copy*irafx + nbghostcellsfine + 1                     
-! J chanut: exclude matching if no open boundaries
-        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
-
-        DO jj=1,nyfin
-          ! West
-          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 
-          ! East
-          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
+     ELSE
+        status = Read_Bathymeter(TRIM(parent_meshmask_file),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
+     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
+     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
         END DO
-        DO ji=1,nxfin
-          ! South 
-          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
-          ! North
-          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
-        
-!        G1%gdepw_ps(1:boundary,:) = gdepw_ps_interp(1:boundary,:)
-!        G1%gdepw_ps(:,1:boundary) = gdepw_ps_interp(:,1:boundary)
-!        G1%gdepw_ps(nxfin-boundary+1:nxfin,:) = gdepw_ps_interp(nxfin-boundary+1:nxfin,:)
-!        G1%gdepw_ps(:,nyfin-boundary+1:nyfin) = gdepw_ps_interp(:,nyfin-boundary+1:nyfin)
-        !                   
-
-        IF(.NOT. smoothing) WRITE(*,*) 'No smoothing process only connection is carried out'
-        WRITE(*,*) ' linear connection on ',npt_connect,'coarse grid points'
-
-        !            
-        gdepw_ps_interp = 0.
-        CALL Check_interp(G0,gdepw_ps_interp)      ! interpolation in connection zone (3 coarse grid cells)
-        !
-        !
-        !
-        !
-        !                    LINEAR CONNECTION
-        !
-        !
-        !
-        !
-        !
+
         wghts = 1.
-        DO ji = boundary + 1 , boundary + npt_connect * irafx
-            wghts = wghts - (1. / (npt_connect*irafx - 1. ) )
-            DO jj=1,nyfin
-               IF (G1%gdepw_ps(nbghostcellsfine+1,jj) > 0.) THEN
-                  G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))  
-               ENDIF
-            END DO
-        END DO 
-      
+        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 ji = nxfin - boundary , nxfin - boundary - npt_connect * irafx +1 , -1
-            wghts = wghts - (1. / (npt_connect*irafx - 1. ) )
-            DO jj=1,nyfin
-               IF (G1%gdepw_ps(nxfin-nbghostcellsfine,jj) > 0.) THEN
-                  G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))
-               ENDIF
-            END DO
-        END DO  
-
-        wghts = 1.
-        DO jj = boundary + 1 , boundary + npt_connect * irafy
-            wghts = wghts - (1. / (npt_connect*irafy - 1. ) )
-            DO ji=1,nxfin
-               IF (G1%gdepw_ps(ji,nbghostcellsfine+1) > 0.) THEN
-                  G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))
-               ENDIF
-            END DO
-        END DO
-
-        wghts = 1.
-        DO jj = nyfin - boundary , nyfin - boundary - npt_connect * irafy +1, -1
-            wghts = wghts - (1. / (npt_connect*irafy - 1. ) )
-            DO ji=1,nxfin
-               IF (G1%gdepw_ps(ji,nyfin-nbghostcellsfine) > 0.) THEN
-                  G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))
-               ENDIF
-            END DO
+        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
@@ -567,72 +542,62 @@
         ENDIF
 
-        G1%bathy_meter = G1%gdepw_ps
-        !                     
-        ! 
-! Chanut: remove smoothing if child grid bathymetry is found to already exist
-        IF(smoothing.AND.(.NOT.identical_grids)) THEN 
-
-           !
-           ! Smoothing to connect the connection zone (3 + npt_connect coarse grid cells) and the interior domain
-           !
-! Chanut: smoothing everywhere then discard result in connection zone
-           CALL smooth_topo(G1%gdepw_ps(1:nxfin,1:nyfin),nbiter)
-           WHERE (G1%wgt(:,:)==0.) G1%bathy_meter(:,:) = G1%gdepw_ps(:,:)
-!           boundary = (npt_copy+npt_connect)*irafx + nbghostcellsfine + 1 
-!           CALL smooth_topo(G1%gdepw_ps(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter)
-!           G1%bathy_meter = G1%gdepw_ps                         
-        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
-           END WHERE
-           WHERE (G1%bathy_meter(nxfin,:).GT.0.)
-              bathy_test(nxfin,:)=1
-           END WHERE
-           WHERE (G1%bathy_meter(:,1).GT.0.)
-              bathy_test(:,1)=1
-           END WHERE
-           WHERE (G1%bathy_meter(:,nyfin).GT.0.)
-              bathy_test(:,nyfin)=1
-           END WHERE
-           nbadd = 1
-           DO WHILE (nbadd.NE.0)
-              nbadd = 0
-              DO j=2,nyfin-1
-                 DO i=2,nxfin-1
-                    IF (G1%bathy_meter(i,j).GT.0.) THEN
-                       IF (MAX(bathy_test(i,j+1),bathy_test(i,j-1), &
-                            bathy_test(i-1,j),bathy_test(i+1,j)).EQ.1) THEN
-                          IF (bathy_test(i,j).NE.1.) nbadd = nbadd + 1
-                          bathy_test(i,j)=1.
-                       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 j=2,nyfin-1
+              DO i=2,nxfin-1
+                 IF (G1%bathy_meter(i,j).GT.0.) THEN
+                    IF (MAX(bathy_test(i,j+1),bathy_test(i,j-1),bathy_test(i-1,j),bathy_test(i+1,j)).EQ.1) THEN
+                       IF (bathy_test(i,j).NE.1.) nbadd = nbadd + 1
+                       bathy_test(i,j)=1.
                     ENDIF
-                 ENDDO
+
+                 ENDIF
               ENDDO
            ENDDO
-           WHERE (bathy_test.EQ.0.)
-              G1%bathy_meter = 0.
-           END WHERE
-           DEALLOCATE(bathy_test)
-        ENDIF
-        !
-        ! Chanut: Compute partial step bathymetry once more
-        CALL get_partial_steps(G1)                 ! compute gdepw_ps for G1
-
-        IF(bathy_update) CALL Update_Parent_Bathy( G0,G1 ) 
-        !
+        ENDDO
+        WHERE (bathy_test.EQ.0.)   G1%bathy_meter = 0.
+        DEALLOCATE(bathy_test)
+     ENDIF
+     !
+     !                              ! ----------------
+     IF( partial_steps ) THEN       ! If partial steps
+        !                           ! ----------------
+        !
+        CALL get_partial_steps(G1)  ! recompute gdepw_ps for G1
+
+        IF(bathy_update)   CALL Update_Parent_Bathy( G0,G1 ) 
         CALL set_child_name(parent_bathy_meter,child_ps)
-        status = Write_Bathy_meter(TRIM(child_ps),G1)        
-
-        IF(bathy_update) status = Write_Bathy_meter(TRIM(updated_parent_file),G0)
+        
+                           status = Write_Bathy_meter(TRIM(child_ps),G1)
+        IF(bathy_update)   status = Write_Bathy_meter(TRIM(updated_parent_file),G0)
+        IF(bathy_update)   status = write_domcfg(TRIM(updated_parent_domcfg),G0)
 
         CALL get_partial_steps(G1)
@@ -640,119 +605,37 @@
         G1%bathy_level=NINT(G1%bathy_level)
         !
-        IF( TRIM(parent_meshmask_file) .NE. '/NULL' ) &
-             status = Write_Bathy_level(TRIM(Childlevel_file),G1)
+        ! store interpolation result in output file
+        CALL levels_to_meter(G1)   ! From levels to meters
+        IF( TRIM(parent_meshmask_file) .NE. '/NULL' )    status = Write_Bathy_level(TRIM(Childlevel_file),G1)
+        IF( TRIM(parent_meshmask_file) .NE. '/NULL' )    status = write_domcfg(TRIM(Child_domcfg),G1)
         !
         WRITE(*,*) '****** Bathymetry successfully created for partial cells ******'
-        WRITE(*,*) ' '
-        !
-        STOP         
-     ENDIF
-     !            
-     !--------------------------------end if partial steps------------------------------------
-     !
-     !
-     status = Read_bathy_level(TRIM(parent_meshmask_file),G0)
-     !            
-     CALL levels_to_meter(G0)
-     ! 
-     ! compute constant bathymetry for Parent-Child bathymetry connection
-     !              
-     WHERE( G0%bathy_meter < 0.000001 ) G0%bathy_meter = 0.
-
-     CALL init_constant_bathy(G0%bathy_meter,bathy_fin_constant)
-     !
-     boundary = npt_copy*irafx + nbghostcellsfine + 1   
-     !             
-     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)
-     !
-     ! bathymetry smoothing
-     !                  
-     CALL smooth_topo(G1%bathy_meter(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter)
-     !
-     ! convert bathymetry from meters to levels
-     !
-     CALL meter_to_levels(G1) 
-     !           
-     DEALLOCATE(G1%bathy_meter)           
-     !	     	     
+        STOP
+        !
+        !                           ! ----------------
+     ELSE                           ! No partial steps
+        !                           ! ----------------
+        !
+        ! convert bathymetry from meters to levels
+        CALL meter_to_levels(G1) 
+
+        IF(bathy_update)   CALL Update_Parent_Bathy( G0,G1 ) 
+        CALL set_child_name(parent_bathy_meter,child_ps)
+        
+                           status = Write_Bathy_meter(TRIM(child_ps),G1)        
+        IF(bathy_update)   status = Write_Bathy_meter(TRIM(updated_parent_file),G0)
+        IF(bathy_update)   status = write_domcfg(TRIM(updated_parent_domcfg),G0)
+        !
+        ! store interpolation result in output file
+        CALL levels_to_meter(G1)   ! From levels to meters
+        status = Write_Bathy_level(TRIM(Childlevel_file),G1)
+        status = write_domcfg(TRIM(Child_domcfg),G1)
+        !
+        WRITE(*,*) '****** Bathymetry successfully created for full cells ******'
+        STOP
+        ! 
+     ENDIF
   ENDIF
   !
-  !
-  ! make connection thanks to constant and interpolated bathymetry
-  !
-  !      
-  G1%bathy_level=NINT(G1%bathy_level)
-  !	    
-  DO j=1,nyfin
-     DO i=1,nxfin
-        IF (g1%bathy_level(i,j).LT.0.) THEN
-           PRINT *,'error in ',i,j,g1%bathy_level(i,j)
-        ENDIF
-     ENDDO
-  ENDDO
-  !	    
-  WHERE ((G1%bathy_level.LT.3.).AND.(G1%bathy_level.GT.0.))
-     G1%bathy_level=3
-  END WHERE
-  !
-  ! possibility to 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
-     END WHERE
-
-     WHERE (G1%bathy_level(nxfin,:).GT.0.)
-        bathy_test(nxfin,:)=1
-     END WHERE
-
-     WHERE (G1%bathy_level(:,1).GT.0.)
-        bathy_test(:,1)=1
-     END WHERE
-
-     WHERE (G1%bathy_level(:,nyfin).GT.0.)
-        bathy_test(:,nyfin)=1
-     END WHERE
-
-     nbadd = 1
-
-     DO WHILE (nbadd.NE.0)
-        nbadd = 0
-        DO j=2,nyfin-1
-           DO i=2,nxfin-1
-              IF (G1%bathy_level(i,j).GT.0.) THEN
-                 IF (MAX(bathy_test(i,j+1),bathy_test(i,j-1),bathy_test(i-1,j),bathy_test(i+1,j)).EQ.1) THEN
-                    IF (bathy_test(i,j).NE.1.) nbadd = nbadd + 1
-                    bathy_test(i,j)=1.
-                 ENDIF
-
-              ENDIF
-           ENDDO
-        ENDDO
-
-     ENDDO
-
-     WHERE (bathy_test.EQ.0.)
-        G1%bathy_level = 0.
-     END WHERE
-     DEALLOCATE(bathy_test)           
-  ENDIF
-
-
-  !
-  ! store interpolation result in output file
-  !
-  status = Write_Bathy_level(TRIM(Childlevel_file),G1)
-
-  WRITE(*,*) '****** Bathymetry successfully created for full cells ******'
-  WRITE(*,*) ' '
-
-  STOP
 END PROGRAM create_bathy
 
Index: /utils/tools/NESTING/src/agrif_create_coordinates.f90
===================================================================
--- /utils/tools/NESTING/src/agrif_create_coordinates.f90	(revision 10024)
+++ /utils/tools/NESTING/src/agrif_create_coordinates.f90	(revision 10025)
@@ -39,15 +39,10 @@
   !
   ! read input file (namelist.input)
-  !
   CALL read_namelist(namelistname)
-  !
-
   !      
   ! read parent coodinates file
-  !
   status = Read_Coordinates(TRIM(parent_coordinate_file),G0)
   ! 
   ! define name of child coordinate file
-  !     
   CALL set_child_name(parent_coordinate_file,Child_filename) 
   !
@@ -64,20 +59,11 @@
   !
   ! allocation of child grid elements
-  !
   CALL agrif_grid_allocate(G1,nxfin,nyfin)
   !
-  !
   !  check potential longitude problems 
-  !      
   IF( G0%glamt(imin,jmin) > G0%glamt(imax,jmax) ) THEN           
-     WRITE(*,*) '    '
-     WHERE ( G0%glamt < 0 )
-        G0%glamt = G0%glamt + 360.
-     END WHERE
-     WHERE ( G0%glamf < 0 )
-        G0%glamf = G0%glamf + 360.
-     END WHERE
-  ENDIF
-
+     WHERE ( G0%glamt < 0 )   G0%glamt = G0%glamt + 360.
+     WHERE ( G0%glamf < 0 )   G0%glamf = G0%glamf + 360.
+  ENDIF
   !      
   ! interpolation from parent grid to child grid for
@@ -89,11 +75,7 @@
   !    gphi = latitude
   !
-
   !      
-  !> M. Dunphy ticket 2082:
   CALL agrif_interp(G0%glamt,G1%glamt,'T')
   CALL agrif_interp(G0%glamf,G1%glamf,'F')
-!  G1%glamu = G1%glamf
-!  G1%glamv = G1%glamt
   CALL agrif_interp(G0%glamu,G1%glamu,'U')
   CALL agrif_interp(G0%glamv,G1%glamv,'V') 
@@ -101,9 +83,6 @@
   CALL agrif_interp(G0%gphit,G1%gphit,'T')
   CALL agrif_interp(G0%gphif,G1%gphif,'F')
-!  G1%gphiu = G1%gphit
-!  G1%gphiv = G1%gphif
   CALL agrif_interp(G0%gphiu,G1%gphiu,'U')
   CALL agrif_interp(G0%gphiv,G1%gphiv,'V')
-  !<  M. Dunphy ticket 2082
   !
   !
@@ -113,103 +92,88 @@
   !
   ! Compute scale factors e1 e2
-  !
-  DO j=1,nyfin
-     DO i=2,nxfin
-        G1%e1t(i,j) = ra * rad * SQRT( (COS(rad*G1%gphit(i,j))*(G1%glamu(i,j)-G1%glamu(i-1,j)))**2 &
-             + (G1%gphiu(i,j)-G1%gphiu(i-1,j))**2)
-        G1%e1v(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiv(i,j))*(G1%glamf(i,j)-G1%glamf(i-1,j)))**2 &
-             + (G1%gphif(i,j)-G1%gphif(i-1,j))**2)                                                             
-     END DO
-  END DO
-  !      
-  G1%e1t(1,:)=G1%e1t(2,:)
-  G1%e1v(1,:)=G1%e1v(2,:)
-  !
-  DO j=1,nyfin
-     DO i=1,nxfin-1
-        G1%e1u(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiu(i,j))*(G1%glamt(i+1,j)-G1%glamt(i,j)))**2 &
-             + (G1%gphit(i+1,j)-G1%gphit(i,j))**2)
-        G1%e1f(i,j) = ra * rad * SQRT( (COS(rad*G1%gphif(i,j))*(G1%glamv(i+1,j)-G1%glamv(i,j)))**2 &
-             + (G1%gphiv(i+1,j)-G1%gphiv(i,j))**2)                                                            
-     END DO
-  END DO
-  !
-  G1%e1u(nxfin,:)=G1%e1u(nxfin-1,:)
-  G1%e1f(nxfin,:)=G1%e1f(nxfin-1,:)
-  !
-  DO j=2,nyfin
-     DO i=1,nxfin                                     
-        G1%e2t(i,j) = ra * rad * SQRT( (COS(rad*G1%gphit(i,j))*(G1%glamv(i,j)-G1%glamv(i,j-1)))**2 &
-             + (G1%gphiv(i,j)-G1%gphiv(i,j-1))**2)
-        G1%e2u(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiu(i,j))*(G1%glamf(i,j)-G1%glamf(i,j-1)))**2 &
-             + (G1%gphif(i,j)-G1%gphif(i,j-1))**2)                                                               
-     END DO
-  END DO
-  !
-  G1%e2t(:,1)=G1%e2t(:,2)
-  G1%e2u(:,1)=G1%e2u(:,2)
-  !
-  DO j=1,nyfin-1
-     DO i=1,nxfin
-        G1%e2v(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiv(i,j))*(G1%glamt(i,j+1)-G1%glamt(i,j)))**2 &
-             + (G1%gphit(i,j+1)-G1%gphit(i,j))**2)
-        G1%e2f(i,j) = ra * rad * SQRT( (COS(rad*G1%gphif(i,j))*(G1%glamu(i,j+1)-G1%glamu(i,j)))**2 &
-             + (G1%gphiu(i,j+1)-G1%gphiu(i,j))**2)
-     END DO
-  END DO
-  !
-  G1%e2v(:,nyfin)=G1%e2v(:,nyfin-1)
-  G1%e2f(:,nyfin)=G1%e2f(:,nyfin-1)
-
-
-  CALL agrif_interp(G0%e1t,G1%e1t,'T')
-  G1%e1t = G1%e1t / REAL(irafx)
-  CALL agrif_interp(G0%e2t,G1%e2t,'T')
-  G1%e2t = G1%e2t / REAL(irafy)
-
-  CALL agrif_interp(G0%e1u,G1%e1u,'U')
-  G1%e1u = G1%e1u / REAL(irafx)
-  CALL agrif_interp(G0%e2u,G1%e2u,'U')
-  G1%e2u = G1%e2u / REAL(irafy)  
-
-  CALL agrif_interp(G0%e1v,G1%e1v,'V')
-  G1%e1v = G1%e1v / REAL(irafx)
-  CALL agrif_interp(G0%e2v,G1%e2v,'V')
-  G1%e2v = G1%e2v / REAL(irafy) 
-
-  CALL agrif_interp(G0%e1f,G1%e1f,'F')
-  G1%e1f = G1%e1f / REAL(irafx)
-  CALL agrif_interp(G0%e2f,G1%e2f,'F')
-  G1%e2f = G1%e2f / REAL(irafy)               
-
-  !
-  WHERE ( G1%glamt > 180 )
-     G1%glamt = G1%glamt - 360.
-  END WHERE
-  WHERE ( G1%glamf > 180 )
-     G1%glamf = G1%glamf - 360.
-  END WHERE
-  WHERE ( G1%glamu > 180 )
-     G1%glamu = G1%glamu - 360.
-  END WHERE
-  WHERE ( G1%glamv > 180 )
-     G1%glamv = G1%glamv - 360.
-  END WHERE
-  !
+!  DO j=1,nyfin
+!     DO i=2,nxfin
+!        G1%e1t(i,j) = ra * rad * SQRT( (COS(rad*G1%gphit(i,j))*(G1%glamu(i,j)-G1%glamu(i-1,j)))**2 &
+!             + (G1%gphiu(i,j)-G1%gphiu(i-1,j))**2)
+!        G1%e1v(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiv(i,j))*(G1%glamf(i,j)-G1%glamf(i-1,j)))**2 &
+!             + (G1%gphif(i,j)-G1%gphif(i-1,j))**2)                                                             
+!     END DO
+!  END DO
+!  !      
+!  G1%e1t(1,:)=G1%e1t(2,:)
+!  G1%e1v(1,:)=G1%e1v(2,:)
+!  !
+!  DO j=1,nyfin
+!     DO i=1,nxfin-1
+!        G1%e1u(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiu(i,j))*(G1%glamt(i+1,j)-G1%glamt(i,j)))**2 &
+!             + (G1%gphit(i+1,j)-G1%gphit(i,j))**2)
+!        G1%e1f(i,j) = ra * rad * SQRT( (COS(rad*G1%gphif(i,j))*(G1%glamv(i+1,j)-G1%glamv(i,j)))**2 &
+!             + (G1%gphiv(i+1,j)-G1%gphiv(i,j))**2)                                                            
+!     END DO
+!  END DO
+!  !
+!  G1%e1u(nxfin,:)=G1%e1u(nxfin-1,:)
+!  G1%e1f(nxfin,:)=G1%e1f(nxfin-1,:)
+!  !
+!  DO j=2,nyfin
+!     DO i=1,nxfin                                     
+!        G1%e2t(i,j) = ra * rad * SQRT( (COS(rad*G1%gphit(i,j))*(G1%glamv(i,j)-G1%glamv(i,j-1)))**2 &
+!             + (G1%gphiv(i,j)-G1%gphiv(i,j-1))**2)
+!        G1%e2u(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiu(i,j))*(G1%glamf(i,j)-G1%glamf(i,j-1)))**2 &
+!             + (G1%gphif(i,j)-G1%gphif(i,j-1))**2)                                                               
+!     END DO
+!  END DO
+!  !
+!  G1%e2t(:,1)=G1%e2t(:,2)
+!  G1%e2u(:,1)=G1%e2u(:,2)
+!  !
+!  DO j=1,nyfin-1
+!     DO i=1,nxfin
+!        G1%e2v(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiv(i,j))*(G1%glamt(i,j+1)-G1%glamt(i,j)))**2 &
+!             + (G1%gphit(i,j+1)-G1%gphit(i,j))**2)
+!        G1%e2f(i,j) = ra * rad * SQRT( (COS(rad*G1%gphif(i,j))*(G1%glamu(i,j+1)-G1%glamu(i,j)))**2 &
+!             + (G1%gphiu(i,j+1)-G1%gphiu(i,j))**2)
+!     END DO
+!  END DO
+!  !
+!  G1%e2v(:,nyfin)=G1%e2v(:,nyfin-1)
+!  G1%e2f(:,nyfin)=G1%e2f(:,nyfin-1)
+
+
+  CALL agrif_interp(G0%e1t,G1%e1t,'T')   ;   G1%e1t = G1%e1t / REAL(irafx)
+  CALL agrif_interp(G0%e2t,G1%e2t,'T')   ;   G1%e2t = G1%e2t / REAL(irafy)
+
+  CALL agrif_interp(G0%e1u,G1%e1u,'U')   ;   G1%e1u = G1%e1u / REAL(irafx)
+  CALL agrif_interp(G0%e2u,G1%e2u,'U')   ;   G1%e2u = G1%e2u / REAL(irafy)  
+
+  CALL agrif_interp(G0%e1v,G1%e1v,'V')   ;   G1%e1v = G1%e1v / REAL(irafx)
+  CALL agrif_interp(G0%e2v,G1%e2v,'V')   ;   G1%e2v = G1%e2v / REAL(irafy) 
+
+  CALL agrif_interp(G0%e1f,G1%e1f,'F')   ;   G1%e1f = G1%e1f / REAL(irafx)
+  CALL agrif_interp(G0%e2f,G1%e2f,'F')   ;   G1%e2f = G1%e2f / REAL(irafy)
+  !
+  WHERE ( G1%glamt > 180 )   G1%glamt = G1%glamt - 360.
+  WHERE ( G1%glamf > 180 )   G1%glamf = G1%glamf - 360.
+  WHERE ( G1%glamu > 180 )   G1%glamu = G1%glamu - 360.
+  WHERE ( G1%glamv > 180 )   G1%glamv = G1%glamv - 360.
   !
   G1%nav_lon=G1%glamt
   G1%nav_lat=G1%gphit 
   !              
-  !
   ! Write interpolation result in child coodinates file
-  !      
   status = Write_Coordinates(TRIM(Child_filename),G1)
-
   !
   WRITE(*,*) 'Child domain position : '
-  WRITE(*,*) 'latmin =',G1%gphit(3,3)
-  WRITE(*,*) 'latmax =',G1%gphit(nxfin-2,nyfin-2)
-  WRITE(*,*) 'lonmin =',G1%glamt(3,3)
-  WRITE(*,*) 'lonmax =',G1%glamt(nxfin-2,nyfin-2)
+  IF( ln_agrif_domain ) THEN
+     WRITE(*,*) 'latmin =',G1%gphit(3,3)
+     WRITE(*,*) 'latmax =',G1%gphit(nxfin-2,nyfin-2)
+     WRITE(*,*) 'lonmin =',G1%glamt(3,3)
+     WRITE(*,*) 'lonmax =',G1%glamt(nxfin-2,nyfin-2)
+  ELSE
+     WRITE(*,*) 'latmin =',G1%gphit(1,1)
+     WRITE(*,*) 'latmax =',G1%gphit(nxfin,nyfin)
+     WRITE(*,*) 'lonmin =',G1%glamt(1,1)
+     WRITE(*,*) 'lonmax =',G1%glamt(nxfin,nyfin)
+  ENDIF
   STOP
 END PROGRAM create_coordinate
Index: /utils/tools/NESTING/src/agrif_create_data.f90
===================================================================
--- /utils/tools/NESTING/src/agrif_create_data.f90	(revision 10024)
+++ /utils/tools/NESTING/src/agrif_create_data.f90	(revision 10025)
@@ -10,16 +10,16 @@
   !************************************************************************
   ! 									*
-  ! PROGRAM  CREATE_DATA							*
+  ! PROGRAM  CREATE_DATA					        *
   !									*
   ! program to implement data interpolation to generate	 		*
   ! child grid forcing files						*
-  !									*									*
+  !									*
   !Interpolation is carried out using bilinear interpolation		*
   !routine from SCRIP package						*		
   !									*
-  !http://climate.lanl.gov/Software/SCRIP/				*									*
+  !http://climate.lanl.gov/Software/SCRIP/				*
   !************************************************************************
   !
-  INTEGER :: narg,iargc,i
+  INTEGER           :: narg, iargc, ji
   CHARACTER(len=80) :: namelistname
 
@@ -33,40 +33,34 @@
 
   ! read input file (namelist.input)
-  !
   CALL read_namelist(namelistname)
-
-  i = 1
   !
   ! Interpolate U grid  data
-  !
-  DO WHILE( TRIM(U_Files(i)) .NE. '/NULL' )
-     PRINT *,'Grid U forcing files = ',u_files(i)
+  ji = 1
+  DO WHILE( TRIM(U_Files(ji)) .NE. '/NULL' )
+     PRINT *,'Grid U forcing files = ',u_files(ji)
      !        
-     CALL Interp_Extrap_var(U_FILES(i), 'U') 
-     i = i+1                
+     CALL Interp_Extrap_var(U_FILES(ji), 'U') 
+     ji = ji+1                
      !             
   END DO
 
-  i = 1
   !
   ! Interpolate V grid  data
-  !
-  DO WHILE( TRIM(V_Files(i)) .NE. '/NULL' )
-     PRINT *,'Grid V forcing files = ',v_files(i)
+  ji = 1
+  DO WHILE( TRIM(V_Files(ji)) .NE. '/NULL' )
+     PRINT *,'Grid V forcing files = ',v_files(ji)
      !        
-     CALL Interp_Extrap_var(V_FILES(i), 'V') 
-     i = i+1                
+     CALL Interp_Extrap_var(V_FILES(ji), 'V') 
+     ji = ji+1                
      !             
   END DO
-
-  i = 1
   !
   ! Interpolate flux data
-  !
-  DO WHILE( TRIM(Flx_Files(i)) .NE. '/NULL' )
-     PRINT *,'flxfiles = ',flx_files(i)
+  ji = 1
+  DO WHILE( TRIM(Flx_Files(ji)) .NE. '/NULL' )
+     PRINT *,'flxfiles = ',flx_files(ji)
      !        
-     CALL Interp_Extrap_var(FLX_FILES(i), 'T') 
-     i = i+1                
+     CALL Interp_Extrap_var(FLX_FILES(ji), 'T') 
+     ji = ji+1                
      !             
   END DO
Index: /utils/tools/NESTING/src/agrif_create_restart.f90
===================================================================
--- /utils/tools/NESTING/src/agrif_create_restart.f90	(revision 10024)
+++ /utils/tools/NESTING/src/agrif_create_restart.f90	(revision 10025)
@@ -389,5 +389,5 @@
            WRITE(*,*) 'copy time_steps'
            CALL Read_Ncdf_var('time_steps',TRIM(restart_file),tabtemp1DInt) 
-           CALL Write_Ncdf_var('time_steps',TRIM(timedimname),Child_file,tabtemp1DInt)
+           CALL Write_Ncdf_var('time_steps',TRIM(timedimname),Child_file,tabtemp1DInt,'integer')
            CALL Copy_Ncdf_att('time_steps',TRIM(restart_file),Child_file) 
            DEALLOCATE(tabtemp1DInt)
Index: /utils/tools/NESTING/src/agrif_interpolation.f90
===================================================================
--- /utils/tools/NESTING/src/agrif_interpolation.f90	(revision 10024)
+++ /utils/tools/NESTING/src/agrif_interpolation.f90	(revision 10025)
@@ -54,12 +54,13 @@
     ! 
     !
-    CALL agrif_base_interp2(tabin,tabout,imin,jmin,typevar)
-
+    IF( ln_agrif_domain ) THEN
+       CALL agrif_base_interp2(tabin,tabout,imin,jmin,typevar)
+    ELSE
+       CALL agrif_base_interp3(tabin,tabout,imin,jmin,typevar)
+    ENDIF
+    
     !
   END SUBROUTINE agrif_interp
   !
-  !********************************************************
-  !   subroutine agrif_base_interp2			*
-  !********************************************************
   !      
   SUBROUTINE agrif_base_interp2(tabin,tabout,i_min,j_min,typevar)
@@ -189,7 +190,70 @@
   END SUBROUTINE agrif_base_interp2
   !
-  !********************************************************
-  !   subroutine polint					*
-  !********************************************************
+  !      
+  SUBROUTINE agrif_base_interp3(tabin,tabout,i_min,j_min,typevar)
+    !
+    IMPLICIT NONE
+    REAL*8,DIMENSION(:,:) :: tabin,tabout
+    INTEGER :: i_min,j_min
+    CHARACTER(*) :: typevar   
+
+    INTEGER :: nxf,nyf,zx,zy
+    INTEGER :: ji,jj,jif,jjf,jic,jjc,jdecx,jdecy
+    REAL*8 :: Ax, Bx, Ay, By
+
+    nxf = SIZE(tabout,1)
+    nyf = SIZE(tabout,2)
+
+    SELECT CASE(typevar)
+    CASE('T')
+       IF(MOD(irafx,2)==1) THEN ! odd
+          zx = 1 ; zy = 1 ; jdecx = FLOOR(irafx/2.) ; jdecy = FLOOR(irafy/2.)
+       ELSE                     ! even
+          zx = 2 ; zy = 2 ; jdecx = FLOOR(irafx/2.) ; jdecy = FLOOR(irafy/2.)
+       ENDIF
+    CASE('U')
+       IF(MOD(irafx,2)==1) THEN ! odd
+          zx = 1 ; zy = 1 ; jdecx = irafx - 1 ; jdecy = FLOOR(irafy/2.)
+       ELSE                     ! even
+          zx = 1 ; zy = 2 ; jdecx = irafx - 1 ; jdecy = FLOOR(irafy/2.)
+       ENDIF
+    CASE('V')
+       IF(MOD(irafx,2)==1) THEN ! odd
+          zx = 1 ; zy = 1 ; jdecx = FLOOR(irafx/2.) ; jdecy = irafy - 1
+       ELSE                     ! even
+          zx = 2 ; zy = 1 ; jdecx = FLOOR(irafx/2.) ; jdecy = irafy - 1
+       ENDIF
+    CASE('F')
+       IF(MOD(irafx,2)==1) THEN ! odd
+          zx = 1 ; zy = 1 ; jdecx = irafx - 1 ; jdecy = irafy - 1
+       ELSE                     ! even
+          zx = 1 ; zy = 1 ; jdecx = irafx - 1 ; jdecy = irafy - 1
+       ENDIF
+    END SELECT
+
+    
+    DO jj = 1, nyf
+
+       jjf = jj - jdecy
+       jjc = j_min + FLOOR((jjf-1.) / 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)
+          Ax = 1. - Bx
+          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
+       END DO
+    END DO
+    
+    !
+  END SUBROUTINE agrif_base_interp3
+
   !       
   SUBROUTINE polint(xin,valin,n,x,val)
@@ -238,7 +302,4 @@
   END SUBROUTINE polint
   !
-  !********************************************************
-  !   subroutine polcoe					*
-  !********************************************************
   !      
   SUBROUTINE polcoe(xin,valin,n,cof)
@@ -709,5 +770,5 @@
           CALL Read_Ncdf_var('time_steps',filename,tabtemp1DInt) 
           !		   print *,'timedeph = ',tabtemp1DInt
-          CALL Write_Ncdf_var('time_steps','time_counter',Child_file,tabtemp1DInt)
+          CALL Write_Ncdf_var('time_steps','time_counter',Child_file,tabtemp1DInt,'integer')
           CALL Copy_Ncdf_att('time_steps',filename,Child_file) 
           DEALLOCATE(tabtemp1DInt)
@@ -984,7 +1045,4 @@
   END SUBROUTINE Interp_Extrap_var
   !
-  !**************************************************************
-  !   end subroutine Interp_var
-  !**************************************************************          
   ! 
 END MODULE agrif_interpolation
Index: /utils/tools/NESTING/src/agrif_partial_steps.f90
===================================================================
--- /utils/tools/NESTING/src/agrif_partial_steps.f90	(revision 10024)
+++ /utils/tools/NESTING/src/agrif_partial_steps.f90	(revision 10025)
@@ -276,6 +276,5 @@
     REAL,DIMENSION(N) :: gdepw,e3t
     REAL :: za0,za1,za2,zsur,zacr,zacr2,zkth,zkth2,zmin,zmax,zdepth
-    INTEGER :: kbathy,jk,diff
-    INTEGER :: bornex,borney,bornex2,borney2
+    INTEGER :: kbathy,jk
     !
     IF ( ( pa0 == 0 .OR. pa1 == 0 .OR. psur == 0 ) &
@@ -348,40 +347,33 @@
        e3t(N) = e3t(N-1)
     END IF
-    !
-    diff = 0      
-    IF ( MOD(irafx,2) .EQ. 0 ) diff = 1
-    !       
-    bornex = nbghostcellsfine + 1 + CEILING(irafx/2.0) + diff - irafx
-    borney = nbghostcellsfine + 1 + CEILING(irafy/2.0) + diff - irafy
-    bornex2 = nxfin - nbghostcellsfine - irafx - CEILING(irafx/2.0) 
-    borney2 = nyfin - nbghostcellsfine - irafy - CEILING(irafy/2.0)                      
-    !
+    !              
     !
     ! west boundary
-    !
-
-    CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,2+nbghostcellsfine+(npt_copy+npt_connect)*irafx-1, &
-         1,nyfin)
-
+    IF( ln_agrif_domain ) THEN
+       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,2+nbghostcellsfine+(npt_copy+npt_connect)*irafx-1,1,nyfin)
+    ELSE
+       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,(npt_copy+npt_connect)*irafx,1,nyfin)
+    ENDIF
     !
     ! east boundary
-    !
-
-    CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,nxfin-1-nbghostcellsfine-((npt_copy+npt_connect)*irafx-1),nxfin, &
-         1,nyfin)
-
+    IF( ln_agrif_domain ) THEN
+       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,nxfin-1-nbghostcellsfine-((npt_copy+npt_connect)*irafx-1),nxfin,1,nyfin)
+    ELSE
+       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,nxfin-((npt_copy+npt_connect)*irafx+1),nxfin,1,nyfin)
+    ENDIF
     !
     ! north boundary
-    !
-
-    CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,nxfin, &
-         nyfin-1 - nbghostcellsfine -((npt_copy+npt_connect)*irafy-1),nyfin )
-
+    IF( ln_agrif_domain ) THEN
+       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,nxfin,nyfin-1-nbghostcellsfine-((npt_copy+npt_connect)*irafy-1),nyfin)
+    ELSE
+       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,nxfin,nyfin-((npt_copy+npt_connect)*irafy+1),nyfin)
+    ENDIF
     !
     ! south boundary
-    !
-    CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,nxfin, &
-         1,2+nbghostcellsfine+(npt_copy+npt_connect)*irafy-1 )
-
+    IF( ln_agrif_domain ) THEN
+       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,nxfin,1,2+nbghostcellsfine+(npt_copy+npt_connect)*irafy-1)
+    ELSE
+       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,nxfin,1,(npt_copy+npt_connect)*irafy)
+    ENDIF
     !       
     !
@@ -398,5 +390,5 @@
     INTEGER :: kbathy,jk,indx,indy,diff
     REAL :: xdiff
-    INTEGER :: i,j,ji,ij,ii,jj,jpt,ipt
+    INTEGER :: i,j,ji,ij,ii,jj,jpt,ipt,i1,i2,j1,j2,ii1,ii2,jj1,jj2
     REAL*8 :: slopex, slopey,val,tmp1,tmp2,tmp3,tmp4
     INTEGER :: parentbathy
@@ -471,8 +463,43 @@
     jend = pty + agrif_int((y-ymin-dyfin/2.)/dyfin)               
 
-
-    ALLOCATE(gdepwtemp(ibegin-irafx:iend+irafx,jbegin-irafy:jend+irafy))
-    ALLOCATE(parentbathytab(ibegin-irafx:iend+irafx,jbegin-irafy:jend+irafy))
-
+    IF( ln_agrif_domain ) THEN
+       ALLOCATE(gdepwtemp(ibegin-irafx:iend+irafx,jbegin-irafy:jend+irafy))
+       ALLOCATE(parentbathytab(ibegin-irafx:iend+irafx,jbegin-irafy:jend+irafy))
+
+       i1 = ibegin
+       i2 = iend
+       j1 = jbegin
+       j2 = jend
+       
+       ii1 = -FLOOR(irafx/2.0)+diff
+       ii2 =  FLOOR(irafx/2.0)
+       jj1 = -FLOOR(irafy/2.0)+diff
+       jj2 =  FLOOR(irafy/2.0)
+    ELSE
+       ibegin = minboundx
+       jbegin = minboundy
+       iend   = maxboundx ! (npt_copy+npt_connect)*irafx
+       jend   = maxboundy ! (npt_copy+npt_connect)*irafy
+       !
+       ipbegin = imin + (ibegin-1)/irafx
+       jpbegin = jmin + (jbegin-1)/irafy
+       ipend   = ipbegin + (npt_copy+npt_connect) - 1
+       jpend   = jpbegin + (npt_copy+npt_connect) - 1
+       !
+       i1 = ibegin
+       i2 = iend
+       j1 = jbegin
+       j2 = jend
+       
+       ii1 = 0
+       ii2 = irafx - 1
+       jj1 = 0
+       jj2 = irafy - 1
+       !
+       ALLOCATE(gdepwtemp(ibegin:iend,jbegin:jend))
+       ALLOCATE(parentbathytab(ibegin:iend,jbegin:jend))
+
+    ENDIF
+    
 
     jpt=jpbegin
@@ -482,5 +509,5 @@
 
 
-       DO i=ibegin,iend,irafx
+       DO i=i1,i2,irafx
 
 
@@ -524,6 +551,6 @@
           ! interpolation on fine grid points (connection zone)
           !
-          DO ii = i-FLOOR(irafx/2.0)+diff,i+FLOOR(irafx/2.0)
-             x = ii-i - xdiff/2.
+          DO ii = i+ii1,i+ii2
+!!             x = ii-i - xdiff/2.
 !!             val = parentgrid%gdepw_ps(ipt,jpt)+slopex * x
 !! chanut: uncomment this to get nearest neighbor interpolation
@@ -563,5 +590,5 @@
 
           slopey = vanleer(gdepwtemp(i,j-irafy:j+irafy:irafy))/REAL(irafy)
-
+          
           tmp1 = (maxdepth - gdepwtemp(i,j)) / REAL(irafy)
           tmp2 = (gdepwtemp(i,j) - mindepth) / REAL(irafy)
@@ -583,6 +610,6 @@
 
 
-          DO jj = j-FLOOR(irafy/2.0)+diff,j+FLOOR(irafy/2.0)
-             y = jj-j - xdiff/2.
+          DO jj = j+jj1,j+jj2
+!!             y = jj-j - xdiff/2.
 !!             val = gdepwtemp(i,j) + slopey*y
 !! chanut: uncomment this to get nearest neighbor interpolation
Index: /utils/tools/NESTING/src/agrif_readwrite.f90
===================================================================
--- /utils/tools/NESTING/src/agrif_readwrite.f90	(revision 10024)
+++ /utils/tools/NESTING/src/agrif_readwrite.f90	(revision 10025)
@@ -40,22 +40,22 @@
     TYPE(Coordinates) :: Grid
     !      
-    CALL Read_Ncdf_var('glamt',name,Grid%glamt)
-    CALL Read_Ncdf_var('glamu',name,Grid%glamu)
-    CALL Read_Ncdf_var('glamv',name,Grid%glamv)
-    CALL Read_Ncdf_var('glamf',name,Grid%glamf)
-    CALL Read_Ncdf_var('gphit',name,Grid%gphit)
-    CALL Read_Ncdf_var('gphiu',name,Grid%gphiu)
-    CALL Read_Ncdf_var('gphiv',name,Grid%gphiv)
-    CALL Read_Ncdf_var('gphif',name,Grid%gphif)
-    CALL Read_Ncdf_var('e1t',name,Grid%e1t)
-    CALL Read_Ncdf_var('e1u',name,Grid%e1u)
-    CALL Read_Ncdf_var('e1v',name,Grid%e1v)
-    CALL Read_Ncdf_var('e1f',name,Grid%e1f)
-    CALL Read_Ncdf_var('e2t',name,Grid%e2t)
-    CALL Read_Ncdf_var('e2u',name,Grid%e2u)
-    CALL Read_Ncdf_var('e2v',name,Grid%e2v)
-    CALL Read_Ncdf_var('e2f',name,Grid%e2f)
-    CALL Read_Ncdf_var('nav_lon',name,Grid%nav_lon)
-    CALL Read_Ncdf_var('nav_lat',name,Grid%nav_lat)       
+    CALL read_ncdf_var('glamt',name,Grid%glamt)
+    CALL read_ncdf_var('glamu',name,Grid%glamu)
+    CALL read_ncdf_var('glamv',name,Grid%glamv)
+    CALL read_ncdf_var('glamf',name,Grid%glamf)
+    CALL read_ncdf_var('gphit',name,Grid%gphit)
+    CALL read_ncdf_var('gphiu',name,Grid%gphiu)
+    CALL read_ncdf_var('gphiv',name,Grid%gphiv)
+    CALL read_ncdf_var('gphif',name,Grid%gphif)
+    CALL read_ncdf_var('e1t',name,Grid%e1t)
+    CALL read_ncdf_var('e1u',name,Grid%e1u)
+    CALL read_ncdf_var('e1v',name,Grid%e1v)
+    CALL read_ncdf_var('e1f',name,Grid%e1f)
+    CALL read_ncdf_var('e2t',name,Grid%e2t)
+    CALL read_ncdf_var('e2u',name,Grid%e2u)
+    CALL read_ncdf_var('e2v',name,Grid%e2v)
+    CALL read_ncdf_var('e2f',name,Grid%e2f)
+    CALL read_ncdf_var('nav_lon',name,Grid%nav_lon)
+    CALL read_ncdf_var('nav_lat',name,Grid%nav_lat)       
     ! 
     IF( PRESENT(Pacifique) )THEN
@@ -103,22 +103,22 @@
     TYPE(Coordinates) :: Grid
     !      
-    CALL Read_Ncdf_var('glamt',name,Grid%glamt,strt,cnt)
-    CALL Read_Ncdf_var('glamu',name,Grid%glamu,strt,cnt)
-    CALL Read_Ncdf_var('glamv',name,Grid%glamv,strt,cnt)
-    CALL Read_Ncdf_var('glamf',name,Grid%glamf,strt,cnt)
-    CALL Read_Ncdf_var('gphit',name,Grid%gphit,strt,cnt)
-    CALL Read_Ncdf_var('gphiu',name,Grid%gphiu,strt,cnt)
-    CALL Read_Ncdf_var('gphiv',name,Grid%gphiv,strt,cnt)
-    CALL Read_Ncdf_var('gphif',name,Grid%gphif,strt,cnt)
-    CALL Read_Ncdf_var('e1t',name,Grid%e1t,strt,cnt)
-    CALL Read_Ncdf_var('e1u',name,Grid%e1u,strt,cnt)
-    CALL Read_Ncdf_var('e1v',name,Grid%e1v,strt,cnt)
-    CALL Read_Ncdf_var('e1f',name,Grid%e1f,strt,cnt)
-    CALL Read_Ncdf_var('e2t',name,Grid%e2t,strt,cnt)
-    CALL Read_Ncdf_var('e2u',name,Grid%e2u,strt,cnt)
-    CALL Read_Ncdf_var('e2v',name,Grid%e2v,strt,cnt)
-    CALL Read_Ncdf_var('e2f',name,Grid%e2f,strt,cnt)
-    CALL Read_Ncdf_var('nav_lon',name,Grid%nav_lon,strt,cnt)
-    CALL Read_Ncdf_var('nav_lat',name,Grid%nav_lat,strt,cnt)       
+    CALL read_ncdf_var('glamt',name,Grid%glamt,strt,cnt)
+    CALL read_ncdf_var('glamu',name,Grid%glamu,strt,cnt)
+    CALL read_ncdf_var('glamv',name,Grid%glamv,strt,cnt)
+    CALL read_ncdf_var('glamf',name,Grid%glamf,strt,cnt)
+    CALL read_ncdf_var('gphit',name,Grid%gphit,strt,cnt)
+    CALL read_ncdf_var('gphiu',name,Grid%gphiu,strt,cnt)
+    CALL read_ncdf_var('gphiv',name,Grid%gphiv,strt,cnt)
+    CALL read_ncdf_var('gphif',name,Grid%gphif,strt,cnt)
+    CALL read_ncdf_var('e1t',name,Grid%e1t,strt,cnt)
+    CALL read_ncdf_var('e1u',name,Grid%e1u,strt,cnt)
+    CALL read_ncdf_var('e1v',name,Grid%e1v,strt,cnt)
+    CALL read_ncdf_var('e1f',name,Grid%e1f,strt,cnt)
+    CALL read_ncdf_var('e2t',name,Grid%e2t,strt,cnt)
+    CALL read_ncdf_var('e2u',name,Grid%e2u,strt,cnt)
+    CALL read_ncdf_var('e2v',name,Grid%e2v,strt,cnt)
+    CALL read_ncdf_var('e2f',name,Grid%e2f,strt,cnt)
+    CALL read_ncdf_var('nav_lon',name,Grid%nav_lon,strt,cnt)
+    CALL read_ncdf_var('nav_lat',name,Grid%nav_lat,strt,cnt)       
     !
     WRITE(*,*) ' '
@@ -146,45 +146,45 @@
     !           
     dimnames = (/ 'x','y' /)
-    CALL Write_Ncdf_dim(dimnames(1),name,nxfin)
-    CALL Write_Ncdf_dim(dimnames(2),name,nyfin)
-    !      
-    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('glamt',dimnames,name,Grid%glamt,'double')
-    CALL Write_Ncdf_var('glamu',dimnames,name,Grid%glamu,'double')
-    CALL Write_Ncdf_var('glamv',dimnames,name,Grid%glamv,'double')
-    CALL Write_Ncdf_var('glamf',dimnames,name,Grid%glamf,'double')
-    CALL Write_Ncdf_var('gphit',dimnames,name,Grid%gphit,'double')
-    CALL Write_Ncdf_var('gphiu',dimnames,name,Grid%gphiu,'double')
-    CALL Write_Ncdf_var('gphiv',dimnames,name,Grid%gphiv,'double')
-    CALL Write_Ncdf_var('gphif',dimnames,name,Grid%gphif,'double')      
-    CALL Write_Ncdf_var('e1t',dimnames,name,Grid%e1t,'double')      
-    CALL Write_Ncdf_var('e1u',dimnames,name,Grid%e1u,'double')     
-    CALL Write_Ncdf_var('e1v',dimnames,name,Grid%e1v,'double')      
-    CALL Write_Ncdf_var('e1f',dimnames,name,Grid%e1f,'double')
-    CALL Write_Ncdf_var('e2t',dimnames,name,Grid%e2t,'double')
-    CALL Write_Ncdf_var('e2u',dimnames,name,Grid%e2u,'double')
-    CALL Write_Ncdf_var('e2v',dimnames,name,Grid%e2v,'double')
-    CALL Write_Ncdf_var('e2f',dimnames,name,Grid%e2f,'double')
-    !      
-    CALL Copy_Ncdf_att('nav_lon',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon))
-    CALL Copy_Ncdf_att('nav_lat',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))
-    CALL Copy_Ncdf_att('glamt',TRIM(parent_coordinate_file),name)
-    CALL Copy_Ncdf_att('glamu',TRIM(parent_coordinate_file),name)
-    CALL Copy_Ncdf_att('glamv',TRIM(parent_coordinate_file),name)
-    CALL Copy_Ncdf_att('glamf',TRIM(parent_coordinate_file),name)
-    CALL Copy_Ncdf_att('gphit',TRIM(parent_coordinate_file),name)
-    CALL Copy_Ncdf_att('gphiu',TRIM(parent_coordinate_file),name)
-    CALL Copy_Ncdf_att('gphiv',TRIM(parent_coordinate_file),name)
-    CALL Copy_Ncdf_att('gphif',TRIM(parent_coordinate_file),name)
-    CALL Copy_Ncdf_att('e1t',TRIM(parent_coordinate_file),name)
-    CALL Copy_Ncdf_att('e1u',TRIM(parent_coordinate_file),name)
-    CALL Copy_Ncdf_att('e1v',TRIM(parent_coordinate_file),name)
-    CALL Copy_Ncdf_att('e1f',TRIM(parent_coordinate_file),name)
-    CALL Copy_Ncdf_att('e2t',TRIM(parent_coordinate_file),name)
-    CALL Copy_Ncdf_att('e2u',TRIM(parent_coordinate_file),name)
-    CALL Copy_Ncdf_att('e2v',TRIM(parent_coordinate_file),name)
-    CALL Copy_Ncdf_att('e2f',TRIM(parent_coordinate_file),name)            
+    CALL write_ncdf_dim(dimnames(1),name,nxfin)
+    CALL write_ncdf_dim(dimnames(2),name,nyfin)
+    !      
+    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('glamt',dimnames,name,Grid%glamt,'double')
+    CALL write_ncdf_var('glamu',dimnames,name,Grid%glamu,'double')
+    CALL write_ncdf_var('glamv',dimnames,name,Grid%glamv,'double')
+    CALL write_ncdf_var('glamf',dimnames,name,Grid%glamf,'double')
+    CALL write_ncdf_var('gphit',dimnames,name,Grid%gphit,'double')
+    CALL write_ncdf_var('gphiu',dimnames,name,Grid%gphiu,'double')
+    CALL write_ncdf_var('gphiv',dimnames,name,Grid%gphiv,'double')
+    CALL write_ncdf_var('gphif',dimnames,name,Grid%gphif,'double')      
+    CALL write_ncdf_var('e1t',dimnames,name,Grid%e1t,'double')      
+    CALL write_ncdf_var('e1u',dimnames,name,Grid%e1u,'double')     
+    CALL write_ncdf_var('e1v',dimnames,name,Grid%e1v,'double')      
+    CALL write_ncdf_var('e1f',dimnames,name,Grid%e1f,'double')
+    CALL write_ncdf_var('e2t',dimnames,name,Grid%e2t,'double')
+    CALL write_ncdf_var('e2u',dimnames,name,Grid%e2u,'double')
+    CALL write_ncdf_var('e2v',dimnames,name,Grid%e2v,'double')
+    CALL write_ncdf_var('e2f',dimnames,name,Grid%e2f,'double')
+    !      
+    CALL copy_ncdf_att('nav_lon',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon))
+    CALL copy_ncdf_att('nav_lat',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))
+    CALL copy_ncdf_att('glamt',TRIM(parent_coordinate_file),name)
+    CALL copy_ncdf_att('glamu',TRIM(parent_coordinate_file),name)
+    CALL copy_ncdf_att('glamv',TRIM(parent_coordinate_file),name)
+    CALL copy_ncdf_att('glamf',TRIM(parent_coordinate_file),name)
+    CALL copy_ncdf_att('gphit',TRIM(parent_coordinate_file),name)
+    CALL copy_ncdf_att('gphiu',TRIM(parent_coordinate_file),name)
+    CALL copy_ncdf_att('gphiv',TRIM(parent_coordinate_file),name)
+    CALL copy_ncdf_att('gphif',TRIM(parent_coordinate_file),name)
+    CALL copy_ncdf_att('e1t',TRIM(parent_coordinate_file),name)
+    CALL copy_ncdf_att('e1u',TRIM(parent_coordinate_file),name)
+    CALL copy_ncdf_att('e1v',TRIM(parent_coordinate_file),name)
+    CALL copy_ncdf_att('e1f',TRIM(parent_coordinate_file),name)
+    CALL copy_ncdf_att('e2t',TRIM(parent_coordinate_file),name)
+    CALL copy_ncdf_att('e2u',TRIM(parent_coordinate_file),name)
+    CALL copy_ncdf_att('e2v',TRIM(parent_coordinate_file),name)
+    CALL copy_ncdf_att('e2f',TRIM(parent_coordinate_file),name)            
     !
     WRITE(*,*) ' '
@@ -209,5 +209,5 @@
     TYPE(Coordinates) :: Grid
     !
-    CALL Read_Ncdf_var('mbathy',name,Grid%Bathy_level)    
+    CALL read_ncdf_var('mbathy',name,Grid%Bathy_level)    
     !
     WRITE(*,*) ' '
@@ -232,18 +232,18 @@
     CHARACTER(len=1),DIMENSION(2) :: dimnames
     !
-    status = nf90_create(name,NF90_NOCLOBBER,ncid)
+    status = nf90_create(name,NF90_WRITE,ncid)
     status = nf90_close(ncid)          
     !
     dimnames = (/ 'x','y' /)
-    CALL Write_Ncdf_dim(dimnames(1),name,nxfin)
-    CALL Write_Ncdf_dim(dimnames(2),name,nyfin)
-    !      
-    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('mbathy' ,dimnames,name,Grid%bathy_level,'float')
-    !
-    CALL Copy_Ncdf_att('nav_lon',TRIM(parent_meshmask_file),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon))
-    CALL Copy_Ncdf_att('nav_lat',TRIM(parent_meshmask_file),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))
-    CALL Copy_Ncdf_att('mbathy' ,TRIM(parent_meshmask_file),name)       
+    CALL write_ncdf_dim(dimnames(1),name,nxfin)
+    CALL write_ncdf_dim(dimnames(2),name,nyfin)
+    !      
+    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('mbathy' ,dimnames,name,Grid%bathy_level,'float')
+    !
+    CALL copy_ncdf_att('nav_lon',TRIM(parent_meshmask_file),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon))
+    CALL copy_ncdf_att('nav_lat',TRIM(parent_meshmask_file),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))
+    CALL copy_ncdf_att('mbathy' ,TRIM(parent_meshmask_file),name)       
     !
     WRITE(*,*) ' '
@@ -276,13 +276,13 @@
        WRITE(*,*) ' etopo format for external high resolution database  '
        WRITE(*,*) '****'
-       CALL Read_Ncdf_var('lon',name,topo_lon)
-       CALL Read_Ncdf_var('lat',name,topo_lat)
+       CALL read_ncdf_var('lon',name,topo_lon)
+       CALL read_ncdf_var('lat',name,topo_lat)
     ELSE IF( Dims_Existence('x',name) .AND. Dims_Existence('y',name) ) THEN
        WRITE(*,*) '****'
        WRITE(*,*) ' OPA format for external high resolution database  '
        WRITE(*,*) '****'
-       CALL Read_Ncdf_var('nav_lon',name,CoarseGrid%nav_lon)
-       CALL Read_Ncdf_var('nav_lat',name,CoarseGrid%nav_lat)
-       CALL Read_Ncdf_var(parent_batmet_name,name,CoarseGrid%Bathy_meter)
+       CALL read_ncdf_var('nav_lon',name,CoarseGrid%nav_lon)
+       CALL read_ncdf_var('nav_lat',name,CoarseGrid%nav_lat)
+       CALL read_ncdf_var(parent_batmet_name,name,CoarseGrid%Bathy_meter)
        !            
        IF ( PRESENT(Pacifique) ) THEN
@@ -303,9 +303,7 @@
     ENDIF
     !
-    IF( MAXVAL(ChildGrid%glamt) > 180 ) THEN                 
+    IF( MAXVAL(ChildGrid%glamt) > 180. ) THEN                 
        !          
-       WHERE( topo_lon < 0 )
-          topo_lon = topo_lon + 360.
-       END WHERE
+       WHERE( topo_lon < 0. )   topo_lon = topo_lon + 360.
        !          
        i_min = MAXLOC(topo_lon,mask = topo_lon < MINVAL(ChildGrid%nav_lon)-zdel)
@@ -316,7 +314,9 @@
        tabdim1 = ( SIZE(topo_lon) - i_min(1) + 1 ) + i_max(1)                    
        !          
-       IF(j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(topo_lat,1) ) THEN
-          j_min(1) = j_min(1)-2
-          j_max(1) = j_max(1)+3
+       IF( ln_agrif_domain ) THEN
+          IF(j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(topo_lat,1) ) THEN
+             j_min(1) = j_min(1)-2
+             j_max(1) = j_max(1)+3
+          ENDIF
        ENDIF
        tabdim2 = j_max(1) - j_min(1) + 1
@@ -352,4 +352,6 @@
     ELSE
        !
+       WHERE( topo_lon > 180. )   topo_lon = topo_lon - 360.
+       !
        i_min = MAXLOC(topo_lon,mask = topo_lon < MINVAL(ChildGrid%nav_lon)-zdel)
        i_max = MINLOC(topo_lon,mask = topo_lon > MAXVAL(ChildGrid%nav_lon)+zdel)
@@ -357,13 +359,17 @@
        j_max = MINLOC(topo_lat,mask = topo_lat > MAXVAL(ChildGrid%nav_lat)+zdel)
        !      
-       IF(i_min(1)-2 >= 1 .AND. i_max(1)+3 <= SIZE(topo_lon,1) ) THEN
-          i_min(1) = i_min(1)-2
-          i_max(1) = i_max(1)+3
+       IF( ln_agrif_domain ) THEN
+          IF(i_min(1)-2 >= 1 .AND. i_max(1)+3 <= SIZE(topo_lon,1) ) THEN
+             i_min(1) = i_min(1)-2
+             i_max(1) = i_max(1)+3
+          ENDIF
        ENDIF
        tabdim1 = i_max(1) - i_min(1) + 1
        !
-       IF(j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(topo_lat,1) ) THEN
-          j_min(1) = j_min(1)-2
-          j_max(1) = j_max(1)+3
+       IF( ln_agrif_domain ) THEN
+          IF(j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(topo_lat,1) ) THEN
+             j_min(1) = j_min(1)-2
+             j_max(1) = j_max(1)+3
+          ENDIF
        ENDIF
        tabdim2 = j_max(1) - j_min(1) + 1
@@ -423,5 +429,5 @@
     TYPE(Coordinates) :: Grid
     !
-    CALL Read_Ncdf_var(parent_batmet_name,name,Grid%Bathy_meter)    
+    CALL read_ncdf_var(parent_batmet_name,name,Grid%Bathy_meter)    
     !
     WRITE(*,*) ' '
@@ -454,15 +460,15 @@
     dimnames = (/ 'x','y' /)
 
-    CALL Write_Ncdf_dim(dimnames(1),name,nx)
-    CALL Write_Ncdf_dim(dimnames(2),name,ny)
+    CALL write_ncdf_dim(dimnames(1),name,nx)
+    CALL write_ncdf_dim(dimnames(2),name,ny)
     !     
-    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_batmet_name,dimnames,name,Grid%bathy_meter,'float')
-    CALL Write_Ncdf_var('weight'          ,dimnames,name,Grid%wgt        ,'float')
-    !
-    CALL Copy_Ncdf_att('nav_lon'         ,TRIM(parent_bathy_meter),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon))
-    CALL Copy_Ncdf_att('nav_lat'         ,TRIM(parent_bathy_meter),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))  
-    CALL Copy_Ncdf_att(parent_batmet_name,TRIM(parent_bathy_meter),name)
+    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_batmet_name,dimnames,name,Grid%bathy_meter,'float')
+    CALL write_ncdf_var('weight'          ,dimnames,name,Grid%wgt        ,'float')
+    !
+    CALL copy_ncdf_att('nav_lon'         ,TRIM(parent_bathy_meter),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon))
+    CALL copy_ncdf_att('nav_lat'         ,TRIM(parent_bathy_meter),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))  
+    CALL copy_ncdf_att(parent_batmet_name,TRIM(parent_bathy_meter),name)
     !
     WRITE(*,*) ' '
@@ -474,4 +480,298 @@
   END FUNCTION Write_Bathy_meter
   !
+  !*****************************************************
+  !   function write_domcfg(name,Grid)
+  !*****************************************************
+
+  INTEGER FUNCTION write_domcfg(name,Grid)
+     !-----------------------------------------
+     ! It creates a domain_cfg.nc used in NEMO4
+     !-----------------------------------------    
+     !
+     USE io_netcdf
+     !
+     CHARACTER(*) name
+     TYPE(Coordinates) :: Grid
+     !
+     INTEGER :: status, ncid
+     INTEGER :: nx, ny, jk
+     INTEGER :: ln_sco, ln_isfcav, ln_zco, ln_zps, jperio
+     REAL*8  :: rpi, rad, rday, rsiyea, rsiday, omega
+     !
+     CHARACTER(len=1),     DIMENSION(3)     ::   dimnames
+     REAL*8 ,              DIMENSION(N)     ::   e3t_1d, e3w_1d, gdept_1d
+     REAL*8 , ALLOCATABLE, DIMENSION(:,:)   ::   ff_t, ff_f
+     INTEGER, ALLOCATABLE, DIMENSION(:,:)   ::   bottom_level, top_level
+     REAL*8 , ALLOCATABLE, DIMENSION(:,:,:) ::   e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0
+     !
+     ! size of the Grid
+     nx = SIZE(Grid%bathy_meter,1)
+     ny = SIZE(Grid%bathy_meter,2)
+
+     ! allocate needed arrays for domain_cfg
+     ALLOCATE( ff_t(nx,ny), ff_f(nx,ny) )
+     ALLOCATE( bottom_level(nx,ny), top_level(nx,ny) )
+     ALLOCATE( e3t_0(nx,ny,N), e3u_0 (nx,ny,N), e3v_0 (nx,ny,N), e3f_0(nx,ny,N),  &
+        &      e3w_0(nx,ny,N), e3uw_0(nx,ny,N), e3vw_0(nx,ny,N) )
+
+     ! some physical parameters
+     rpi    = 3.141592653589793
+     rad    = 3.141592653589793 / 180.
+     rday   = 24.*60.*60.
+     rsiyea = 365.25 * rday * 2. * rpi / 6.283076
+     rsiday = rday / ( 1. + rday / rsiyea )
+     omega  = 2. * rpi / rsiday
+
+     ! Coriolis
+     ff_f(:,:) = 2. * omega * SIN( rad * Grid%gphif(:,:) )     ! compute it on the sphere at f-point
+     ff_t(:,:) = 2. * omega * SIN( rad * Grid%gphit(:,:) )     !    -        -       -    at t-point
+
+     ! top/bottom levels
+     bottom_level(:,:) = Grid%bathy_level(:,:)
+     top_level(:,:)    = MIN( 1, bottom_level(:,:) )
+
+     ! vertical scale factors
+     CALL zgr_z( e3t_1d, e3w_1d, gdept_1d )    
+     DO jk = 1, N
+        e3t_0  (:,:,jk) = e3t_1d  (jk)
+        e3u_0  (:,:,jk) = e3t_1d  (jk)
+        e3v_0  (:,:,jk) = e3t_1d  (jk)
+        e3f_0  (:,:,jk) = e3t_1d  (jk)
+        e3w_0  (:,:,jk) = e3w_1d  (jk)
+        e3uw_0 (:,:,jk) = e3w_1d  (jk)
+        e3vw_0 (:,:,jk) = e3w_1d  (jk)
+     END DO
+
+     ! logicals and others
+     ln_sco = 0
+     ln_isfcav = 0
+     IF( partial_steps ) THEN
+        ln_zps = 1
+        ln_zco = 0
+     ELSE
+        ln_zps = 0
+        ln_zco = 1
+     ENDIF
+
+     ! closed domain (agrif)
+     jperio = 0
+
+     ! -------------------
+     ! write domain_cfg.nc
+     ! -------------------
+     status = nf90_create(name,NF90_WRITE,ncid)
+     status = nf90_close(ncid)
+     !
+     ! dimensions
+     dimnames = (/'x','y','z'/)
+     CALL write_ncdf_dim(dimnames(1),name,nx)
+     CALL write_ncdf_dim(dimnames(2),name,ny)
+     CALL write_ncdf_dim(dimnames(3),name,N)
+     !      
+     ! variables
+     CALL write_ncdf_var('nav_lon',dimnames(1:2),name,Grid%nav_lon,'float')      
+     CALL write_ncdf_var('nav_lat',dimnames(1:2),name,Grid%nav_lat,'float')
+     CALL write_ncdf_var('nav_lev',dimnames(3)  ,name,gdept_1d    ,'float')
+     !
+     CALL write_ncdf_var('jpiglo',name,nx    ,'integer')
+     CALL write_ncdf_var('jpjglo',name,ny    ,'integer')
+     CALL write_ncdf_var('jpkglo',name,N     ,'integer')
+     CALL write_ncdf_var('jperio',name,jperio,'integer')
+     !
+     CALL write_ncdf_var('ln_zco'   ,name,ln_zco   ,'integer')
+     CALL write_ncdf_var('ln_zps'   ,name,ln_zps   ,'integer')
+     CALL write_ncdf_var('ln_sco'   ,name,ln_sco   ,'integer')
+     CALL write_ncdf_var('ln_isfcav',name,ln_isfcav,'integer')
+
+     CALL write_ncdf_var('glamt',dimnames(1:2),name,Grid%glamt,'double')
+     CALL write_ncdf_var('glamu',dimnames(1:2),name,Grid%glamu,'double')
+     CALL write_ncdf_var('glamv',dimnames(1:2),name,Grid%glamv,'double')
+     CALL write_ncdf_var('glamf',dimnames(1:2),name,Grid%glamf,'double')
+     CALL write_ncdf_var('gphit',dimnames(1:2),name,Grid%gphit,'double')
+     CALL write_ncdf_var('gphiu',dimnames(1:2),name,Grid%gphiu,'double')
+     CALL write_ncdf_var('gphiv',dimnames(1:2),name,Grid%gphiv,'double')
+     CALL write_ncdf_var('gphif',dimnames(1:2),name,Grid%gphif,'double')      
+
+     CALL write_ncdf_var('e1t',dimnames(1:2),name,Grid%e1t,'double')      
+     CALL write_ncdf_var('e1u',dimnames(1:2),name,Grid%e1u,'double')     
+     CALL write_ncdf_var('e1v',dimnames(1:2),name,Grid%e1v,'double')      
+     CALL write_ncdf_var('e1f',dimnames(1:2),name,Grid%e1f,'double')
+     CALL write_ncdf_var('e2t',dimnames(1:2),name,Grid%e2t,'double')
+     CALL write_ncdf_var('e2u',dimnames(1:2),name,Grid%e2u,'double')
+     CALL write_ncdf_var('e2v',dimnames(1:2),name,Grid%e2v,'double')
+     CALL write_ncdf_var('e2f',dimnames(1:2),name,Grid%e2f,'double')
+
+     CALL write_ncdf_var('ff_f',dimnames(1:2),name,ff_f,'double')
+     CALL write_ncdf_var('ff_t',dimnames(1:2),name,ff_t,'double')
+
+     CALL write_ncdf_var('e3t_1d',dimnames(3),name,e3t_1d,'double')
+     CALL write_ncdf_var('e3w_1d',dimnames(3),name,e3w_1d,'double')
+
+     CALL write_ncdf_var('e3t_0' ,dimnames(:),name,e3t_0 ,'double')
+     CALL write_ncdf_var('e3w_0' ,dimnames(:),name,e3w_0 ,'double')
+     CALL write_ncdf_var('e3u_0' ,dimnames(:),name,e3u_0 ,'double')
+     CALL write_ncdf_var('e3v_0' ,dimnames(:),name,e3v_0 ,'double')
+     CALL write_ncdf_var('e3f_0' ,dimnames(:),name,e3f_0 ,'double')
+     CALL write_ncdf_var('e3w_0' ,dimnames(:),name,e3w_0 ,'double')
+     CALL write_ncdf_var('e3uw_0',dimnames(:),name,e3uw_0,'double')
+     CALL write_ncdf_var('e3vw_0',dimnames(:),name,e3vw_0,'double')
+
+     CALL write_ncdf_var('bottom_level',dimnames(1:2),name,bottom_level,'integer')
+     CALL write_ncdf_var('top_level'   ,dimnames(1:2),name,top_level   ,'integer')
+
+     CALL write_ncdf_var('bathy_meter' ,dimnames(1:2),name,Grid%bathy_meter,'float')
+
+     ! some attributes      
+     CALL copy_ncdf_att('nav_lon',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon))
+     CALL copy_ncdf_att('nav_lat',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))
+
+     CALL copy_ncdf_att('glamt',TRIM(parent_coordinate_file),name)
+     CALL copy_ncdf_att('glamu',TRIM(parent_coordinate_file),name)
+     CALL copy_ncdf_att('glamv',TRIM(parent_coordinate_file),name)
+     CALL copy_ncdf_att('glamf',TRIM(parent_coordinate_file),name)
+     CALL copy_ncdf_att('gphit',TRIM(parent_coordinate_file),name)
+     CALL copy_ncdf_att('gphiu',TRIM(parent_coordinate_file),name)
+     CALL copy_ncdf_att('gphiv',TRIM(parent_coordinate_file),name)
+     CALL copy_ncdf_att('gphif',TRIM(parent_coordinate_file),name)
+
+     CALL copy_ncdf_att('e1t',TRIM(parent_coordinate_file),name)
+     CALL copy_ncdf_att('e1u',TRIM(parent_coordinate_file),name)
+     CALL copy_ncdf_att('e1v',TRIM(parent_coordinate_file),name)
+     CALL copy_ncdf_att('e1f',TRIM(parent_coordinate_file),name)
+     CALL copy_ncdf_att('e2t',TRIM(parent_coordinate_file),name)
+     CALL copy_ncdf_att('e2u',TRIM(parent_coordinate_file),name)
+     CALL copy_ncdf_att('e2v',TRIM(parent_coordinate_file),name)
+     CALL copy_ncdf_att('e2f',TRIM(parent_coordinate_file),name)            
+     !
+     ! control print
+     WRITE(*,*) ' '
+     WRITE(*,*) 'Writing domcfg file: ',name
+     WRITE(*,*) ' '
+     !
+     DEALLOCATE( ff_t, ff_f )
+     DEALLOCATE( bottom_level, top_level )
+     DEALLOCATE( e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 )
+     !      
+     write_domcfg = 1
+
+  END FUNCTION write_domcfg
+  !
+  SUBROUTINE zgr_z(e3t_1d, e3w_1d, gdept_1d)
+     !!----------------------------------------------------------------------
+     !!                   ***  ROUTINE zgr_z (from NEMO4) ***
+     !!                   
+     !! ** Purpose :   set the depth of model levels and the resulting 
+     !!      vertical scale factors.
+     !!
+     !! ** Method  :   z-coordinate system (use in all type of coordinate)
+     !!        The depth of model levels is defined from an analytical
+     !!      function the derivative of which gives the scale factors.
+     !!        both depth and scale factors only depend on k (1d arrays).
+     !!              w-level: gdepw_1d  = gdep(k)
+     !!                       e3w_1d(k) = dk(gdep)(k)     = e3(k)
+     !!              t-level: gdept_1d  = gdep(k+0.5)
+     !!                       e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5)
+     !!
+     !! ** Action  : - gdept_1d, gdepw_1d : depth of T- and W-point (m)
+     !!              - e3t_1d  , e3w_1d   : scale factors at T- and W-levels (m)
+     !!
+     !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
+     !!----------------------------------------------------------------------
+     INTEGER  ::   jk                   ! dummy loop indices
+     INTEGER  ::   jpk
+     REAL*8 ::   zt, zw                 ! temporary scalars
+     REAL*8 ::   zsur, za0, za1, zkth   ! Values set from parameters in
+     REAL*8 ::   zacr, zdzmin, zhmax    ! par_CONFIG_Rxx.h90
+     REAL*8 ::   za2, zkth2, zacr2      ! Values for optional double tanh function set from parameters 
+
+     REAL*8, DIMENSION(:), INTENT(inout) ::  e3t_1d, e3w_1d, gdept_1d
+     REAL*8, DIMENSION(N)                ::  gdepw_1d
+     !!----------------------------------------------------------------------
+     !
+     !
+     ! Set variables from parameters
+     ! ------------------------------
+     zkth = ppkth       ;   zacr = ppacr
+     zdzmin = ppdzmin   ;   zhmax = pphmax
+     zkth2 = ppkth2     ;   zacr2 = ppacr2   ! optional (ldbletanh=T) double tanh parameters
+
+     ! If pa1 and pa0 and psur are et to pp_to_be_computed
+     !  za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
+    IF ( ( pa0 == 0 .OR. pa1 == 0 .OR. psur == 0 ) &
+         .AND. ppdzmin.NE.0 .AND. pphmax.NE.0 ) THEN
+        !
+        za1  = (  ppdzmin - pphmax / FLOAT(N-1)  )                                                      &
+           & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(N-1) * (  LOG( COSH( (N - ppkth) / ppacr) )      &
+           &                                                   - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
+        za0  = ppdzmin - za1 *              TANH( (1-ppkth) / ppacr )
+        zsur =   - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
+     ELSE
+        za1 = pa1 ;       za0 = pa0 ;          zsur = psur
+        za2 = pa2                            ! optional (ldbletanh=T) double tanh parameter
+     ENDIF
+
+     ! Reference z-coordinate (depth - scale factor at T- and W-points)
+     ! ======================
+     IF( ppkth == 0. ) THEN            !  uniform vertical grid 
+
+        za1 = zhmax / FLOAT(N-1)
+
+        DO jk = 1, N
+           zw = FLOAT( jk )
+           zt = FLOAT( jk ) + 0.5
+           gdepw_1d(jk) = ( zw - 1 ) * za1
+           gdept_1d(jk) = ( zt - 1 ) * za1
+           e3w_1d  (jk) =  za1
+           e3t_1d  (jk) =  za1
+        END DO
+     ELSE                                ! Madec & Imbard 1996 function
+        IF( .NOT. ldbletanh ) THEN
+           DO jk = 1, N
+              zw = REAL( jk )
+              zt = REAL( jk ) + 0.5
+              gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )  )
+              gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )  )
+              e3w_1d  (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   )
+              e3t_1d  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   )
+           END DO
+        ELSE
+           DO jk = 1, N
+              zw = FLOAT( jk )
+              zt = FLOAT( jk ) + 0.5
+              ! Double tanh function
+              gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr  ) )    &
+                 &                             + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) )  )
+              gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr  ) )    &
+                 &                             + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) )  )
+              e3w_1d  (jk) =          za0      + za1        * TANH(       (zw-zkth ) / zacr  )      &
+                 &                             + za2        * TANH(       (zw-zkth2) / zacr2 )
+              e3t_1d  (jk) =          za0      + za1        * TANH(       (zt-zkth ) / zacr  )      &
+                 &                             + za2        * TANH(       (zt-zkth2) / zacr2 )
+           END DO
+        ENDIF
+        gdepw_1d(1) = 0.                    ! force first w-level to be exactly at zero
+     ENDIF
+
+     IF ( ln_e3_dep ) THEN      ! e3. = dk[gdep]   
+        !
+        !==>>>   need to be like this to compute the pressure gradient with ISF. 
+        !        If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
+        !        define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively
+        !
+        DO jk = 1, N-1
+           e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)
+        END DO
+        e3t_1d(N) = e3t_1d(N-1)   ! we don't care because this level is masked in NEMO
+
+        DO jk = 2, N
+           e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)
+        END DO
+        e3w_1d(1  ) = 2. * (gdept_1d(1) - gdepw_1d(1))
+     END IF
+
+     !
+  END SUBROUTINE zgr_z
+
+
   !*****************************************************
   !   function set_child_name(Parentname,Childname)
@@ -501,8 +801,4 @@
     !   
   END SUBROUTINE set_child_name
-  !
-  !*****************************************************
-  !   function set_child_name(Parentname,Childname)
-  !*****************************************************
   !
   !*****************************************************
@@ -558,8 +854,4 @@
   END SUBROUTINE get_interptype
   !      
-  !*****************************************************              
-  !   end subroutine get_interptype
-  !*****************************************************
-  !           
   !*****************************************************
   !   subroutine Init_mask(name,Grid)
@@ -576,7 +868,7 @@
     !
     IF(jpiglo == 1 .AND. jpjglo == 1) THEN
-       CALL Read_Ncdf_var('Bathy_level',name,Grid%Bathy_level)
+       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level)
     ELSE
-       CALL Read_Ncdf_var('Bathy_level',name,Grid%Bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) )
+       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) )
     ENDIF
 
@@ -655,8 +947,4 @@
   !
   !*****************************************************
-  !   end subroutine Init_mask
-  !*****************************************************
-  ! 
-  !*****************************************************
   !   subroutine Init_Tmask(name,Grid)
   !*****************************************************
@@ -672,7 +960,7 @@
     !
     IF(jpiglo == 1 .AND. jpjglo == 1) THEN
-       CALL Read_Ncdf_var('Bathy_level',name,Grid%Bathy_level)
+       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level)
     ELSE
-       CALL Read_Ncdf_var('Bathy_level',name,Grid%Bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) )
+       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) )
     ENDIF
     !
@@ -712,5 +1000,5 @@
     TYPE(Coordinates) :: Grid
     !
-    CALL Read_Ncdf_var('Bathy_level',filename,Grid%Bathy_level)
+    CALL read_ncdf_var('Bathy_level',filename,Grid%Bathy_level)
     ! 
     nx = SIZE(Grid%Bathy_level,1)
@@ -764,8 +1052,4 @@
   END SUBROUTINE get_mask
   !
-  !*****************************************************
-  !   end subroutine get_mask
-  !*****************************************************
-  !       
   !
   !*****************************************************
@@ -803,5 +1087,4 @@
     !
   END SUBROUTINE write_dimg_var
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!             
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
 END MODULE agrif_readwrite
Index: /utils/tools/NESTING/src/agrif_types.f90
===================================================================
--- /utils/tools/NESTING/src/agrif_types.f90	(revision 10024)
+++ /utils/tools/NESTING/src/agrif_types.f90	(revision 10025)
@@ -15,17 +15,17 @@
   TYPE Coordinates
      !
-     REAL*8, DIMENSION(:,:), POINTER :: nav_lon,nav_lat => NULL()
-     REAL*8, DIMENSION(:,:), POINTER :: glamv, glamu, glamt, glamf => NULL()
-     REAL*8, DIMENSION(:,:), POINTER :: gphit, gphiu, gphiv, gphif => NULL()
-     REAL*8, DIMENSION(:,:), POINTER :: e1t, e1u, e1v, e1f => NULL()
-     REAL*8, DIMENSION(:,:), POINTER :: e2t, e2u, e2v, e2f => NULL()
-     REAL*8, DIMENSION(:,:), POINTER :: bathy_level => NULL()
-     REAL*8, DIMENSION(:,:), POINTER :: bathy_meter => NULL()
-     REAL*8, DIMENSION(:,:), POINTER :: wgt => NULL()
-     REAL*8, DIMENSION(:,:,:),POINTER :: fmask,umask,vmask,tmask => NULL()
-     REAL*8, DIMENSION(:,:,:),POINTER :: e3t_ps,e3w_ps,gdept_ps,gdepwps => NULL()
-     REAL*8, DIMENSION(:,:),POINTER :: gdepw_ps => NULL()
-     REAL*8, DIMENSION(:), POINTER :: gdeptht => NULL()
-     INTEGER, DIMENSION(:) , POINTER :: time_steps => NULL()
+     REAL*8,  DIMENSION(:,:),  POINTER :: nav_lon, nav_lat => NULL()
+     REAL*8,  DIMENSION(:,:),  POINTER :: glamv, glamu, glamt, glamf => NULL()
+     REAL*8,  DIMENSION(:,:),  POINTER :: gphit, gphiu, gphiv, gphif => NULL()
+     REAL*8,  DIMENSION(:,:),  POINTER :: e1t, e1u, e1v, e1f => NULL()
+     REAL*8,  DIMENSION(:,:),  POINTER :: e2t, e2u, e2v, e2f => NULL()
+     REAL*8,  DIMENSION(:,:),  POINTER :: bathy_level => NULL()
+     REAL*8,  DIMENSION(:,:),  POINTER :: bathy_meter => NULL()
+     REAL*8,  DIMENSION(:,:),  POINTER :: wgt => NULL()
+     REAL*8,  DIMENSION(:,:,:),POINTER :: fmask, umask, vmask, tmask => NULL()
+     REAL*8,  DIMENSION(:,:,:),POINTER :: e3t_ps, e3w_ps, gdept_ps, gdepwps => NULL()
+     REAL*8,  DIMENSION(:,:),  POINTER :: gdepw_ps => NULL()
+     REAL*8,  DIMENSION(:),    POINTER :: gdeptht => NULL()
+     INTEGER, DIMENSION(:) ,   POINTER :: time_steps => NULL()
      !     
   END TYPE Coordinates
@@ -33,5 +33,5 @@
   !
   !
-  CHARACTER*8,DIMENSION(10) :: flxtab = (/'socliot1','socliot2','socliopl',&
+  CHARACTER*8,DIMENSION(10) :: flxtab = (/'socliot1','socliot2','socliopl', &
        'socliocl','socliohu','socliowi','soshfldo','sohefldo','sowaflup','sofbt   '/)
   !
@@ -41,53 +41,52 @@
   !**************************************************************
   !
-  INTEGER irafx,irafy
-  INTEGER nxfin,nyfin
-  INTEGER, PARAMETER :: nbghostcellsfine = 3 
-  INTEGER, PARAMETER :: nbghostcellscoarse = 1
-  !      
-  INTEGER imin,jmin,imax,jmax,rho,rhot
-  INTEGER shlat
-  INTEGER N,type_bathy_interp
+  INTEGER ::   irafx, irafy
+  INTEGER ::   nxfin, nyfin
+  !      
+  INTEGER ::   nbghostcellsfine, imin, jmin, imax, jmax, rho, rhot
+  INTEGER ::   shlat
+  INTEGER ::   N, type_bathy_interp
   ! 
-  INTEGER jpizoom,jpjzoom,npt_connect,npt_copy
-  !      
-  REAL*8 rn_hmin
-  REAL*8 ppkth2, ppacr2, ppkth,ppacr,ppdzmin,pphmax,smoothing_factor,e3zps_min,e3zps_rat
-  REAL*8 psur,pa0,pa1,pa2,adatrj
+  INTEGER ::   jpizoom, jpjzoom, npt_connect, npt_copy
+  !      
+  REAL*8 ::   rn_hmin
+  REAL*8 ::   ppkth2, ppacr2, ppkth, ppacr, ppdzmin, pphmax, smoothing_factor, e3zps_min, e3zps_rat
+  REAL*8 ::   psur, pa0, pa1, pa2, adatrj
   !       
-  LOGICAL ldbletanh,ln_e3_dep
-  LOGICAL partial_steps,smoothing,bathy_update
-  LOGICAL new_topo,removeclosedseas,dimg,iom_activated
+  LOGICAL ::   ldbletanh, ln_e3_dep
+  LOGICAL ::   partial_steps, smoothing, bathy_update
+  LOGICAL ::   new_topo, removeclosedseas, dimg, iom_activated
+  LOGICAL ::   ln_agrif_domain
   !       
-  CHARACTER*100 parent_meshmask_file,elevation_database,parent_bathy_meter
-  CHARACTER*100 elevation_name,parent_batmet_name
-  CHARACTER*100 parent_coordinate_file,restart_file,updated_parent_file,restart_trc_file
-  CHARACTER*100 dimg_output_file,interp_type
-  !      
-  CHARACTER(len=80),DIMENSION(20) :: flx_Files, u_files, v_files
-  CHARACTER(len=255),DIMENSION(20) :: VAR_INTERP
+  CHARACTER*100 ::   parent_meshmask_file, elevation_database, parent_bathy_meter, parent_domcfg_output
+  CHARACTER*100 ::   elevation_name, parent_batmet_name
+  CHARACTER*100 ::   parent_coordinate_file, restart_file, updated_parent_file, updated_parent_domcfg, restart_trc_file
+  CHARACTER*100 ::   dimg_output_file, interp_type
+  !      
+  CHARACTER(len=80) , DIMENSION(20) :: flx_Files, u_files, v_files
+  CHARACTER(len=255), DIMENSION(20) :: VAR_INTERP
   !
   NAMELIST /input_output/iom_activated
   !
-  NAMELIST /coarse_grid_files/parent_coordinate_file,parent_meshmask_file
-  !      
-  NAMELIST /bathymetry/new_topo,elevation_database,elevation_name,smoothing,smoothing_factor, &
-       npt_connect,npt_copy,removeclosedseas,type_bathy_interp,rn_hmin      
-  !      
-  NAMELIST /nesting/imin,imax,jmin,jmax,rho,rhot,bathy_update,updated_parent_file      
-  !
-  NAMELIST /vertical_grid/ppkth,ppacr,ppdzmin,pphmax,psur,pa0,pa1,N,ldbletanh,ln_e3_dep,pa2,ppkth2,ppacr2
+  NAMELIST /coarse_grid_files/parent_coordinate_file, parent_meshmask_file, parent_domcfg_output
+  !      
+  NAMELIST /bathymetry/new_topo, elevation_database, elevation_name, smoothing, smoothing_factor,  &
+                       ln_agrif_domain, npt_connect, npt_copy, removeclosedseas, type_bathy_interp, rn_hmin      
+  !      
+  NAMELIST /nesting/nbghostcellsfine, imin, imax, jmin, jmax, rho, rhot, bathy_update, updated_parent_file, updated_parent_domcfg      
+  !
+  NAMELIST /vertical_grid/ppkth, ppacr, ppdzmin, pphmax, psur, pa0, pa1, N, ldbletanh, ln_e3_dep, pa2, ppkth2, ppacr2
   ! 
-  NAMELIST /partial_cells/partial_steps,parent_bathy_meter,parent_batmet_name,e3zps_min,e3zps_rat      
-  !
-  NAMELIST /nemo_coarse_grid/ jpizoom,jpjzoom 
+  NAMELIST /partial_cells/partial_steps, parent_bathy_meter, parent_batmet_name, e3zps_min, e3zps_rat      
+  !
+  NAMELIST /nemo_coarse_grid/ jpizoom, jpjzoom 
   !          
-  NAMELIST /forcing_files/ flx_files, u_files, v_files 
+  NAMELIST /forcing_files/ flx_files,  u_files,  v_files 
   !           
   NAMELIST /interp/ VAR_INTERP
   !      
-  NAMELIST /restart/ restart_file,shlat,dimg,dimg_output_file,adatrj,interp_type 
-
-  NAMELIST /restart_trc/ restart_trc_file,interp_type 
+  NAMELIST /restart/ restart_file, shlat, dimg, dimg_output_file, adatrj, interp_type 
+
+  NAMELIST /restart_trc/ restart_trc_file, interp_type 
   !
 CONTAINS
@@ -102,8 +101,8 @@
   !********************************************************
   !       
-  SUBROUTINE agrif_grid_allocate(Grid,nx,ny)
+  SUBROUTINE agrif_grid_allocate(Grid, nx, ny)
     !
     TYPE(Coordinates) :: Grid
-    INTEGER :: nx,ny
+    INTEGER :: nx, ny
     !
     ALLOCATE(Grid%nav_lon(nx,ny),Grid%nav_lat(nx,ny))
@@ -175,11 +174,17 @@
        jmax = jmax + jpjzoom - 1
        !
-       nxfin = (imax-imin)*irafx+2*nbghostcellsfine+2
-       nyfin = (jmax-jmin)*irafy+2*nbghostcellsfine+2
+       IF( ln_agrif_domain ) THEN
+          nxfin = (imax-imin)*irafx+2*nbghostcellsfine+2
+          nyfin = (jmax-jmin)*irafy+2*nbghostcellsfine+2
+       ELSE
+          nbghostcellsfine = 0
+          nxfin = (imax-imin+1)*irafx
+          nyfin = (jmax-jmin+1)*irafy
+       ENDIF
        !
     ELSE
        !
        PRINT *,'namelist file ''',TRIM(namelistname),''' not found'
-       STOP	 
+       STOP 
        !
     END IF
Index: /utils/tools/NESTING/src/io_netcdf.f90
===================================================================
--- /utils/tools/NESTING/src/io_netcdf.f90	(revision 10024)
+++ /utils/tools/NESTING/src/io_netcdf.f90	(revision 10025)
@@ -25,39 +25,21 @@
   USE agrif_types
   !      
-  INTERFACE Read_Ncdf_var
-     MODULE PROCEDURE Read_Ncdf_var1d_Real,   &
-          Read_Ncdf_var2d_Real,   &
-          Read_Ncdf_var2d_Real_bis,   &
-          Read_Ncdf_var3d_Real,   &
-          Read_Ncdf_var4d_Real,   &
-          Read_Ncdf_var3d_Real_t, &
-          Read_Ncdf_var4d_Real_t, &
-          Read_Ncdf_var4d_Real_nt, &
-          Read_Ncdf_var1d_Int,    &
-          Read_Ncdf_var2d_Int,    &
-          Read_Ncdf_var3d_Int,    &
-          Read_Ncdf_var4d_Int,    &
-          Read_Ncdf_var0d_Int,    &
-          Read_Ncdf_var0d_Real
+  INTERFACE read_ncdf_var
+     MODULE PROCEDURE  &
+        read_ncdf_var0d_real, read_ncdf_var1d_real, read_ncdf_var2d_real  , read_ncdf_var2d_real_bis,                        &
+        read_ncdf_var3d_real, read_ncdf_var4d_real, read_ncdf_var3d_real_t, read_ncdf_var4d_real_t, read_ncdf_var4d_real_nt, &
+        read_ncdf_var0d_int,  read_ncdf_var1d_int , read_ncdf_var2d_int   , read_ncdf_var3d_int   , read_ncdf_var4d_int
   END INTERFACE
   !
-  INTERFACE Write_Ncdf_var
-     MODULE PROCEDURE Write_Ncdf_var1d_Real,   &
-          Write_Ncdf_var2d_Real,   &
-          Write_Ncdf_var3d_Real,   &
-          Write_Ncdf_var4d_Real,   &
-          Write_Ncdf_var3d_Real_t, &
-          Write_Ncdf_var4d_Real_t, &
-          Write_Ncdf_var4d_Real_nt, &                           
-          Write_Ncdf_var2d_Real_bis,&			   
-          Write_Ncdf_var1d_Int,    &
-          Write_Ncdf_var2d_Int,    &
-          Write_Ncdf_var3d_Int,    &
-          Write_Ncdf_var4d_Int,    &
-          Write_Ncdf_var0d_Real
+  INTERFACE write_ncdf_var
+     MODULE PROCEDURE &
+        write_ncdf_var0d_real, write_ncdf_var1d_real  , write_ncdf_var2d_real  , write_ncdf_var3d_real,    &
+        write_ncdf_var4d_real, write_ncdf_var3d_real_t, write_ncdf_var4d_real_t, write_ncdf_var4d_real_nt, &                           
+        write_ncdf_var2d_real_bis ,                                                                        &			   
+        write_ncdf_var0d_int, write_ncdf_var1d_int, write_ncdf_var2d_int, write_ncdf_var3d_int, write_ncdf_var4d_int
   END INTERFACE
   ! 
-  INTERFACE Copy_Ncdf_att
-     MODULE PROCEDURE Copy_Ncdf_att_latlon,Copy_Ncdf_att_var
+  INTERFACE copy_ncdf_att
+     MODULE PROCEDURE copy_ncdf_att_latlon,copy_ncdf_att_var
   END INTERFACE
   !
@@ -65,5 +47,5 @@
   !
   !****************************************************************
-  !   subroutine Read_Ncdf_dim					*
+  !   subroutine read_ncdf_dim					*
   !								*
   ! subroutine to retrieve value of a given dimension 		*
@@ -75,5 +57,5 @@
   !****************************************************************
   !
-  SUBROUTINE Read_Ncdf_dim(dimname,file,dimval)   
+  SUBROUTINE read_ncdf_dim(dimname,file,dimval)   
     !      
     IMPLICIT NONE
@@ -97,12 +79,12 @@
     status = nf90_close(ncid)
     !      
-  END SUBROUTINE Read_Ncdf_dim
+  END SUBROUTINE read_ncdf_dim
   !
   !**************************************************************
-  ! end subroutine Read_Ncdf_dim
+  ! end subroutine read_ncdf_dim
   !**************************************************************      
   !
   !****************************************************************
-  !   subroutine Write_Ncdf_dim					*
+  !   subroutine write_ncdf_dim					*
   !								*
   ! subroutine to write a dimension in a given file 		*
@@ -114,5 +96,5 @@
   !****************************************************************
   !
-  SUBROUTINE Write_Ncdf_dim(dimname,file,dimval)   
+  SUBROUTINE write_ncdf_dim(dimname,file,dimval)   
     !      
     IMPLICIT NONE
@@ -143,12 +125,12 @@
     status = nf90_close(ncid)
     !      
-  END SUBROUTINE Write_Ncdf_dim
+  END SUBROUTINE write_ncdf_dim
   !
   !**************************************************************
-  ! end subroutine Write_Ncdf_dim
+  ! end subroutine write_ncdf_dim
   !**************************************************************            
   ! 
   !****************************************************************
-  !   subroutine Read_Ncdf_var					*
+  !   subroutine read_ncdf_var					*
   !								*
   ! subroutine to retrieve values of a given variable 		*
@@ -159,10 +141,6 @@
   !								*
   !****************************************************************
-  !          
-  !**************************************************************
-  !   subroutine Read_Ncdf_var1d_real
-  !**************************************************************          
-  !      
-  SUBROUTINE Read_Ncdf_var1d_Real(varname,file,tabvar)
+  !      
+  SUBROUTINE read_ncdf_var1d_real(varname,file,tabvar)
     !      
     IMPLICIT NONE
@@ -202,16 +180,8 @@
     status = nf90_close(ncid)
     ! 
-  END SUBROUTINE Read_Ncdf_var1d_Real
+  END SUBROUTINE read_ncdf_var1d_real
   !           
-  !**************************************************************
-  ! end subroutine Read_Ncdf_var1d_real
-  !**************************************************************          
-  !            
-  !           
-  !**************************************************************
-  !   subroutine Read_Ncdf_var2d_real
-  !**************************************************************          
-  !      
-  SUBROUTINE Read_Ncdf_var2d_Real(varname,file,tabvar)
+  !      
+  SUBROUTINE read_ncdf_var2d_real(varname,file,tabvar)
     !      
     IMPLICIT NONE
@@ -250,17 +220,8 @@
     status = nf90_close(ncid)
     !     
-  END SUBROUTINE Read_Ncdf_var2d_Real
+  END SUBROUTINE read_ncdf_var2d_real
   !           
-  !**************************************************************
-  ! end subroutine Read_Ncdf_var2d_real
-  !**************************************************************          
-  !    
-
-  !           
-  !**************************************************************
-  !   subroutine Read_Ncdf_var2d_real_bis
-  !**************************************************************          
-  !      
-  SUBROUTINE Read_Ncdf_var2d_Real_bis(varname,file,tabvar,strt,cnt)
+  !      
+  SUBROUTINE read_ncdf_var2d_real_bis(varname,file,tabvar,strt,cnt)
     !      
     IMPLICIT NONE
@@ -301,18 +262,8 @@
     status = nf90_close(ncid)
     !     
-  END SUBROUTINE Read_Ncdf_var2d_Real_bis
+  END SUBROUTINE read_ncdf_var2d_real_bis
   !           
-  !**************************************************************
-  ! end subroutine Read_Ncdf_var2d_real_bis
-  !**************************************************************          
-  !  
-
-
-  !           
-  !**************************************************************
-  !   subroutine Read_Ncdf_var3d_real
-  !**************************************************************          
-  !      
-  SUBROUTINE Read_Ncdf_var3d_Real(varname,file,tabvar)
+  !      
+  SUBROUTINE read_ncdf_var3d_real(varname,file,tabvar)
     !      
     IMPLICIT NONE
@@ -358,16 +309,8 @@
     status = nf90_close(ncid)
     !     
-  END SUBROUTINE Read_Ncdf_var3d_Real
+  END SUBROUTINE read_ncdf_var3d_real
   !           
-  !**************************************************************
-  ! end subroutine Read_Ncdf_var3d_real
-  !**************************************************************          
-  ! 
-  !           
-  !**************************************************************
-  !   subroutine Read_Ncdf_var4d_real
-  !**************************************************************          
-  !      
-  SUBROUTINE Read_Ncdf_var4d_Real(varname,file,tabvar)
+  !      
+  SUBROUTINE read_ncdf_var4d_real(varname,file,tabvar)
     !      
     IMPLICIT NONE
@@ -410,7 +353,7 @@
     status = nf90_close(ncid)
     !     
-  END SUBROUTINE Read_Ncdf_var4d_Real
-
-  SUBROUTINE Read_Ncdf_var0d_Real(varname,file,tabvar)
+  END SUBROUTINE read_ncdf_var4d_real
+
+  SUBROUTINE read_ncdf_var0d_real(varname,file,tabvar)
     !      
     IMPLICIT NONE
@@ -438,7 +381,7 @@
     status = nf90_close(ncid)
     !     
-  END SUBROUTINE Read_Ncdf_var0d_Real
-
-  SUBROUTINE Read_Ncdf_var0d_Int(varname,file,tabvar)
+  END SUBROUTINE read_ncdf_var0d_real
+
+  SUBROUTINE read_ncdf_var0d_int(varname,file,tabvar)
     !      
     IMPLICIT NONE
@@ -466,16 +409,8 @@
     status = nf90_close(ncid)
     !     
-  END SUBROUTINE Read_Ncdf_var0d_Int
+  END SUBROUTINE read_ncdf_var0d_int
   !           
-  !**************************************************************
-  ! end subroutine Read_Ncdf_var4d_real
-  !**************************************************************          
-  !             
-  !           
-  !**************************************************************
-  !   subroutine Read_Ncdf_var1d_int
-  !**************************************************************          
-  !      
-  SUBROUTINE Read_Ncdf_var1d_Int(varname,file,tabvar)
+  !      
+  SUBROUTINE read_ncdf_var1d_int(varname,file,tabvar)
     !      
     IMPLICIT NONE
@@ -515,16 +450,8 @@
     status = nf90_close(ncid)
     !      
-  END SUBROUTINE Read_Ncdf_var1d_Int
+  END SUBROUTINE read_ncdf_var1d_int
   !           
-  !**************************************************************
-  ! end subroutine Read_Ncdf_var1d_int
-  !**************************************************************          
-  !            
-  !           
-  !**************************************************************
-  !   subroutine Read_Ncdf_var2d_int
-  !**************************************************************          
-  !      
-  SUBROUTINE Read_Ncdf_var2d_Int(varname,file,tabvar)
+  !      
+  SUBROUTINE read_ncdf_var2d_int(varname,file,tabvar)
     !      
     IMPLICIT NONE
@@ -563,17 +490,8 @@
     status = nf90_close(ncid)
     !      
-  END SUBROUTINE Read_Ncdf_var2d_Int
+  END SUBROUTINE read_ncdf_var2d_int
   !           
-  !**************************************************************
-  ! end subroutine Read_Ncdf_var2d_Int
-  !**************************************************************          
-  !    
-
-  !           
-  !**************************************************************
-  !   subroutine Read_Ncdf_var3d_Int
-  !**************************************************************          
-  !      
-  SUBROUTINE Read_Ncdf_var3d_Int(varname,file,tabvar)
+  !      
+  SUBROUTINE read_ncdf_var3d_int(varname,file,tabvar)
     !      
     IMPLICIT NONE
@@ -615,17 +533,8 @@
     status = nf90_close(ncid)
     !      
-  END SUBROUTINE Read_Ncdf_var3d_Int
+  END SUBROUTINE read_ncdf_var3d_int
   !           
-  !**************************************************************
-  ! end subroutine Read_Ncdf_var3d_int
-  !**************************************************************  
-  !
-  !
-  !
-  !**************************************************************
-  !   subroutine Read_Ncdf_var4d_int
-  !**************************************************************          
-  !      
-  SUBROUTINE Read_Ncdf_var4d_Int(varname,file,tabvar)
+  !      
+  SUBROUTINE read_ncdf_var4d_int(varname,file,tabvar)
     !      
     IMPLICIT NONE
@@ -668,13 +577,9 @@
     status = nf90_close(ncid)
     !     
-  END SUBROUTINE Read_Ncdf_var4d_Int
+  END SUBROUTINE read_ncdf_var4d_int
   !           
-  !**************************************************************
-  ! end subroutine Read_Ncdf_var4d_real
-  !**************************************************************         
-  !    
-  !
-  !****************************************************************
-  !   subroutine Write_Ncdf_var					*
+  !
+  !****************************************************************
+  !   subroutine write_ncdf_var					*
   !								*
   ! subroutine to write a variable in a given file 		*
@@ -688,14 +593,10 @@
   !    
   !      
-  !**************************************************************
-  !   subroutine Write_Ncdf_var1d_real
-  !**************************************************************    
-  !      
-  SUBROUTINE Write_Ncdf_var1d_Real(varname,dimname,file,tabvar,typevar)
+  SUBROUTINE write_ncdf_var1d_real(varname,dimname,file,tabvar,typevar)
     !      
     IMPLICIT NONE
     !       
     CHARACTER(*),INTENT(in) :: varname,file,dimname,typevar
-    REAL*8, DIMENSION(:), POINTER :: tabvar
+    REAL*8, DIMENSION(:), INTENT(in) :: tabvar
     !
     ! local variables
@@ -711,5 +612,5 @@
     ENDIF
     !     
-    status = nf90_inq_dimid(ncid,dimname, dimid)
+    status = nf90_inq_dimid(ncid,dimname,dimid)
     status = nf90_inq_varid(ncid,varname,varid)
     status = nf90_redef(ncid)
@@ -725,15 +626,8 @@
     status = nf90_close(ncid)
     ! 
-  END SUBROUTINE Write_Ncdf_var1d_Real
-  !      
-  !**************************************************************
-  !   end subroutine Write_Ncdf_var1d_real
-  !**************************************************************
-  !      
-  !**************************************************************
-  !   subroutine Write_Ncdf_var2d_real
-  !**************************************************************    
-  !      
-  SUBROUTINE Write_Ncdf_var2d_Real_bis(varname,dimname,file,tabvar,nbdim,typevar)
+  END SUBROUTINE write_ncdf_var1d_real
+  !      
+  !      
+  SUBROUTINE write_ncdf_var2d_real_bis(varname,dimname,file,tabvar,nbdim,typevar)
     !      
     IMPLICIT NONE
@@ -743,6 +637,6 @@
     CHARACTER(*), DIMENSION(4) :: dimname
     REAL*8, DIMENSION(:,:) :: tabvar
-    REAL*8, DIMENSION(:,:,:),POINTER :: tabtemp3d
-    REAL*8, DIMENSION(:,:,:,:),POINTER :: tabtemp4d
+    REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: tabtemp3d
+    REAL*8, ALLOCATABLE, DIMENSION(:,:,:,:) :: tabtemp4d
     !
     ! local variables
@@ -796,21 +690,13 @@
     IF(nbdim==3) status = nf90_put_var(ncid,varid,tabtemp3d)
     !      
-    IF(ASSOCIATED( tabtemp3d ) ) DEALLOCATE( tabtemp3d )   
-    IF(ASSOCIATED( tabtemp4d ) ) DEALLOCATE( tabtemp4d )     
-    !
-    status = nf90_close(ncid)
-    ! 
-  END SUBROUTINE Write_Ncdf_var2d_Real_bis
-  !      
-  !**************************************************************
-  !   end subroutine Write_Ncdf_var2d_real
-  !**************************************************************  
-  !   
-  !      
-  !**************************************************************
-  !   subroutine Write_Ncdf_var2d_real
-  !**************************************************************    
-  !      
-  SUBROUTINE Write_Ncdf_var2d_Real(varname,dimname,file,tabvar,typevar)
+    IF(ALLOCATED( tabtemp3d ) ) DEALLOCATE( tabtemp3d )   
+    IF(ALLOCATED( tabtemp4d ) ) DEALLOCATE( tabtemp4d )     
+    !
+    status = nf90_close(ncid)
+    ! 
+  END SUBROUTINE write_ncdf_var2d_real_bis
+  !      
+  !      
+  SUBROUTINE write_ncdf_var2d_real(varname,dimname,file,tabvar,typevar)
     !      
     !      implicit none
@@ -818,5 +704,5 @@
     CHARACTER(*),INTENT(in) :: varname,file,typevar
     CHARACTER(*), DIMENSION(2) :: dimname
-    REAL*8, DIMENSION(:,:), POINTER :: tabvar
+    REAL*8, DIMENSION(:,:), INTENT(in) :: tabvar
     !
     ! local variables
@@ -851,16 +737,8 @@
     status = nf90_close(ncid)
     ! 
-  END SUBROUTINE Write_Ncdf_var2d_Real
-  !      
-  !**************************************************************
-  !   end subroutine Write_Ncdf_var2d_real
-  !**************************************************************  
-  !
-  !      
-  !**************************************************************
-  !   subroutine Write_Ncdf_var3d_real
-  !**************************************************************
-  !      
-  SUBROUTINE Write_Ncdf_var3d_Real(varname,dimname,file,tabvar,typevar)
+  END SUBROUTINE write_ncdf_var2d_real
+  !      
+  !      
+  SUBROUTINE write_ncdf_var3d_real(varname,dimname,file,tabvar,typevar)
     !      
     IMPLICIT NONE
@@ -868,5 +746,5 @@
     CHARACTER(*),INTENT(in) :: varname,file,typevar
     CHARACTER(*),DIMENSION(3),INTENT(in) :: dimname
-    REAL*8, DIMENSION(:,:,:), POINTER :: tabvar
+    REAL*8, DIMENSION(:,:,:), INTENT(in) :: tabvar
     !
     ! local variables
@@ -902,16 +780,8 @@
     status = nf90_close(ncid)
     ! 
-  END SUBROUTINE Write_Ncdf_var3d_Real
-  !      
-  !**************************************************************
-  !   end subroutine Write_Ncdf_var3d_real
-  !**************************************************************
-  ! 
-  !
-  !**************************************************************
-  !   subroutine Write_Ncdf_var4d_real
-  !**************************************************************
-  !      
-  SUBROUTINE Write_Ncdf_var4d_Real(varname,dimname,file,tabvar,typevar)
+  END SUBROUTINE write_ncdf_var3d_real
+  !      
+  !      
+  SUBROUTINE write_ncdf_var4d_real(varname,dimname,file,tabvar,typevar)
     !      
     IMPLICIT NONE
@@ -919,5 +789,5 @@
     CHARACTER(*),INTENT(in) :: varname,file,typevar
     CHARACTER(*),DIMENSION(4),INTENT(in) :: dimname
-    REAL*8, DIMENSION(:,:,:,:), POINTER :: tabvar
+    REAL*8, DIMENSION(:,:,:,:), INTENT(in) :: tabvar
     !
     ! local variables
@@ -959,22 +829,13 @@
     status = nf90_close(ncid)
     ! 
-  END SUBROUTINE Write_Ncdf_var4d_Real
-  !      
-  !**************************************************************
-  !   end subroutine Write_Ncdf_var4d_real
-  !**************************************************************
-  !
-  !
-  !  
-  !**************************************************************
-  !   subroutine Write_Ncdf_var1d_Int
-  !**************************************************************    
-  !      
-  SUBROUTINE Write_Ncdf_var1d_Int(varname,dimname,file,tabvar)
-    !      
-    IMPLICIT NONE
-    !       
-    CHARACTER(*),INTENT(in) :: varname,file,dimname
-    INTEGER, DIMENSION(:), POINTER :: tabvar
+  END SUBROUTINE write_ncdf_var4d_real
+  !      
+  !      
+  SUBROUTINE write_ncdf_var1d_int(varname,dimname,file,tabvar,typevar)
+    !      
+    IMPLICIT NONE
+    !       
+    CHARACTER(*),INTENT(in) :: varname,file,dimname,typevar
+    INTEGER, DIMENSION(:), INTENT(in) :: tabvar
     !
     ! local variables
@@ -1000,22 +861,14 @@
     status = nf90_close(ncid)
     ! 
-  END SUBROUTINE Write_Ncdf_var1d_Int
-  !      
-  !**************************************************************
-  !   end subroutine Write_Ncdf_var1d_Int
-  !**************************************************************
-  !  
-  !      
-  !**************************************************************
-  !   subroutine Write_Ncdf_var2d_Int
-  !**************************************************************    
-  !      
-  SUBROUTINE Write_Ncdf_var2d_Int(varname,dimname,file,tabvar)
-    !      
-    IMPLICIT NONE
-    !       
-    CHARACTER(*),INTENT(in) :: varname,file
-    CHARACTER(*), DIMENSION(2) :: dimname
-    INTEGER, DIMENSION(:,:), POINTER :: tabvar
+  END SUBROUTINE write_ncdf_var1d_int
+  !      
+  !      
+  SUBROUTINE write_ncdf_var2d_int(varname,dimname,file,tabvar,typevar)
+    !      
+    IMPLICIT NONE
+    !       
+    CHARACTER(*), INTENT(in) :: varname,file,typevar
+    CHARACTER(*), DIMENSION(2), INTENT(in) :: dimname
+    INTEGER, DIMENSION(:,:), INTENT(in) :: tabvar
     !
     ! local variables
@@ -1031,10 +884,9 @@
     ENDIF
     !     
-    status = nf90_inq_dimid(ncid,dimname(1), dimid1)
-    status = nf90_inq_dimid(ncid,dimname(2), dimid2)
+    status = nf90_inq_dimid(ncid,dimname(1),dimid1)
+    status = nf90_inq_dimid(ncid,dimname(2),dimid2)
     status = nf90_inq_varid(ncid,varname,varid)
     status = nf90_redef(ncid)
-    status = nf90_def_var(ncid,varname,nf90_int,     &
-         (/dimid1,dimid2/),varid)
+    status = nf90_def_var(ncid,varname,nf90_int,(/dimid1,dimid2/),varid)
     status = nf90_enddef(ncid)
     status = nf90_put_var(ncid,varid,tabvar)     
@@ -1042,22 +894,14 @@
     status = nf90_close(ncid)
     ! 
-  END SUBROUTINE Write_Ncdf_var2d_Int
-  !      
-  !**************************************************************
-  !   end subroutine Write_Ncdf_var2d_Int
-  !**************************************************************  
-  ! 
-  !      
-  !**************************************************************
-  !   subroutine Write_Ncdf_var3d_Int
-  !**************************************************************
-  !      
-  SUBROUTINE Write_Ncdf_var3d_Int(varname,dimname,file,tabvar)
-    !      
-    IMPLICIT NONE
-    !       
-    CHARACTER(*),INTENT(in) :: varname,file
+  END SUBROUTINE write_ncdf_var2d_int
+  !      
+  !      
+  SUBROUTINE write_ncdf_var3d_int(varname,dimname,file,tabvar,typevar)
+    !      
+    IMPLICIT NONE
+    !       
+    CHARACTER(*),INTENT(in) :: varname,file,typevar
     CHARACTER(*),DIMENSION(3),INTENT(in) :: dimname
-    INTEGER, DIMENSION(:,:,:), POINTER :: tabvar
+    INTEGER, DIMENSION(:,:,:), INTENT(in) :: tabvar
     !
     ! local variables
@@ -1085,22 +929,14 @@
     status = nf90_close(ncid)
     ! 
-  END SUBROUTINE Write_Ncdf_var3d_Int
-  !      
-  !**************************************************************
-  !   end subroutine Write_Ncdf_var3d_Int
-  !**************************************************************
-  ! 
-  !
-  !**************************************************************
-  !   subroutine Write_Ncdf_var4d_Int
-  !**************************************************************
-  !      
-  SUBROUTINE Write_Ncdf_var4d_Int(varname,dimname,file,tabvar)
-    !      
-    IMPLICIT NONE
-    !       
-    CHARACTER(*),INTENT(in) :: varname,file
+  END SUBROUTINE write_ncdf_var3d_int
+  !      
+  !      
+  SUBROUTINE write_ncdf_var4d_int(varname,dimname,file,tabvar,typevar)
+    !      
+    IMPLICIT NONE
+    !       
+    CHARACTER(*),INTENT(in) :: varname,file,typevar
     CHARACTER(*),DIMENSION(4),INTENT(in) :: dimname
-    INTEGER, DIMENSION(:,:,:,:), POINTER :: tabvar
+    INTEGER, DIMENSION(:,:,:,:), INTENT(in) :: tabvar
     !
     ! local variables
@@ -1129,13 +965,9 @@
     status = nf90_close(ncid)
     ! 
-  END SUBROUTINE Write_Ncdf_var4d_Int
-  !      
-  !**************************************************************
-  !   end subroutine Write_Ncdf_var4d_Int
-  !************************************************************** 
-  !
-  !
-  !****************************************************************
-  !   subroutine Read_Ncdf_var_t					*
+  END SUBROUTINE write_ncdf_var4d_int
+  !      
+  !
+  !****************************************************************
+  !   subroutine read_ncdf_var_t					*
   !								*
   ! subroutine to read a variable in a given file for time t	*
@@ -1148,10 +980,6 @@
   !****************************************************************
   !
-  !
-  !**************************************************************
-  !   subroutine Read_Ncdf_var3d_real_t
-  !**************************************************************          
-  !      
-  SUBROUTINE Read_Ncdf_var3d_Real_t(varname,file,tabvar,time)
+  !      
+  SUBROUTINE read_ncdf_var3d_real_t(varname,file,tabvar,time)
     !      
     USE agrif_types     
@@ -1201,14 +1029,9 @@
     status = nf90_close(ncid)
     !     
-  END SUBROUTINE Read_Ncdf_var3d_Real_t
+  END SUBROUTINE read_ncdf_var3d_real_t
   !           
-  !**************************************************************
-  ! end subroutine Read_Ncdf_var3d_real_t
-  !**************************************************************          
-  !
-  !
-  !
-  !****************************************************************
-  !   subroutine Write_Ncdf_var_t					*
+  !
+  !****************************************************************
+  !   subroutine write_ncdf_var_t					*
   !								*
   ! subroutine to write a variable in a given file for time t	*
@@ -1223,9 +1046,5 @@
   !
   !      
-  !**************************************************************
-  !   subroutine Write_Ncdf_var3d_real_t
-  !**************************************************************
-  !      
-  SUBROUTINE Write_Ncdf_var3d_Real_t(varname,dimname,file,tabvar,time,typevar)
+  SUBROUTINE write_ncdf_var3d_real_t(varname,dimname,file,tabvar,time,typevar)
     !      
     IMPLICIT NONE
@@ -1234,5 +1053,5 @@
     CHARACTER(*),DIMENSION(3),INTENT(in) :: dimname
     INTEGER :: time
-    REAL*8, DIMENSION(:,:,:), POINTER :: tabvar
+    REAL*8, DIMENSION(:,:,:), INTENT(in) :: tabvar
     !
     ! local variables
@@ -1280,13 +1099,9 @@
     status = nf90_close(ncid)
     ! 
-  END SUBROUTINE Write_Ncdf_var3d_Real_t
-  !      
-  !**************************************************************
-  !   end subroutine Write_Ncdf_var3d_real_t
-  !**************************************************************
-  !
-  !
-  !****************************************************************
-  !   subroutine Read_Ncdf_var_t					*
+  END SUBROUTINE write_ncdf_var3d_real_t
+  !      
+  !
+  !****************************************************************
+  !   subroutine read_ncdf_var_t					*
   !								*
   ! subroutine to read a variable in a given file for time t	*
@@ -1300,10 +1115,6 @@
   !****************************************************************
   !
-  !
-  !**************************************************************
-  !   subroutine Read_Ncdf_var4d_real_nt
-  !**************************************************************          
-  !      
-  SUBROUTINE Read_Ncdf_var4d_Real_nt(varname,file,tabvar,time,level)
+  !      
+  SUBROUTINE read_ncdf_var4d_real_nt(varname,file,tabvar,time,level)
     !      
     USE agrif_types     
@@ -1353,16 +1164,8 @@
     status = nf90_close(ncid)
     !     
-  END SUBROUTINE Read_Ncdf_var4d_Real_nt
+  END SUBROUTINE read_ncdf_var4d_real_nt
   !           
-  !**************************************************************
-  ! end subroutine Read_Ncdf_var4d_real_nt
-  !**************************************************************          
-  !
-  !
-  !**************************************************************
-  !   subroutine Read_Ncdf_var4d_real_t
-  !**************************************************************          
-  !      
-  SUBROUTINE Read_Ncdf_var4d_Real_t(varname,file,tabvar,time)
+  !      
+  SUBROUTINE read_ncdf_var4d_real_t(varname,file,tabvar,time)
     !      
     USE agrif_types     
@@ -1405,13 +1208,8 @@
     status = nf90_close(ncid)
     !     
-  END SUBROUTINE Read_Ncdf_var4d_Real_t
+  END SUBROUTINE read_ncdf_var4d_real_t
   !           
-  !**************************************************************
-  ! end subroutine Read_Ncdf_var4d_real_t
-  !**************************************************************          
-  !
-  !
-  !****************************************************************
-  !   subroutine Write_Ncdf_var_t					*
+  !****************************************************************
+  !   subroutine write_ncdf_var_t					*
   !								*
   ! subroutine to write a variable in a given file for time t	*
@@ -1427,9 +1225,5 @@
   !
   !      
-  !**************************************************************
-  !   subroutine Write_Ncdf_var4d_real_t
-  !**************************************************************
-  !      
-  SUBROUTINE Write_Ncdf_var4d_Real_t(varname,dimname,file,tabvar,time,typevar)
+  SUBROUTINE write_ncdf_var4d_real_t(varname,dimname,file,tabvar,time,typevar)
     !      
     IMPLICIT NONE
@@ -1438,5 +1232,5 @@
     CHARACTER(*),DIMENSION(4),INTENT(in) :: dimname
     INTEGER :: time,level
-    REAL*8, DIMENSION(:,:,:,:), POINTER :: tabvar
+    REAL*8, DIMENSION(:,:,:,:), INTENT(in) :: tabvar
     !
     ! local variables
@@ -1484,16 +1278,8 @@
     status = nf90_close(ncid)
     ! 
-  END SUBROUTINE Write_Ncdf_var4d_Real_t
-  !      
-  !**************************************************************
-  !   end subroutine Write_Ncdf_var4d_real_t
-  !**************************************************************
-  !
-  !      
-  !**************************************************************
-  !   subroutine Write_Ncdf_var4d_real_nt
-  !**************************************************************
-  !      
-  SUBROUTINE Write_Ncdf_var4d_Real_nt(varname,dimname,file,tabvar,time,level,typevar)
+  END SUBROUTINE write_ncdf_var4d_real_t
+  !      
+  !      
+  SUBROUTINE write_ncdf_var4d_real_nt(varname,dimname,file,tabvar,time,level,typevar)
     !      
     IMPLICIT NONE
@@ -1502,5 +1288,5 @@
     CHARACTER(*),DIMENSION(4),INTENT(in) :: dimname
     INTEGER :: time,level
-    REAL*8, DIMENSION(:,:,:,:), POINTER :: tabvar
+    REAL*8, DIMENSION(:,:,:,:), INTENT(in) :: tabvar
     !
     ! local variables
@@ -1548,7 +1334,7 @@
     status = nf90_close(ncid)
     ! 
-  END SUBROUTINE Write_Ncdf_var4d_Real_nt
-
-  SUBROUTINE Write_Ncdf_var0d_Real(varname,file,tabvar,typevar)
+  END SUBROUTINE write_ncdf_var4d_real_nt
+
+  SUBROUTINE write_ncdf_var0d_real(varname,file,tabvar,typevar)
     !      
     IMPLICIT NONE
@@ -1593,12 +1379,41 @@
     status = nf90_close(ncid)
     ! 
-  END SUBROUTINE Write_Ncdf_var0d_Real
-  !      
-  !**************************************************************
-  !   end subroutine Write_Ncdf_var4d_real_nt
-  !**************************************************************
-  !
-  !****************************************************************
-  !   subroutine Read_Ncdf_VarName				*
+  END SUBROUTINE write_ncdf_var0d_real
+
+  SUBROUTINE write_ncdf_var0d_int(varname,file,tabvar,typevar)
+    !      
+    IMPLICIT NONE
+    !       
+    CHARACTER(*),INTENT(in) :: varname,file,typevar
+    INTEGER :: tabvar
+    !
+    ! local variables
+    !
+    INTEGER :: status,ncid
+    INTEGER :: varid
+    !
+    status = nf90_open(file,NF90_WRITE,ncid)
+    IF (status/=nf90_noerr) THEN
+       WRITE(*,*)"unable to open netcdf file : ",file
+       STOP
+    ENDIF
+    !     
+    status = nf90_redef(ncid)
+    status = nf90_def_var(ncid,TRIM(varname),nf90_int,varid)
+    status = nf90_enddef(ncid)
+    status = nf90_put_var(ncid,varid,tabvar)
+    !     
+    IF (status/=nf90_noerr) THEN    
+       WRITE(*,*)"unable to store variable ",varname, &
+	    " in file ",file
+       STOP
+    ENDIF
+    status = nf90_close(ncid)
+    ! 
+  END SUBROUTINE write_ncdf_var0d_int
+
+  !
+  !****************************************************************
+  !   subroutine read_ncdf_VarName				*
   !								*
   ! subroutine to retrieve of all variables 			*
@@ -1610,9 +1425,6 @@
   !****************************************************************
   !
-  !**************************************************************
-  !   subroutine Read_Ncdf_VarName
-  !**************************************************************          
-  !
-  SUBROUTINE Read_Ncdf_VarName(filename,tabvarname)
+  !
+  SUBROUTINE read_ncdf_VarName(filename,tabvarname)
     !      
     CHARACTER(*),INTENT(in) :: filename
@@ -1637,15 +1449,8 @@
     END DO
 
-  END SUBROUTINE Read_Ncdf_Varname
-  !
-  !**************************************************************
-  !   end subroutine Read_Ncdf_VarName
-  !**************************************************************          
-  !
-  !**************************************************************
-  !   subroutine Copy_Ncdf_att
-  !**************************************************************          
-  !
-  SUBROUTINE Copy_Ncdf_att_var(varname,filein,fileout)
+  END SUBROUTINE read_ncdf_Varname
+  !
+  !
+  SUBROUTINE copy_ncdf_att_var(varname,filein,fileout)
     !      
     CHARACTER(*),INTENT(in) :: filein,fileout
@@ -1703,14 +1508,8 @@
     !      print *,'ici2'
     ! 
-  END SUBROUTINE Copy_Ncdf_att_var
-  !
-  !**************************************************************
-  !   end subroutine Copy_Ncdf_att
-  !**************************************************************         
-  !**************************************************************
-  !   subroutine Copy_Ncdf_att
-  !**************************************************************          
-  !
-  SUBROUTINE Copy_Ncdf_att_latlon(varname,filein,fileout,min,max)
+  END SUBROUTINE copy_ncdf_att_var
+  !
+  !
+  SUBROUTINE copy_ncdf_att_latlon(varname,filein,fileout,min,max)
     !      
     CHARACTER(*),INTENT(in) :: filein,fileout
@@ -1760,10 +1559,8 @@
     status = nf90_close(ncid_in)
     status = nf90_close(ncid_out)
-  END SUBROUTINE Copy_Ncdf_att_latlon
+  END SUBROUTINE copy_ncdf_att_latlon
+
   !*************************************************************
   !**************************************************************
-  !**************************************************************
-  !   function Get_NbDims
-  !**************************************************************          
   !
   INTEGER FUNCTION Get_NbDims( varname , filename )
@@ -1784,8 +1581,4 @@
   END FUNCTION Get_NbDims
   !
-  !*************************************************************
-  !**************************************************************
-  !   function Get_NbDims
-  !**************************************************************          
   !
   LOGICAL FUNCTION Dims_Existence( dimname , filename )
@@ -1811,5 +1604,3 @@
   END FUNCTION Dims_Existence
   !
-  !*************************************************************
-  !**************************************************************
 END MODULE io_netcdf
