Index: NEMO/trunk/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg
===================================================================
--- NEMO/trunk/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg	(revision 13215)
+++ NEMO/trunk/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg	(revision 13216)
@@ -158,6 +158,4 @@
 !-----------------------------------------------------------------------
    ln_spc_dyn    = .true.  !  use 0 as special value for dynamics
-   rn_sponge_tra = 1440.   !  coefficient for tracer   sponge layer [m2/s]
-   rn_sponge_dyn = 1440.   !  coefficient for dynamics sponge layer [m2/s]
    ln_chk_bathy  = .true.  !  =T  check the parent bathymetry
 /
Index: NEMO/trunk/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg
===================================================================
--- NEMO/trunk/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg	(revision 13215)
+++ NEMO/trunk/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg	(revision 13216)
@@ -158,6 +158,4 @@
 !-----------------------------------------------------------------------
    ln_spc_dyn    = .true.  !  use 0 as special value for dynamics
-   rn_sponge_tra =  480.   !  coefficient for tracer   sponge layer [m2/s]
-   rn_sponge_dyn =  480.   !  coefficient for dynamics sponge layer [m2/s]
    ln_chk_bathy  = .true.  !  =T  check the parent bathymetry
 /
Index: NEMO/trunk/cfgs/SHARED/namelist_ref
===================================================================
--- NEMO/trunk/cfgs/SHARED/namelist_ref	(revision 13215)
+++ NEMO/trunk/cfgs/SHARED/namelist_ref	(revision 13216)
@@ -643,11 +643,12 @@
 &namagrif      !  AGRIF zoom                                            ("key_agrif")
 !-----------------------------------------------------------------------
-   ln_agrif_2way = .true.  !  activate two way nesting
-   ln_spc_dyn    = .true.  !  use 0 as special value for dynamics
-   rn_sponge_tra = 2880.   !  coefficient for tracer   sponge layer [m2/s]
-   rn_sponge_dyn = 2880.   !  coefficient for dynamics sponge layer [m2/s]
-   rn_trelax_tra = 0.01    !  inverse of relaxation time (in steps) for tracers []
-   rn_trelax_dyn = 0.01    !  inverse of relaxation time (in steps) for dynamics []
-   ln_chk_bathy  = .false. !  =T  check the parent bathymetry
+   ln_agrif_2way   = .true.  !  activate two way nesting
+   ln_init_chfrpar = .false. !  initialize child grids from parent
+   ln_spc_dyn      = .true.  !  use 0 as special value for dynamics
+   rn_sponge_tra   = 0.002   !  coefficient for tracer   sponge layer []
+   rn_sponge_dyn   = 0.002   !  coefficient for dynamics sponge layer []
+   rn_trelax_tra   = 0.01    !  inverse of relaxation time (in steps) for tracers []
+   rn_trelax_dyn   = 0.01    !  inverse of relaxation time (in steps) for dynamics []
+   ln_chk_bathy    = .false. !  =T  check the parent bathymetry
 /
 !-----------------------------------------------------------------------
Index: NEMO/trunk/src/ICE/iceistate.F90
===================================================================
--- NEMO/trunk/src/ICE/iceistate.F90	(revision 13215)
+++ NEMO/trunk/src/ICE/iceistate.F90	(revision 13216)
@@ -32,4 +32,10 @@
    USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
    USE fldread        ! read input fields
+
+# if defined key_agrif
+   USE agrif_oce
+   USE agrif_ice
+   USE agrif_ice_interp 
+# endif   
 
    IMPLICIT NONE
@@ -168,198 +174,242 @@
       ! 2) overwrite some of the fields with namelist parameters or netcdf file
       !------------------------------------------------------------------------
+
+
       IF( ln_iceini ) THEN
          !                             !---------------!
-         IF( ln_iceini_file )THEN      ! Read a file   !
-            !                          !---------------!
-            WHERE( ff_t(:,:) >= 0._wp )   ;   zswitch(:,:) = 1._wp
-            ELSEWHERE                     ;   zswitch(:,:) = 0._wp
+         
+         IF( Agrif_Root() ) THEN
+
+            IF( ln_iceini_file )THEN      ! Read a file   !
+               !                          !---------------!
+               WHERE( ff_t(:,:) >= 0._wp )   ;   zswitch(:,:) = 1._wp
+               ELSEWHERE                     ;   zswitch(:,:) = 0._wp
+               END WHERE
+               !
+               CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step
+               !
+               ! -- mandatory fields -- !
+               zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) * tmask(:,:,1)
+               zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) * tmask(:,:,1)
+               zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) * tmask(:,:,1)
+
+               ! -- optional fields -- !
+               !    if fields do not exist then set them to the values present in the namelist (except for snow and surface temperature)
+               !
+               ! ice salinity
+               IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) &
+                  &     si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
+               !
+               ! temperatures
+               IF    ( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. &
+                  &    TRIM(si(jp_tms)%clrootname) == 'NOT USED' ) THEN
+                  si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
+                  si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
+                  si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
+               ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2
+                  si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 )
+               ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2
+                  si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tsu)%fnow(:,:,1) + 271.15 )
+               ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_su, set T_su = T_s
+                  si(jp_tsu)%fnow(:,:,1) = si(jp_tms)%fnow(:,:,1)
+               ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_su, set T_su = T_i
+                  si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
+               ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_s, set T_s = T_su
+                  si(jp_tms)%fnow(:,:,1) = si(jp_tsu)%fnow(:,:,1)
+               ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_s, set T_s = T_i
+                  si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
+               ENDIF
+               !
+               ! pond concentration
+               IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) &
+                  &     si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc.
+                  &                              * si(jp_ati)%fnow(:,:,1) 
+               !
+               ! pond depth
+               IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) &
+                  &     si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
+               !
+               zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) * tmask(:,:,1)
+               ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) * tmask(:,:,1)
+               zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) * tmask(:,:,1)
+               ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) * tmask(:,:,1)
+               zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) * tmask(:,:,1)
+               zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) * tmask(:,:,1)
+               !
+               ! change the switch for the following
+               WHERE( zat_i_ini(:,:) > 0._wp )   ;   zswitch(:,:) = tmask(:,:,1) 
+               ELSEWHERE                         ;   zswitch(:,:) = 0._wp
+               END WHERE
+
+               !                          !---------------!
+            ELSE                          ! Read namelist !
+               !                          !---------------!
+               ! no ice if (sst - Tfreez) >= thresold
+               WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst )   ;   zswitch(:,:) = 0._wp 
+               ELSEWHERE                                                                    ;   zswitch(:,:) = tmask(:,:,1)
+               END WHERE
+               !
+               ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
+               WHERE( ff_t(:,:) >= 0._wp )
+                  zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:)
+                  zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:)
+                  zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:)
+                  zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:)
+                  ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
+                  zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:)
+                  ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:)
+                  zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 
+                  zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:)
+               ELSEWHERE
+                  zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
+                  zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:)
+                  zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:)
+                  zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:)
+                  ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
+                  zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:)
+                  ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:)
+                  zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
+                  zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:)
+               END WHERE
+               !
+            ENDIF
+
+
+
+            ! make sure ponds = 0 if no ponds scheme
+            IF ( .NOT.ln_pnd ) THEN
+               zapnd_ini(:,:) = 0._wp
+               zhpnd_ini(:,:) = 0._wp
+            ENDIF
+            
+            !-------------!
+            ! fill fields !
+            !-------------!
+            ! select ice covered grid points
+            npti = 0 ; nptidx(:) = 0
+            DO_2D_11_11
+               IF ( zht_i_ini(ji,jj) > 0._wp ) THEN
+                  npti         = npti  + 1
+                  nptidx(npti) = (jj - 1) * jpi + ji
+               ENDIF
+            END_2D
+
+            ! move to 1D arrays: (jpi,jpj) -> (jpi*jpj)
+            CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti)  , zht_i_ini )
+            CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti)  , zht_s_ini )
+            CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti)  , zat_i_ini )
+            CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini )
+            CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini )
+            CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti)  , zt_su_ini )
+            CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti)  , zsm_i_ini )
+            CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti)  , zapnd_ini )
+            CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti)  , zhpnd_ini )
+
+            ! allocate temporary arrays
+            ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), &
+               &      zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) )
+            
+            ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl)
+            CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                                                   &
+               &              zhi_2d          , zhs_2d          , zai_2d         ,                                                   &
+               &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti), a_ip_1d(1:npti), h_ip_1d(1:npti), &
+               &              zti_2d          , zts_2d          , ztsu_2d        , zsi_2d        , zaip_2d        , zhip_2d )
+
+            ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl)
+            DO jl = 1, jpl
+               zti_3d(:,:,jl) = rt0 * tmask(:,:,1)
+               zts_3d(:,:,jl) = rt0 * tmask(:,:,1)
+            END DO
+            CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d   , h_i    )
+            CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d   , h_s    )
+            CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d   , a_i    )
+            CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d   , zti_3d )
+            CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d   , zts_3d )
+            CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d  , t_su   )
+            CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d   , s_i    )
+            CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d  , a_ip   )
+            CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d  , h_ip   )
+
+            ! deallocate temporary arrays
+            DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
+               &        zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d )
+
+            ! calculate extensive and intensive variables
+            CALL ice_var_salprof ! for sz_i
+            DO jl = 1, jpl
+               DO_2D_11_11
+                  v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)
+                  v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)
+                  sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl)
+               END_2D
+            END DO
+            !
+            DO jl = 1, jpl
+               DO_3D_11_11( 1, nlay_s )
+                  t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl)
+                  e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * &
+                     &               rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )
+               END_3D
+            END DO
+            !
+            DO jl = 1, jpl
+               DO_3D_11_11( 1, nlay_i )
+                  t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 
+                  ztmelts          = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K
+                  e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * &
+                     &               rhoi * (  rcpi  * ( ztmelts - t_i(ji,jj,jk,jl) ) + &
+                     &                         rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) &
+                     &                       - rcp   * ( ztmelts - rt0 ) )
+               END_3D
+            END DO
+
+            ! Melt ponds
+            WHERE( a_i > epsi10 )
+               a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:)
+            ELSEWHERE
+               a_ip_frac(:,:,:) = 0._wp
             END WHERE
+            v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
+             
+            ! specific temperatures for coupled runs
+            tn_ice(:,:,:) = t_su(:,:,:)
+            t1_ice(:,:,:) = t_i (:,:,1,:)
             !
-            CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step
-            !
-            ! -- mandatory fields -- !
-            zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) * tmask(:,:,1)
-            zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) * tmask(:,:,1)
-            zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) * tmask(:,:,1)
-
-            ! -- optional fields -- !
-            !    if fields do not exist then set them to the values present in the namelist (except for snow and surface temperature)
-            !
-            ! ice salinity
-            IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) &
-               &     si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
-            !
-            ! temperatures
-            IF    ( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. &
-               &    TRIM(si(jp_tms)%clrootname) == 'NOT USED' ) THEN
-               si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
-               si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
-               si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
-            ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2
-               si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 )
-            ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2
-               si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tsu)%fnow(:,:,1) + 271.15 )
-            ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_su, set T_su = T_s
-               si(jp_tsu)%fnow(:,:,1) = si(jp_tms)%fnow(:,:,1)
-            ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_su, set T_su = T_i
-               si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
-            ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_s, set T_s = T_su
-               si(jp_tms)%fnow(:,:,1) = si(jp_tsu)%fnow(:,:,1)
-            ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_s, set T_s = T_i
-               si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
-            ENDIF
-            !
-            ! pond concentration
-            IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) &
-               &     si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc.
-               &                              * si(jp_ati)%fnow(:,:,1) 
-            !
-            ! pond depth
-            IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) &
-               &     si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
-            !
-            zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) * tmask(:,:,1)
-            ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) * tmask(:,:,1)
-            zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) * tmask(:,:,1)
-            ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) * tmask(:,:,1)
-            zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) * tmask(:,:,1)
-            zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) * tmask(:,:,1)
-            !
-            ! change the switch for the following
-            WHERE( zat_i_ini(:,:) > 0._wp )   ;   zswitch(:,:) = tmask(:,:,1) 
-            ELSEWHERE                         ;   zswitch(:,:) = 0._wp
+         
+#if  defined key_agrif
+         ELSE
+ 
+            Agrif_SpecialValue    = -9999.
+            Agrif_UseSpecialValue = .TRUE.
+            CALL Agrif_init_variable(tra_iceini_id,procname=interp_tra_ice)
+            use_sign_north = .TRUE.
+            sign_north = -1.
+            CALL Agrif_init_variable(u_iceini_id  ,procname=interp_u_ice)
+            CALL Agrif_init_variable(v_iceini_id  ,procname=interp_v_ice)
+            Agrif_SpecialValue    = 0._wp
+            use_sign_north = .FALSE.
+            Agrif_UseSpecialValue = .FALSE.
+        ! lbc ???? 
+   ! Here we know : a_i, v_i, v_s, sv_i, oa_i, a_ip, v_ip, t_su, e_s, e_i
+            CALL ice_var_glo2eqv
+            CALL ice_var_zapsmall
+            CALL ice_var_agg(2)
+
+            ! Melt ponds
+            WHERE( a_i > epsi10 )
+               a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:)
+            ELSEWHERE
+               a_ip_frac(:,:,:) = 0._wp
             END WHERE
-            !                          !---------------!
-         ELSE                          ! Read namelist !
-            !                          !---------------!
-            ! no ice if (sst - Tfreez) >= thresold
-            WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst )   ;   zswitch(:,:) = 0._wp 
-            ELSEWHERE                                                                    ;   zswitch(:,:) = tmask(:,:,1)
-            END WHERE
-            !
-            ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
-            WHERE( ff_t(:,:) >= 0._wp )
-               zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:)
-               zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:)
-               zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:)
-               zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:)
-               ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
-               zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:)
-               ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:)
-               zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 
-               zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:)
+            WHERE( a_ip > 0._wp )       ! ???????    
+               h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:)
             ELSEWHERE
-               zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
-               zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:)
-               zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:)
-               zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:)
-               ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
-               zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:)
-               ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:)
-               zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
-               zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:)
-            END WHERE
-            !
-         ENDIF
-
-         ! make sure ponds = 0 if no ponds scheme
-         IF ( .NOT.ln_pnd ) THEN
-            zapnd_ini(:,:) = 0._wp
-            zhpnd_ini(:,:) = 0._wp
-         ENDIF
-         
-         !-------------!
-         ! fill fields !
-         !-------------!
-         ! select ice covered grid points
-         npti = 0 ; nptidx(:) = 0
-         DO_2D_11_11
-            IF ( zht_i_ini(ji,jj) > 0._wp ) THEN
-               npti         = npti  + 1
-               nptidx(npti) = (jj - 1) * jpi + ji
-            ENDIF
-         END_2D
-
-         ! move to 1D arrays: (jpi,jpj) -> (jpi*jpj)
-         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti)  , zht_i_ini )
-         CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti)  , zht_s_ini )
-         CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti)  , zat_i_ini )
-         CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini )
-         CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini )
-         CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti)  , zt_su_ini )
-         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti)  , zsm_i_ini )
-         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti)  , zapnd_ini )
-         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti)  , zhpnd_ini )
-
-         ! allocate temporary arrays
-         ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), &
-            &      zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) )
-         
-         ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl)
-         CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                                                   &
-            &              zhi_2d          , zhs_2d          , zai_2d         ,                                                   &
-            &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti), a_ip_1d(1:npti), h_ip_1d(1:npti), &
-            &              zti_2d          , zts_2d          , ztsu_2d        , zsi_2d        , zaip_2d        , zhip_2d )
-
-         ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl)
-         DO jl = 1, jpl
-            zti_3d(:,:,jl) = rt0 * tmask(:,:,1)
-            zts_3d(:,:,jl) = rt0 * tmask(:,:,1)
-         END DO
-         CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d   , h_i    )
-         CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d   , h_s    )
-         CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d   , a_i    )
-         CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d   , zti_3d )
-         CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d   , zts_3d )
-         CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d  , t_su   )
-         CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d   , s_i    )
-         CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d  , a_ip   )
-         CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d  , h_ip   )
-
-         ! deallocate temporary arrays
-         DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
-            &        zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d )
-
-         ! calculate extensive and intensive variables
-         CALL ice_var_salprof ! for sz_i
-         DO jl = 1, jpl
-            DO_2D_11_11
-               v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)
-               v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)
-               sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl)
-            END_2D
-         END DO
-         !
-         DO jl = 1, jpl
-            DO_3D_11_11( 1, nlay_s )
-               t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl)
-               e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * &
-                  &               rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )
-            END_3D
-         END DO
-         !
-         DO jl = 1, jpl
-            DO_3D_11_11( 1, nlay_i )
-               t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 
-               ztmelts          = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K
-               e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * &
-                  &               rhoi * (  rcpi  * ( ztmelts - t_i(ji,jj,jk,jl) ) + &
-                  &                         rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) &
-                  &                       - rcp   * ( ztmelts - rt0 ) )
-            END_3D
-         END DO
-
-         ! Melt ponds
-         WHERE( a_i > epsi10 )
-            a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:)
-         ELSEWHERE
-            a_ip_frac(:,:,:) = 0._wp
-         END WHERE
-         v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
-          
-         ! specific temperatures for coupled runs
-         tn_ice(:,:,:) = t_su(:,:,:)
-         t1_ice(:,:,:) = t_i (:,:,1,:)
-         !
+               h_ip(:,:,:) = 0._wp
+            END WHERE   
+
+            tn_ice(:,:,:) = t_su(:,:,:)
+            t1_ice(:,:,:) = t_i (:,:,1,:)
+#endif
+          ENDIF ! Agrif_Root
       ENDIF ! ln_iceini
       !
Index: NEMO/trunk/src/ICE/icestp.F90
===================================================================
--- NEMO/trunk/src/ICE/icestp.F90	(revision 13215)
+++ NEMO/trunk/src/ICE/icestp.F90	(revision 13216)
@@ -240,4 +240,8 @@
       CALL par_init                ! set some ice run parameters
       !
+#if defined key_agrif
+      CALL Agrif_Declare_Var_ice  !  "      "   "   "      "  Sea ice
+#endif
+      !
       !                                ! Allocate the ice arrays (sbc_ice already allocated in sbc_init)
       ierr =        ice_alloc        ()      ! ice variables
Index: NEMO/trunk/src/NST/agrif_ice.F90
===================================================================
--- NEMO/trunk/src/NST/agrif_ice.F90	(revision 13215)
+++ NEMO/trunk/src/NST/agrif_ice.F90	(revision 13216)
@@ -16,4 +16,5 @@
 
    INTEGER, PUBLIC ::  u_ice_id, v_ice_id, tra_ice_id
+   INTEGER, PUBLIC ::  u_iceini_id, v_iceini_id, tra_iceini_id
    INTEGER, PUBLIC ::  nbstep_ice = 0    ! child time position in sea-ice model
 
Index: NEMO/trunk/src/NST/agrif_ice_interp.F90
===================================================================
--- NEMO/trunk/src/NST/agrif_ice_interp.F90	(revision 13215)
+++ NEMO/trunk/src/NST/agrif_ice_interp.F90	(revision 13216)
@@ -14,7 +14,7 @@
    !!----------------------------------------------------------------------
    !!  agrif_interp_ice    : interpolation of ice at "after" sea-ice time step
-   !!  agrif_interp_u_ice   : atomic routine to interpolate u_ice 
-   !!  agrif_interp_v_ice   : atomic routine to interpolate v_ice 
-   !!  agrif_interp_tra_ice : atomic routine to interpolate ice properties 
+   !!  interp_u_ice   : atomic routine to interpolate u_ice 
+   !!  interp_v_ice   : atomic routine to interpolate v_ice 
+   !!  interp_tra_ice : atomic routine to interpolate ice properties 
    !!----------------------------------------------------------------------
    USE par_oce
@@ -23,4 +23,5 @@
    USE ice
    USE agrif_ice
+   USE agrif_oce
    USE phycst , ONLY: rt0
    
@@ -29,4 +30,5 @@
 
    PUBLIC   agrif_interp_ice   ! called by agrif_user.F90
+   PUBLIC   interp_tra_ice, interp_u_ice, interp_v_ice  ! called by iceistate.F90
 
    !!----------------------------------------------------------------------
@@ -68,4 +70,9 @@
       Agrif_SpecialValue    = -9999.
       Agrif_UseSpecialValue = .TRUE.
+
+      use_sign_north = .TRUE.
+      sign_north = -1.
+      if (cd_type == 'T') use_sign_north = .FALSE.
+
       SELECT CASE( cd_type )
       CASE('U')   ;   CALL Agrif_Bc_variable( u_ice_id  , procname=interp_u_ice  , calledweight=zbeta )
@@ -75,4 +82,6 @@
       Agrif_SpecialValue    = 0._wp
       Agrif_UseSpecialValue = .FALSE.
+      
+      use_sign_north = .FALSE.
       !
    END SUBROUTINE agrif_interp_ice
@@ -156,5 +165,5 @@
       ! and it is ok since we conserve tracers (same as in the ocean).
       ALLOCATE( ztab(SIZE(a_i,1),SIZE(a_i,2),SIZE(ptab,3)) )
-     
+
       IF( before ) THEN  ! parent grid
          jm = 1
Index: NEMO/trunk/src/NST/agrif_ice_update.F90
===================================================================
--- NEMO/trunk/src/NST/agrif_ice_update.F90	(revision 13215)
+++ NEMO/trunk/src/NST/agrif_ice_update.F90	(revision 13216)
@@ -66,4 +66,7 @@
       CALL Agrif_Update_Variable( tra_ice_id , locupdate=(/1,0/), procname = update_tra_ice  )
 #endif
+      use_sign_north = .TRUE.
+      sign_north = -1.
+
 # if ! defined DECAL_FEEDBACK
       CALL Agrif_Update_Variable( u_ice_id   , procname = update_u_ice    )
@@ -73,4 +76,5 @@
       CALL Agrif_Update_Variable( v_ice_id   , locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname=update_v_ice)
 #endif
+      use_sign_north = .FALSE.
 !      CALL Agrif_Update_Variable( tra_ice_id , locupdate=(/0,2/), procname = update_tra_ice  )
 !      CALL Agrif_Update_Variable( u_ice_id   , locupdate=(/0,1/), procname = update_u_ice    )
Index: NEMO/trunk/src/NST/agrif_oce.F90
===================================================================
--- NEMO/trunk/src/NST/agrif_oce.F90	(revision 13215)
+++ NEMO/trunk/src/NST/agrif_oce.F90	(revision 13216)
@@ -19,4 +19,5 @@
   
    !                                              !!* Namelist namagrif: AGRIF parameters
+   LOGICAL , PUBLIC ::   ln_init_chfrpar = .FALSE. !: set child grids initial state from parent
    LOGICAL , PUBLIC ::   ln_agrif_2way = .TRUE.    !: activate two way nesting 
    LOGICAL , PUBLIC ::   ln_spc_dyn    = .FALSE.   !: use zeros (.false.) or not (.true.) in
@@ -29,4 +30,5 @@
    !
    INTEGER , PUBLIC, PARAMETER ::   nn_sponge_len = 2  !: Sponge width (in number of parent grid points)
+
    LOGICAL , PUBLIC :: spongedoneT = .FALSE.       !: tracer   sponge layer indicator
    LOGICAL , PUBLIC :: spongedoneU = .FALSE.       !: dynamics sponge layer indicator
@@ -49,8 +51,6 @@
    INTEGER , PUBLIC,              SAVE                 ::  Kbb_a, Kmm_a, Krhs_a   !: AGRIF module-specific copies of time-level indices
 
-# if defined key_vertical
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht0_parent, hu0_parent, hv0_parent
    INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt_parent, mbku_parent, mbkv_parent
-# endif
 
    INTEGER, PUBLIC :: tsn_id                                                  ! AGRIF profile for tracers interpolation and update
@@ -58,4 +58,5 @@
    INTEGER, PUBLIC :: un_update_id, vn_update_id                              ! AGRIF profiles for udpates
    INTEGER, PUBLIC :: tsn_sponge_id, un_sponge_id, vn_sponge_id               ! AGRIF profiles for sponge layers
+   INTEGER, PUBLIC :: tsini_id, uini_id, vini_id, sshini_id                   ! AGRIF profile for initialization
 # if defined key_top
    INTEGER, PUBLIC :: trn_id, trn_sponge_id
@@ -68,4 +69,16 @@
    INTEGER, PUBLIC :: mbkt_id, ht0_id
    INTEGER, PUBLIC :: kindic_agr
+
+   ! North fold
+!$AGRIF_DO_NOT_TREAT
+   LOGICAL, PUBLIC :: use_sign_north
+   REAL, PUBLIC :: sign_north
+   LOGICAL, PUBLIC :: l_ini_child = .FALSE.
+# if defined key_vertical
+   LOGICAL, PUBLIC :: l_vremap    = .TRUE.
+# else
+   LOGICAL, PUBLIC :: l_vremap    = .FALSE.
+# endif
+!$AGRIF_END_DO_NOT_TREAT
    
    !!----------------------------------------------------------------------
@@ -91,9 +104,7 @@
          &      tabspongedone_trn(jpi,jpj),           &
 # endif   
-# if defined key_vertical
          &      ht0_parent(jpi,jpj), mbkt_parent(jpi,jpj),  &
          &      hu0_parent(jpi,jpj), mbku_parent(jpi,jpj),  &
          &      hv0_parent(jpi,jpj), mbkv_parent(jpi,jpj),  &
-# endif      
          &      tabspongedone_u  (jpi,jpj),           &
          &      tabspongedone_v  (jpi,jpj), STAT = ierr(1) )
Index: NEMO/trunk/src/NST/agrif_oce_interp.F90
===================================================================
--- NEMO/trunk/src/NST/agrif_oce_interp.F90	(revision 13215)
+++ NEMO/trunk/src/NST/agrif_oce_interp.F90	(revision 13216)
@@ -34,4 +34,5 @@
    USE lib_mpp
    USE vremap
+   USE lbclnk
  
    IMPLICIT NONE
@@ -44,7 +45,7 @@
    PUBLIC   interpunb, interpvnb , interpub2b, interpvb2b
    PUBLIC   interpe3t
-#if defined key_vertical
    PUBLIC   interpht0, interpmbkt
-# endif
+   PUBLIC   agrif_initts, agrif_initssh
+
    INTEGER ::   bdy_tinterp = 0
 
@@ -89,251 +90,262 @@
       Agrif_UseSpecialValue = ln_spc_dyn
       !
+      use_sign_north = .TRUE.
+      sign_north = -1.
       CALL Agrif_Bc_variable( un_interp_id, procname=interpun )
       CALL Agrif_Bc_variable( vn_interp_id, procname=interpvn )
+      use_sign_north = .FALSE.
       !
       Agrif_UseSpecialValue = .FALSE.
       !
       ! --- West --- !
-      ibdy1 = 2
-      ibdy2 = 1+nbghostcells 
-      !
-      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
+      IF( lk_west ) THEN
+         ibdy1 = 2
+         ibdy2 = 1+nbghostcells 
+         !
+         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
+            DO ji = mi0(ibdy1), mi1(ibdy2)
+               uu_b(ji,:,Krhs_a) = 0._wp
+
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpj
+                     uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
+                  END DO
+               END DO
+
+               DO jj = 1, jpj
+                  uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) * r1_hu(ji,jj,Krhs_a)
+               END DO
+            END DO
+         ENDIF
+         !
          DO ji = mi0(ibdy1), mi1(ibdy2)
-            uu_b(ji,:,Krhs_a) = 0._wp
-
+            zub(ji,:) = 0._wp    ! Correct transport
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
-               END DO
-            END DO
-
-            DO jj = 1, jpj
-               uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) * r1_hu(ji,jj,Krhs_a)
-            END DO
-         END DO
-      ENDIF
-      !
-      DO ji = mi0(ibdy1), mi1(ibdy2)
-         zub(ji,:) = 0._wp    ! Correct transport
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               zub(ji,jj) = zub(ji,jj) & 
-                  & + e3u(ji,jj,jk,Krhs_a)  * uu(ji,jj,jk,Krhs_a)*umask(ji,jj,jk)
-            END DO
-         END DO
-         DO jj=1,jpj
-            zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
-         END DO
-            
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) + uu_b(ji,jj,Krhs_a)-zub(ji,jj)) * umask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-            
-      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
-         DO ji = mi0(ibdy1), mi1(ibdy2)
-            zvb(ji,:) = 0._wp
+                  zub(ji,jj) = zub(ji,jj) & 
+                     & + e3u(ji,jj,jk,Krhs_a)  * uu(ji,jj,jk,Krhs_a)*umask(ji,jj,jk)
+               END DO
+            END DO
+            DO jj=1,jpj
+               zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
+            END DO
+               
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
-               END DO
-            END DO
-            DO jj = 1, jpj
-               zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a)
-            END DO
+                  uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) + uu_b(ji,jj,Krhs_a)-zub(ji,jj)) * umask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+               
+         IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
+            DO ji = mi0(ibdy1), mi1(ibdy2)
+               zvb(ji,:) = 0._wp
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpj
+                     zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
+                  END DO
+               END DO
+               DO jj = 1, jpj
+                  zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a)
+               END DO
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpj
+                     vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) + vv_b(ji,jj,Krhs_a)-zvb(ji,jj))*vmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         ENDIF
+      ENDIF
+
+      ! --- East --- !
+      IF( lk_east) THEN
+         ibdy1 = jpiglo-1-nbghostcells
+         ibdy2 = jpiglo-2 
+         !
+         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
+            DO ji = mi0(ibdy1), mi1(ibdy2)
+               uu_b(ji,:,Krhs_a) = 0._wp
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpj
+                     uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) & 
+                         & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
+                  END DO
+               END DO
+               DO jj = 1, jpj
+                  uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) * r1_hu(ji,jj,Krhs_a)
+               END DO
+            END DO
+         ENDIF
+         !
+         DO ji = mi0(ibdy1), mi1(ibdy2)
+            zub(ji,:) = 0._wp    ! Correct transport
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) + vv_b(ji,jj,Krhs_a)-zvb(ji,jj))*vmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-      ENDIF
-
-      ! --- East --- !
-      ibdy1 = jpiglo-1-nbghostcells
-      ibdy2 = jpiglo-2 
-      !
-      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
-         DO ji = mi0(ibdy1), mi1(ibdy2)
-            uu_b(ji,:,Krhs_a) = 0._wp
+                  zub(ji,jj) = zub(ji,jj) & 
+                     & + e3u(ji,jj,jk,Krhs_a)  * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
+               END DO
+            END DO
+            DO jj=1,jpj
+               zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
+            END DO
+               
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) & 
-                      & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
-               END DO
-            END DO
-            DO jj = 1, jpj
-               uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) * r1_hu(ji,jj,Krhs_a)
-            END DO
-         END DO
-      ENDIF
-      !
-      DO ji = mi0(ibdy1), mi1(ibdy2)
-         zub(ji,:) = 0._wp    ! Correct transport
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               zub(ji,jj) = zub(ji,jj) & 
-                  & + e3u(ji,jj,jk,Krhs_a)  * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
-            END DO
-         END DO
-         DO jj=1,jpj
-            zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
-         END DO
-            
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) & 
-                 & + uu_b(ji,jj,Krhs_a)-zub(ji,jj))*umask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-            
-      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
-         ibdy1 = jpiglo-nbghostcells
-         ibdy2 = jpiglo-1 
-         DO ji = mi0(ibdy1), mi1(ibdy2)
-            zvb(ji,:) = 0._wp
-            DO jk = 1, jpkm1
+                  uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) & 
+                    & + uu_b(ji,jj,Krhs_a)-zub(ji,jj))*umask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+               
+         IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
+            ibdy1 = jpiglo-nbghostcells
+            ibdy2 = jpiglo-1 
+            DO ji = mi0(ibdy1), mi1(ibdy2)
+               zvb(ji,:) = 0._wp
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpj
+                     zvb(ji,jj) = zvb(ji,jj) &
+                        & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
+                  END DO
+               END DO
                DO jj = 1, jpj
-                  zvb(ji,jj) = zvb(ji,jj) &
+                  zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a)
+               END DO
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpj
+                     vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) & 
+                         & + vv_b(ji,jj,Krhs_a)-zvb(ji,jj)) * vmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         ENDIF
+      ENDIF
+
+      ! --- South --- !
+      IF( lk_south ) THEN
+         jbdy1 = 2
+         jbdy2 = 1+nbghostcells 
+         !
+         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
+            DO jj = mj0(jbdy1), mj1(jbdy2)
+               vv_b(:,jj,Krhs_a) = 0._wp
+               DO jk = 1, jpkm1
+                  DO ji = 1, jpi
+                     vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) & 
+                         & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
+                  END DO
+               END DO
+               DO ji=1,jpi
+                  vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) * r1_hv(ji,jj,Krhs_a)     
+               END DO
+            END DO
+         ENDIF
+         !
+         DO jj = mj0(jbdy1), mj1(jbdy2)
+            zvb(:,jj) = 0._wp    ! Correct transport
+            DO jk=1,jpkm1
+               DO ji=1,jpi
+                  zvb(ji,jj) = zvb(ji,jj) & 
                      & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
                END DO
             END DO
-            DO jj = 1, jpj
+            DO ji = 1, jpi
                zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a)
             END DO
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) & 
-                      & + vv_b(ji,jj,Krhs_a)-zvb(ji,jj)) * vmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-      ENDIF
-
-      ! --- South --- !
-      jbdy1 = 2
-      jbdy2 = 1+nbghostcells 
-      !
-      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
-         DO jj = mj0(jbdy1), mj1(jbdy2)
-            vv_b(:,jj,Krhs_a) = 0._wp
+
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) & 
-                      & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
-               END DO
-            END DO
-            DO ji=1,jpi
-               vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) * r1_hv(ji,jj,Krhs_a)     
-            END DO
-         END DO
-      ENDIF
-      !
-      DO jj = mj0(jbdy1), mj1(jbdy2)
-         zvb(:,jj) = 0._wp    ! Correct transport
-         DO jk=1,jpkm1
-            DO ji=1,jpi
-               zvb(ji,jj) = zvb(ji,jj) & 
-                  & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
-            END DO
-         END DO
-         DO ji = 1, jpi
-            zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a)
-         END DO
-
-         DO jk = 1, jpkm1
+                  vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) & 
+                    & + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) ) * vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+               
+         IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
+            DO jj = mj0(jbdy1), mj1(jbdy2)
+               zub(:,jj) = 0._wp
+               DO jk = 1, jpkm1
+                  DO ji = 1, jpi
+                     zub(ji,jj) = zub(ji,jj) & 
+                        & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
+                  END DO
+               END DO
+               DO ji = 1, jpi
+                  zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
+               END DO
+                  
+               DO jk = 1, jpkm1
+                  DO ji = 1, jpi
+                     uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) & 
+                       & + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         ENDIF
+      ENDIF
+
+      ! --- North --- !
+      IF( lk_north ) THEN
+         jbdy1 = jpjglo-1-nbghostcells
+         jbdy2 = jpjglo-2 
+         !
+         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
+            DO jj = mj0(jbdy1), mj1(jbdy2)
+               vv_b(:,jj,Krhs_a) = 0._wp
+               DO jk = 1, jpkm1
+                  DO ji = 1, jpi
+                     vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) & 
+                         & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
+                  END DO
+               END DO
+               DO ji=1,jpi
+                  vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) * r1_hv(ji,jj,Krhs_a)
+               END DO
+            END DO
+         ENDIF
+         !
+         DO jj = mj0(jbdy1), mj1(jbdy2)
+            zvb(:,jj) = 0._wp    ! Correct transport
+            DO jk=1,jpkm1
+               DO ji=1,jpi
+                  zvb(ji,jj) = zvb(ji,jj) & 
+                     & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
+               END DO
+            END DO
             DO ji = 1, jpi
-               vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) & 
-                 & + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) ) * vmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-            
-      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
-         DO jj = mj0(jbdy1), mj1(jbdy2)
-            zub(:,jj) = 0._wp
+               zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a)
+            END DO
+
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  zub(ji,jj) = zub(ji,jj) & 
-                     & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
-               END DO
-            END DO
-            DO ji = 1, jpi
-               zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
-            END DO
+                  vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) & 
+                    & + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) ) * vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
                
-            DO jk = 1, jpkm1
+         IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
+            jbdy1 = jpjglo-nbghostcells
+            jbdy2 = jpjglo-1
+            DO jj = mj0(jbdy1), mj1(jbdy2)
+               zub(:,jj) = 0._wp
+               DO jk = 1, jpkm1
+                  DO ji = 1, jpi
+                     zub(ji,jj) = zub(ji,jj) & 
+                        & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
+                  END DO
+               END DO
                DO ji = 1, jpi
-                  uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) & 
-                    & + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-      ENDIF
-
-      ! --- North --- !
-      jbdy1 = jpjglo-1-nbghostcells
-      jbdy2 = jpjglo-2 
-      !
-      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
-         DO jj = mj0(jbdy1), mj1(jbdy2)
-            vv_b(:,jj,Krhs_a) = 0._wp
-            DO jk = 1, jpkm1
-               DO ji = 1, jpi
-                  vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) & 
-                      & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
-               END DO
-            END DO
-            DO ji=1,jpi
-               vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) * r1_hv(ji,jj,Krhs_a)
-            END DO
-         END DO
-      ENDIF
-      !
-      DO jj = mj0(jbdy1), mj1(jbdy2)
-         zvb(:,jj) = 0._wp    ! Correct transport
-         DO jk=1,jpkm1
-            DO ji=1,jpi
-               zvb(ji,jj) = zvb(ji,jj) & 
-                  & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
-            END DO
-         END DO
-         DO ji = 1, jpi
-            zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a)
-         END DO
-
-         DO jk = 1, jpkm1
-            DO ji = 1, jpi
-               vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) & 
-                 & + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) ) * vmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-            
-      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
-         jbdy1 = jpjglo-nbghostcells
-         jbdy2 = jpjglo-1
-         DO jj = mj0(jbdy1), mj1(jbdy2)
-            zub(:,jj) = 0._wp
-            DO jk = 1, jpkm1
-               DO ji = 1, jpi
-                  zub(ji,jj) = zub(ji,jj) & 
-                     & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
-               END DO
-            END DO
-            DO ji = 1, jpi
-               zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
-            END DO
-               
-            DO jk = 1, jpkm1
-               DO ji = 1, jpi
-                  uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) & 
-                    & + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+                  zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
+               END DO
+                  
+               DO jk = 1, jpkm1
+                  DO ji = 1, jpi
+                     uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) & 
+                       & + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         ENDIF
       ENDIF
       !
@@ -354,54 +366,64 @@
       !
       !--- West ---!
-      istart = 2
-      iend   = nbghostcells+1
-      DO ji = mi0(istart), mi1(iend)
-         DO jj=1,jpj
-            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
-            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
-         END DO
-      END DO
+      IF( lk_west ) THEN
+         istart = 2
+         iend   = nbghostcells+1
+         DO ji = mi0(istart), mi1(iend)
+            DO jj=1,jpj
+               va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
+               ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
+            END DO
+         END DO
+      ENDIF
       !
       !--- East ---!
-      istart = jpiglo-nbghostcells
-      iend   = jpiglo-1
-      DO ji = mi0(istart), mi1(iend)
-         DO jj=1,jpj
-            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
-         END DO
-      END DO
-      istart = jpiglo-nbghostcells-1
-      iend   = jpiglo-2
-      DO ji = mi0(istart), mi1(iend)
-         DO jj=1,jpj
-            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
-         END DO
-      END DO
+      IF( lk_east ) THEN
+         istart = jpiglo-nbghostcells
+         iend   = jpiglo-1
+         DO ji = mi0(istart), mi1(iend)
+
+            DO jj=1,jpj
+               va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
+            END DO
+         END DO
+         istart = jpiglo-nbghostcells-1
+         iend   = jpiglo-2
+         DO ji = mi0(istart), mi1(iend)
+            DO jj=1,jpj
+               ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
+            END DO
+         END DO
+      ENDIF 
       !
       !--- South ---!
-      jstart = 2
-      jend   = nbghostcells+1
-      DO jj = mj0(jstart), mj1(jend)
-         DO ji=1,jpi
-            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
-            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
-         END DO
-      END DO
+      IF( lk_south ) THEN
+         jstart = 2
+         jend   = nbghostcells+1
+         DO jj = mj0(jstart), mj1(jend)
+
+            DO ji=1,jpi
+               ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
+               va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
+            END DO
+         END DO
+      ENDIF       
       !
       !--- North ---!
-      jstart = jpjglo-nbghostcells
-      jend   = jpjglo-1
-      DO jj = mj0(jstart), mj1(jend)
-         DO ji=1,jpi
-            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
-         END DO
-      END DO
-      jstart = jpjglo-nbghostcells-1
-      jend   = jpjglo-2
-      DO jj = mj0(jstart), mj1(jend)
-         DO ji=1,jpi
-            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
-         END DO
-      END DO
+      IF( lk_north ) THEN
+         jstart = jpjglo-nbghostcells
+         jend   = jpjglo-1
+         DO jj = mj0(jstart), mj1(jend)
+            DO ji=1,jpi
+               ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
+            END DO
+         END DO
+         jstart = jpjglo-nbghostcells-1
+         jend   = jpjglo-2
+         DO jj = mj0(jstart), mj1(jend)
+            DO ji=1,jpi
+               va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
+            END DO
+         END DO
+      ENDIF 
       !
    END SUBROUTINE Agrif_dyn_ts
@@ -421,54 +443,62 @@
       !
       !--- West ---!
-      istart = 2
-      iend   = nbghostcells+1
-      DO ji = mi0(istart), mi1(iend)
-         DO jj=1,jpj
-            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
-            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
-         END DO
-      END DO
+      IF( lk_west ) THEN
+         istart = 2
+         iend   = nbghostcells+1
+         DO ji = mi0(istart), mi1(iend)
+            DO jj=1,jpj
+               zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
+               zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
+            END DO
+         END DO
+      ENDIF
       !
       !--- East ---!
-      istart = jpiglo-nbghostcells
-      iend   = jpiglo-1
-      DO ji = mi0(istart), mi1(iend)
-         DO jj=1,jpj
-            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
-         END DO
-      END DO
-      istart = jpiglo-nbghostcells-1
-      iend   = jpiglo-2
-      DO ji = mi0(istart), mi1(iend)
-         DO jj=1,jpj
-            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
-         END DO
-      END DO
+      IF( lk_east ) THEN
+         istart = jpiglo-nbghostcells
+         iend   = jpiglo-1
+         DO ji = mi0(istart), mi1(iend)
+            DO jj=1,jpj
+               zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
+            END DO
+         END DO
+         istart = jpiglo-nbghostcells-1
+         iend   = jpiglo-2
+         DO ji = mi0(istart), mi1(iend)
+            DO jj=1,jpj
+               zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
+            END DO
+         END DO
+      ENDIF
       !
       !--- South ---!
-      jstart = 2
-      jend   = nbghostcells+1
-      DO jj = mj0(jstart), mj1(jend)
-         DO ji=1,jpi
-            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
-            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
-         END DO
-      END DO
+      IF( lk_south ) THEN
+         jstart = 2
+         jend   = nbghostcells+1
+         DO jj = mj0(jstart), mj1(jend)
+            DO ji=1,jpi
+               zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
+               zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
+            END DO
+         END DO
+      ENDIF
       !
       !--- North ---!
-      jstart = jpjglo-nbghostcells
-      jend   = jpjglo-1
-      DO jj = mj0(jstart), mj1(jend)
-         DO ji=1,jpi
-            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
-         END DO
-      END DO
-      jstart = jpjglo-nbghostcells-1
-      jend   = jpjglo-2
-      DO jj = mj0(jstart), mj1(jend)
-         DO ji=1,jpi
-            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
-         END DO
-      END DO
+      IF( lk_north ) THEN
+         jstart = jpjglo-nbghostcells
+         jend   = jpjglo-1
+         DO jj = mj0(jstart), mj1(jend)
+            DO ji=1,jpi
+               zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
+            END DO
+         END DO
+         jstart = jpjglo-nbghostcells-1
+         jend   = jpjglo-2
+         DO jj = mj0(jstart), mj1(jend)
+            DO ji=1,jpi
+               zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
+            END DO
+         END DO
+      ENDIF
       !
    END SUBROUTINE Agrif_dyn_ts_flux
@@ -494,4 +524,8 @@
       Agrif_SpecialValue = 0._wp
       Agrif_UseSpecialValue = ln_spc_dyn
+
+      use_sign_north = .TRUE.
+      sign_north = -1.
+
       !
       ! Set bdy time interpolation stage to 0 (latter incremented locally do deal with corners)
@@ -518,4 +552,5 @@
       ENDIF
       Agrif_UseSpecialValue = .FALSE.
+      use_sign_north = .FALSE.
       ! 
    END SUBROUTINE Agrif_dta_ts
@@ -542,38 +577,46 @@
       !
       ! --- West --- !
-      istart = 2
-      iend   = 1 + nbghostcells
-      DO ji = mi0(istart), mi1(iend)
-         DO jj = 1, jpj
-            ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
+      IF(lk_west) THEN
+         istart = 2
+         iend   = 1 + nbghostcells
+         DO ji = mi0(istart), mi1(iend)
+            DO jj = 1, jpj
+               ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
+            ENDDO
          ENDDO
-      ENDDO
+      ENDIF
       !
       ! --- East --- !
-      istart = jpiglo - nbghostcells
-      iend   = jpiglo - 1
-      DO ji = mi0(istart), mi1(iend)
-         DO jj = 1, jpj
-            ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
+      IF(lk_east) THEN
+         istart = jpiglo - nbghostcells
+         iend   = jpiglo - 1
+         DO ji = mi0(istart), mi1(iend)
+            DO jj = 1, jpj
+               ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
+            ENDDO
          ENDDO
-      ENDDO
+      ENDIF
       !
       ! --- South --- !
-      jstart = 2
-      jend   = 1 + nbghostcells
-      DO jj = mj0(jstart), mj1(jend)
-         DO ji = 1, jpi
-            ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
+      IF(lk_south) THEN
+         jstart = 2
+         jend   = 1 + nbghostcells
+         DO jj = mj0(jstart), mj1(jend)
+            DO ji = 1, jpi
+               ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
+            ENDDO
          ENDDO
-      ENDDO
+      ENDIF
       !
       ! --- North --- !
-      jstart = jpjglo - nbghostcells
-      jend   = jpjglo - 1
-      DO jj = mj0(jstart), mj1(jend)
-         DO ji = 1, jpi
-            ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
+      IF(lk_north) THEN
+         jstart = jpjglo - nbghostcells
+         jend   = jpjglo - 1
+         DO jj = mj0(jstart), mj1(jend)
+            DO ji = 1, jpi
+               ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
+            ENDDO
          ENDDO
-      ENDDO
+      ENDIF
       !
    END SUBROUTINE Agrif_ssh
@@ -593,38 +636,46 @@
       !
       ! --- West --- !
-      istart = 2
-      iend   = 1+nbghostcells
-      DO ji = mi0(istart), mi1(iend)
-         DO jj = 1, jpj
-            ssha_e(ji,jj) = hbdy(ji,jj)
+      IF(lk_west) THEN
+         istart = 2
+         iend   = 1+nbghostcells
+         DO ji = mi0(istart), mi1(iend)
+            DO jj = 1, jpj
+               ssha_e(ji,jj) = hbdy(ji,jj)
+            ENDDO
          ENDDO
-      ENDDO
+      ENDIF
       !
       ! --- East --- !
-      istart = jpiglo - nbghostcells
-      iend   = jpiglo - 1
-      DO ji = mi0(istart), mi1(iend)
-         DO jj = 1, jpj
-            ssha_e(ji,jj) = hbdy(ji,jj)
+      IF(lk_east) THEN
+         istart = jpiglo - nbghostcells
+         iend   = jpiglo - 1
+         DO ji = mi0(istart), mi1(iend)
+            DO jj = 1, jpj
+               ssha_e(ji,jj) = hbdy(ji,jj)
+            ENDDO
          ENDDO
-      ENDDO
+      ENDIF
       !
       ! --- South --- !
-      jstart = 2
-      jend   = 1+nbghostcells
-      DO jj = mj0(jstart), mj1(jend)
-         DO ji = 1, jpi
-            ssha_e(ji,jj) = hbdy(ji,jj)
+      IF(lk_south) THEN
+         jstart = 2
+         jend   = 1+nbghostcells
+         DO jj = mj0(jstart), mj1(jend)
+            DO ji = 1, jpi
+               ssha_e(ji,jj) = hbdy(ji,jj)
+            ENDDO
          ENDDO
-      ENDDO
+      ENDIF
       !
       ! --- North --- !
-      jstart = jpjglo - nbghostcells
-      jend   = jpjglo - 1
-      DO jj = mj0(jstart), mj1(jend)
-         DO ji = 1, jpi
-            ssha_e(ji,jj) = hbdy(ji,jj)
+      IF(lk_north) THEN
+         jstart = jpjglo - nbghostcells
+         jend   = jpjglo - 1
+         DO jj = mj0(jstart), mj1(jend)
+            DO ji = 1, jpi
+               ssha_e(ji,jj) = hbdy(ji,jj)
+            ENDDO
          ENDDO
-      ENDDO
+      ENDIF
       !
    END SUBROUTINE Agrif_ssh_ts
@@ -662,12 +713,17 @@
       INTEGER  ::   ji, jj, jk, jn  ! dummy loop indices
       INTEGER  ::   N_in, N_out
+      INTEGER  :: item
       ! vertical interpolation:
       REAL(wp) :: zhtot
       REAL(wp), DIMENSION(k1:k2,1:jpts) :: tabin
-      REAL(wp), DIMENSION(k1:k2) :: h_in
-      REAL(wp), DIMENSION(1:jpk) :: h_out
-      !!----------------------------------------------------------------------
-
-      IF( before ) THEN         
+      REAL(wp), DIMENSION(k1:k2) :: h_in, z_in
+      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
+      !!----------------------------------------------------------------------
+
+      IF( before ) THEN
+
+         item = Kmm_a
+         IF( l_ini_child )   Kmm_a = Kbb_a  
+
          DO jn = 1,jpts
             DO jk=k1,k2
@@ -678,71 +734,96 @@
               END DO
            END DO
-        END DO
-
-# if defined key_vertical
-        ! Interpolate thicknesses
-        ! Warning: these are masked, hence extrapolated prior interpolation.
-        DO jk=k1,k2
-           DO jj=j1,j2
-              DO ji=i1,i2
-                  ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a)
-              END DO
-           END DO
-        END DO
-
-        ! Extrapolate thicknesses in partial bottom cells:
-        ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
-        IF (ln_zps) THEN
-           DO jj=j1,j2
-              DO ji=i1,i2
-                  jk = mbkt(ji,jj)
-                  ptab(ji,jj,jk,jpts+1) = 0._wp
-              END DO
-           END DO           
-        END IF
-     
-        ! Save ssh at last level:
-        IF (.NOT.ln_linssh) THEN
-           ptab(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 
-        ELSE
-           ptab(i1:i2,j1:j2,k2,jpts+1) = 0._wp
-        END IF      
-# endif
+         END DO
+
+         IF( l_vremap .OR. l_ini_child) THEN
+            ! Interpolate thicknesses
+            ! Warning: these are masked, hence extrapolated prior interpolation.
+            DO jk=k1,k2
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                      ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a)
+
+                  END DO
+               END DO
+            END DO
+
+            ! Extrapolate thicknesses in partial bottom cells:
+            ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
+            IF (ln_zps) THEN
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                      jk = mbkt(ji,jj)
+                      ptab(ji,jj,jk,jpts+1) = 0._wp
+                  END DO
+               END DO           
+            END IF
+        
+            ! Save ssh at last level:
+            IF (.NOT.ln_linssh) THEN
+               ptab(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 
+            ELSE
+               ptab(i1:i2,j1:j2,k2,jpts+1) = 0._wp
+            END IF      
+         ENDIF
+         Kmm_a = item
+
       ELSE 
-
-# if defined key_vertical 
-         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp 
-            
-         DO jj=j1,j2
-            DO ji=i1,i2
-               ts(ji,jj,:,:,Krhs_a) = 0._wp
-               N_in = mbkt_parent(ji,jj)
-               zhtot = 0._wp
-               DO jk=1,N_in !k2 = jpk of parent grid
-                  IF (jk==N_in) THEN
-                     h_in(jk) = ht0_parent(ji,jj) + ptab(ji,jj,k2,n2) - zhtot
-                  ELSE
-                     h_in(jk) = ptab(ji,jj,jk,n2)
+         item = Krhs_a
+         IF( l_ini_child )   Krhs_a = Kbb_a  
+
+         IF( l_vremap .OR. l_ini_child ) THEN
+            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp 
+               
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  ts(ji,jj,:,:,Krhs_a) = 0.                  
+               !   IF( l_ini_child) ts(ji,jj,:,:,Krhs_a) = ptab(ji,jj,:,1:jpts)
+                  N_in = mbkt_parent(ji,jj)
+                  zhtot = 0._wp
+                  DO jk=1,N_in !k2 = jpk of parent grid
+                     IF (jk==N_in) THEN
+                        h_in(jk) = ht0_parent(ji,jj) + ptab(ji,jj,k2,n2) - zhtot
+                     ELSE
+                        h_in(jk) = ptab(ji,jj,jk,n2)
+                     ENDIF
+                     zhtot = zhtot + h_in(jk)
+                     tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)
+                  END DO
+                  z_in(1) = 0.5_wp * h_in(1) - zhtot + ht0_parent(ji,jj)
+                  DO jk=2,N_in
+                     z_in(jk) = z_in(jk-1) + 0.5_wp * h_in(jk)
+                  ENDDO
+
+                  N_out = 0
+                  DO jk=1,jpk ! jpk of child grid
+                     IF (tmask(ji,jj,jk) == 0._wp) EXIT 
+                     N_out = N_out + 1
+                     h_out(jk) = e3t(ji,jj,jk,Krhs_a)
+                  ENDDO
+
+                  z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + ht_0(ji,jj)
+                  DO jk=2,N_out
+                     z_out(jk) = z_out(jk-1) + 0.5_wp * h_out(jk)
+                  ENDDO
+
+                  IF (N_in*N_out > 0) THEN
+                     IF( l_ini_child ) THEN
+                        CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),          &
+                                      &   z_out(1:N_out),N_in,N_out,jpts)  
+                     ELSE 
+                        CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),   &
+                                      &   h_out(1:N_out),N_in,N_out,jpts)  
+                     ENDIF
                   ENDIF
-                  zhtot = zhtot + h_in(jk)
-                  tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)
-               END DO
-               N_out = 0
-               DO jk=1,jpk ! jpk of child grid
-                  IF (tmask(ji,jj,jk) == 0._wp) EXIT 
-                  N_out = N_out + 1
-                  h_out(jk) = e3t(ji,jj,jk,Krhs_a)
                ENDDO
-               IF (N_in*N_out > 0) THEN
-                  CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),h_out(1:N_out),N_in,N_out,jpts)
-               ENDIF
             ENDDO
-         ENDDO
-# else
-         !
-         DO jn=1, jpts
-            ts(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
-         END DO
-# endif
+            Krhs_a = item
+ 
+         ELSE
+         
+            DO jn=1, jpts
+                ts(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
+            END DO
+         ENDIF
 
       ENDIF
@@ -780,85 +861,111 @@
       REAL(wp) :: zrhoy, zhtot
       ! vertical interpolation:
-      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
-      REAL(wp), DIMENSION(1:jpk) :: h_out
-      INTEGER  :: N_in, N_out
+      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in, z_in
+      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
+      INTEGER  :: N_in, N_out,item
       REAL(wp) :: h_diff
       !!---------------------------------------------    
       !
       IF (before) THEN 
+
+         item = Kmm_a
+         IF( l_ini_child )   Kmm_a = Kbb_a     
+
          DO jk=1,jpk
             DO jj=j1,j2
                DO ji=i1,i2
                   ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)*umask(ji,jj,jk)) 
-# if defined key_vertical
-                  ! Interpolate thicknesses (masked for subsequent extrapolation)
-                  ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
-# endif
-               END DO
-            END DO
-         END DO
-# if defined key_vertical
+                  IF( l_vremap .OR. l_ini_child) THEN
+                     ! Interpolate thicknesses (masked for subsequent extrapolation)
+                     ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
+                  ENDIF
+               END DO
+            END DO
+         END DO
+
+        IF( l_vremap .OR. l_ini_child) THEN
          ! Extrapolate thicknesses in partial bottom cells:
          ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
-         IF (ln_zps) THEN
-            DO jj=j1,j2
-               DO ji=i1,i2
-                  jk = mbku(ji,jj)
-                  ptab(ji,jj,jk,2) = 0._wp
-               END DO
-            END DO           
-         END IF
-        ! Save ssh at last level:
-        ptab(i1:i2,j1:j2,k2,2) = 0._wp
-        IF (.NOT.ln_linssh) THEN
-           ! This vertical sum below should be replaced by the sea-level at U-points (optimization):
-           DO jk=1,jpk
-              ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3u(i1:i2,j1:j2,jk,Kmm_a) * umask(i1:i2,j1:j2,jk)
-           END DO
-           ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hu_0(i1:i2,j1:j2)
-        END IF 
-# endif
+            IF (ln_zps) THEN
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                     jk = mbku(ji,jj)
+                     ptab(ji,jj,jk,2) = 0._wp
+                  END DO
+               END DO           
+            END IF
+
+           ! Save ssh at last level:
+           ptab(i1:i2,j1:j2,k2,2) = 0._wp
+           IF (.NOT.ln_linssh) THEN
+              ! This vertical sum below should be replaced by the sea-level at U-points (optimization):
+              DO jk=1,jpk
+                 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3u(i1:i2,j1:j2,jk,Kmm_a) * umask(i1:i2,j1:j2,jk)
+              END DO
+              ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hu_0(i1:i2,j1:j2)
+           END IF
+        ENDIF
+
+         Kmm_a = item
          !
       ELSE
          zrhoy = Agrif_rhoy()
-# if defined key_vertical
+
+        IF( l_vremap .OR. l_ini_child) THEN
 ! VERTICAL REFINEMENT BEGIN
 
-         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
-
-         DO ji=i1,i2
-            DO jj=j1,j2
-               uu(ji,jj,:,Krhs_a) = 0._wp
-               N_in = mbku_parent(ji,jj)
-               zhtot = 0._wp
-               DO jk=1,N_in
-                  IF (jk==N_in) THEN
-                     h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot
-                  ELSE
-                     h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 
-                  ENDIF
-                  zhtot = zhtot + h_in(jk)
-                  tabin(jk) = ptab(ji,jj,jk,1)/(e2u(ji,jj)*zrhoy*h_in(jk))
-              ENDDO
-                  
-              N_out = 0
-              DO jk=1,jpk
-                 if (umask(ji,jj,jk) == 0) EXIT
-                 N_out = N_out + 1
-                 h_out(N_out) = e3u(ji,jj,jk,Krhs_a)
-              ENDDO
-              IF (N_in*N_out > 0) THEN
-                 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1)
-              ENDIF
+            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
+
+            DO ji=i1,i2
+               DO jj=j1,j2
+                  uu(ji,jj,:,Krhs_a) = 0._wp
+                  N_in = mbku_parent(ji,jj)
+                  zhtot = 0._wp
+                  DO jk=1,N_in
+                     IF (jk==N_in) THEN
+                        h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot
+                     ELSE
+                        h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 
+                     ENDIF
+                     zhtot = zhtot + h_in(jk)
+                     IF( h_in(jk) .GT. 0. ) THEN
+                     tabin(jk) = ptab(ji,jj,jk,1)/(e2u(ji,jj)*zrhoy*h_in(jk))
+                     ELSE
+                     tabin(jk) = 0.
+                     ENDIF
+                 ENDDO
+                 z_in(1) = 0.5_wp * h_in(1) - zhtot + hu0_parent(ji,jj) 
+                 DO jk=2,N_in
+                    z_in(jk) = z_in(jk-1) + 0.5_wp * h_in(jk)
+                 ENDDO
+                     
+                 N_out = 0
+                 DO jk=1,jpk
+                    IF (umask(ji,jj,jk) == 0) EXIT
+                    N_out = N_out + 1
+                    h_out(N_out) = e3u(ji,jj,jk,Krhs_a)
+                 ENDDO
+
+                 z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + hu_0(ji,jj)
+                 DO jk=2,N_out
+                    z_out(jk) = z_out(jk-1) + 0.5_wp * h_out(jk) 
+                 ENDDO  
+
+                 IF (N_in*N_out > 0) THEN
+                     IF( l_ini_child ) THEN
+                        CALL remap_linear       (tabin(1:N_in),z_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),z_out(1:N_out),N_in,N_out,1)
+                     ELSE
+                        CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1)
+                     ENDIF   
+                 ENDIF
+               ENDDO
             ENDDO
-         ENDDO
-
-# else
-         DO jk = 1, jpkm1
-            DO jj=j1,j2
-               uu(i1:i2,jj,jk,Krhs_a) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u(i1:i2,jj,jk,Krhs_a) )
-            END DO
-         END DO
-# endif
+         ELSE
+            DO jk = 1, jpkm1
+               DO jj=j1,j2
+                  uu(i1:i2,jj,jk,Krhs_a) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u(i1:i2,jj,jk,Krhs_a) )
+               END DO
+            END DO
+         ENDIF
 
       ENDIF
@@ -878,80 +985,107 @@
       REAL(wp) :: zrhox
       ! vertical interpolation:
-      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
-      REAL(wp), DIMENSION(1:jpk) :: h_out
-      INTEGER  :: N_in, N_out
+      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in, z_in
+      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
+      INTEGER  :: N_in, N_out, item
       REAL(wp) :: h_diff, zhtot
       !!---------------------------------------------    
       !      
-      IF (before) THEN          
+      IF (before) THEN   
+
+         item = Kmm_a
+         IF( l_ini_child )   Kmm_a = Kbb_a     
+       
          DO jk=k1,k2
             DO jj=j1,j2
                DO ji=i1,i2
                   ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)*vmask(ji,jj,jk))
-# if defined key_vertical
-                  ! Interpolate thicknesses (masked for subsequent extrapolation)
-                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
-# endif
-               END DO
-            END DO
-         END DO
-# if defined key_vertical
+                  IF( l_vremap .OR. l_ini_child) THEN
+                     ! Interpolate thicknesses (masked for subsequent extrapolation)
+                     ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
+                  ENDIF
+               END DO
+            END DO
+         END DO
+
+         IF( l_vremap .OR. l_ini_child) THEN
          ! Extrapolate thicknesses in partial bottom cells:
          ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
-         IF (ln_zps) THEN
+            IF (ln_zps) THEN
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                     jk = mbkv(ji,jj)
+                     ptab(ji,jj,jk,2) = 0._wp
+                  END DO
+               END DO           
+            END IF
+            ! Save ssh at last level:
+            ptab(i1:i2,j1:j2,k2,2) = 0._wp
+            IF (.NOT.ln_linssh) THEN
+               ! This vertical sum below should be replaced by the sea-level at V-points (optimization):
+               DO jk=1,jpk
+                  ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3v(i1:i2,j1:j2,jk,Kmm_a) * vmask(i1:i2,j1:j2,jk)
+               END DO
+               ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hv_0(i1:i2,j1:j2)
+            END IF 
+         ENDIF
+         item = Kmm_a
+
+      ELSE       
+         zrhox = Agrif_rhox()
+
+         IF( l_vremap .OR. l_ini_child ) THEN
+
+            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
+
             DO jj=j1,j2
                DO ji=i1,i2
-                  jk = mbkv(ji,jj)
-                  ptab(ji,jj,jk,2) = 0._wp
-               END DO
-            END DO           
-         END IF
-        ! Save ssh at last level:
-        ptab(i1:i2,j1:j2,k2,2) = 0._wp
-        IF (.NOT.ln_linssh) THEN
-           ! This vertical sum below should be replaced by the sea-level at V-points (optimization):
-           DO jk=1,jpk
-              ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3v(i1:i2,j1:j2,jk,Kmm_a) * vmask(i1:i2,j1:j2,jk)
-           END DO
-           ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hv_0(i1:i2,j1:j2)
-        END IF 
-# endif
-      ELSE       
-         zrhox = Agrif_rhox()
-# if defined key_vertical
-
-         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
-
-         DO jj=j1,j2
-            DO ji=i1,i2
-               vv(ji,jj,:,Krhs_a) = 0._wp
-               N_in = mbkv_parent(ji,jj)
-               zhtot = 0._wp
-               DO jk=1,N_in
-                  IF (jk==N_in) THEN
-                     h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot
-                  ELSE
-                     h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 
+                  vv(ji,jj,:,Krhs_a) = 0._wp
+                  N_in = mbkv_parent(ji,jj)
+                  zhtot = 0._wp
+                  DO jk=1,N_in
+                     IF (jk==N_in) THEN
+                        h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot
+                     ELSE
+                        h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 
+                     ENDIF
+                     zhtot = zhtot + h_in(jk)
+                     IF( h_in(jk) .GT. 0. ) THEN
+                       tabin(jk) = ptab(ji,jj,jk,1)/(e1v(ji,jj)*zrhox*h_in(jk))
+                     ELSE
+                       tabin(jk)  = 0.
+                     ENDIF 
+                  ENDDO
+
+                  z_in(1) = 0.5_wp * h_in(1) - zhtot + hv0_parent(ji,jj)
+                  DO jk=2,N_in
+                     z_in(jk) = z_in(jk-1) + 0.5_wp * h_in(jk)
+                  ENDDO
+
+                  N_out = 0
+                  DO jk=1,jpk
+                     IF (vmask(ji,jj,jk) == 0) EXIT
+                     N_out = N_out + 1
+                     h_out(N_out) = e3v(ji,jj,jk,Krhs_a)
+                  ENDDO
+
+                  z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + hv_0(ji,jj)
+                  DO jk=2,N_out
+                     z_out(jk) = z_out(jk-1) + 0.5_wp * h_out(jk)
+                  ENDDO
+ 
+                  IF (N_in*N_out > 0) THEN
+                     IF( l_ini_child ) THEN
+                        CALL remap_linear       (tabin(1:N_in),z_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),z_out(1:N_out),N_in,N_out,1)
+                     ELSE
+                        CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1)
+                     ENDIF   
                   ENDIF
-                  zhtot = zhtot + h_in(jk)
-                  tabin(jk) = ptab(ji,jj,jk,1)/(e1v(ji,jj)*zrhox*h_in(jk))
-              ENDDO
-         
-               N_out = 0
-               DO jk=1,jpk
-                  if (vmask(ji,jj,jk) == 0) EXIT
-                  N_out = N_out + 1
-                  h_out(N_out) = e3v(ji,jj,jk,Krhs_a)
-               END DO
-               IF (N_in*N_out > 0) THEN
-                  call reconstructandremap(tabin(1:N_in),h_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1)
-               ENDIF
-            END DO
-         END DO
-# else
-         DO jk = 1, jpkm1
-            vv(i1:i2,j1:j2,jk,Krhs_a) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Krhs_a) )
-         END DO
-# endif
+               END DO
+            END DO
+         ELSE
+            DO jk = 1, jpkm1
+               vv(i1:i2,j1:j2,jk,Krhs_a) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Krhs_a) )
+            END DO
+         ENDIF
       ENDIF
       !        
@@ -1163,5 +1297,4 @@
    END SUBROUTINE interpe3t
 
-
    SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
       !!----------------------------------------------------------------------
@@ -1185,66 +1318,67 @@
               END DO
            END DO
-        END DO
-
-# if defined key_vertical
-        ! Interpolate thicknesses
-        ! Warning: these are masked, hence extrapolated prior interpolation.
-        DO jk=k1,k2
-           DO jj=j1,j2
-              DO ji=i1,i2
-                  ptab(ji,jj,jk,2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a)
-              END DO
-           END DO
-        END DO
-
-        ! Extrapolate thicknesses in partial bottom cells:
-        ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
-        IF (ln_zps) THEN
-           DO jj=j1,j2
-              DO ji=i1,i2
-                  jk = mbkt(ji,jj)
-                  ptab(ji,jj,jk,2) = 0._wp
-              END DO
-           END DO           
-        END IF
-     
-        ! Save ssh at last level:
-        IF (.NOT.ln_linssh) THEN
-           ptab(i1:i2,j1:j2,k2,2) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 
-        ELSE
-           ptab(i1:i2,j1:j2,k2,2) = 0._wp
-        END IF      
-# endif
+         END DO
+
+         IF( l_vremap ) THEN
+            ! Interpolate thicknesses
+            ! Warning: these are masked, hence extrapolated prior interpolation.
+            DO jk=k1,k2
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                      ptab(ji,jj,jk,2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a)
+                  END DO
+               END DO
+            END DO
+
+            ! Extrapolate thicknesses in partial bottom cells:
+            ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
+            IF (ln_zps) THEN
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                      jk = mbkt(ji,jj)
+                      ptab(ji,jj,jk,2) = 0._wp
+                  END DO
+               END DO           
+            END IF
+        
+           ! Save ssh at last level:
+            IF (.NOT.ln_linssh) THEN
+               ptab(i1:i2,j1:j2,k2,2) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 
+            ELSE
+               ptab(i1:i2,j1:j2,k2,2) = 0._wp
+            END IF      
+          ENDIF
+
       ELSE 
-#ifdef key_vertical         
-         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
-         avm_k(i1:i2,j1:j2,k1:k2) = 0._wp
-            
-         DO jj = j1, j2
-            DO ji =i1, i2
-               N_in = mbkt_parent(ji,jj)
-               IF ( tmask(ji,jj,1) == 0._wp) N_in = 0
-               z_in(N_in+1) = ht0_parent(ji,jj) + ptab(ji,jj,k2,2)
-               DO jk = N_in, 1, -1  ! Parent vertical grid               
-                     z_in(jk) = z_in(jk+1) - ptab(ji,jj,jk,2)
-                    tabin(jk) = ptab(ji,jj,jk,1)
-               END DO
-               N_out = mbkt(ji,jj) 
-               DO jk = 1, N_out        ! Child vertical grid
-                  z_out(jk) = gdepw(ji,jj,jk,Kmm_a)
+
+         IF( l_vremap ) THEN
+            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
+            avm_k(i1:i2,j1:j2,k1:k2) = 0._wp
+               
+            DO jj = j1, j2
+               DO ji =i1, i2
+                  N_in = mbkt_parent(ji,jj)
+                  IF ( tmask(ji,jj,1) == 0._wp) N_in = 0
+                  z_in(N_in+1) = ht0_parent(ji,jj) + ptab(ji,jj,k2,2)
+                  DO jk = N_in, 1, -1  ! Parent vertical grid               
+                        z_in(jk) = z_in(jk+1) - ptab(ji,jj,jk,2)
+                       tabin(jk) = ptab(ji,jj,jk,1)
+                  END DO
+                  N_out = mbkt(ji,jj) 
+                  DO jk = 1, N_out        ! Child vertical grid
+                     z_out(jk) = gdepw(ji,jj,jk,Kmm_a)
+                  ENDDO
+                  IF (N_in*N_out > 0) THEN
+                     CALL remap_linear(tabin(1:N_in),z_in(1:N_in),avm_k(ji,jj,1:N_out),z_out(1:N_out),N_in,N_out,1)
+                  ENDIF
                ENDDO
-               IF (N_in*N_out > 0) THEN
-                  CALL remap_linear(tabin(1:N_in),z_in(1:N_in),avm_k(ji,jj,1:N_out),z_out(1:N_out),N_in,N_out,1)
-               ENDIF
             ENDDO
-         ENDDO
-#else
-         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1)
-#endif
+         ELSE
+            avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1)
+         ENDIF
       ENDIF
       !
    END SUBROUTINE interpavm
 
-# if defined key_vertical
    SUBROUTINE interpmbkt( ptab, i1, i2, j1, j2, before )
       !!----------------------------------------------------------------------
@@ -1282,6 +1416,41 @@
       !
    END SUBROUTINE interpht0
-#endif
-
+
+   SUBROUTINE agrif_initts(tabres,i1,i2,j1,j2,k1,k2,m1,m2,before)
+       INTEGER :: i1, i2, j1, j2, k1, k2, m1, m2
+       REAL(wp):: tabres(i1:i2,j1:j2,k1:k2,m1:m2)
+       LOGICAL :: before
+
+       INTEGER :: jm
+
+       IF (before) THEN
+         DO jm=1,jpts
+             tabres(i1:i2,j1:j2,k1:k2,jm) = ts(i1:i2,j1:j2,k1:k2,jm,Kbb_a)
+         END DO
+       ELSE
+         DO jm=1,jpts
+             ts(i1:i2,j1:j2,k1:k2,jm,Kbb_a)=tabres(i1:i2,j1:j2,k1:k2,jm)
+         END DO
+       ENDIF
+   END SUBROUTINE agrif_initts 
+
+   SUBROUTINE agrif_initssh( ptab, i1, i2, j1, j2, before )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE interpsshn  ***
+      !!----------------------------------------------------------------------  
+      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
+      LOGICAL                         , INTENT(in   ) ::   before
+      !
+      !!----------------------------------------------------------------------  
+      !
+      IF( before) THEN
+         ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kbb_a)
+      ELSE
+         ssh(i1:i2,j1:j2,Kbb_a) = ptab(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1)
+      ENDIF
+      !
+   END SUBROUTINE agrif_initssh
+   
 #else
    !!----------------------------------------------------------------------
Index: NEMO/trunk/src/NST/agrif_oce_sponge.F90
===================================================================
--- NEMO/trunk/src/NST/agrif_oce_sponge.F90	(revision 13215)
+++ NEMO/trunk/src/NST/agrif_oce_sponge.F90	(revision 13216)
@@ -80,4 +80,6 @@
       Agrif_SpecialValue=0.
       Agrif_UseSpecialValue = ln_spc_dyn
+      use_sign_north = .TRUE.
+      sign_north = -1.
       !
       tabspongedone_u = .FALSE.
@@ -90,4 +92,5 @@
       !
       Agrif_UseSpecialValue = .FALSE.
+      use_sign_north = .FALSE.
 #endif
       !
@@ -127,42 +130,51 @@
 
          ! --- West --- !
-         ztabramp(:,:) = 0._wp
-         ind1 = 1+nbghostcells
-         DO ji = mi0(ind1), mi1(ind1)                
-            ztabramp(ji,:) = ssumask(ji,:)
-         END DO
-         !
-         zmskwest(:) = 0._wp
-         zmskwest(1:jpj) = MAXVAL(ztabramp(:,:), dim=1)
+         IF( lk_west) THEN
+            ztabramp(:,:) = 0._wp
+            ind1 = 1+nbghostcells
+            DO ji = mi0(ind1), mi1(ind1)                
+               ztabramp(ji,:) = ssumask(ji,:)
+            END DO
+            !
+            zmskwest(:) = 0._wp
+            zmskwest(1:jpj) = MAXVAL(ztabramp(:,:), dim=1)
+         ENDIF
 
          ! --- East --- !
-         ztabramp(:,:) = 0._wp
-         ind1 = jpiglo - nbghostcells - 1
-         DO ji = mi0(ind1), mi1(ind1)                 
-            ztabramp(ji,:) = ssumask(ji,:)
-         END DO
-         !
-         zmskeast(:) = 0._wp
-         zmskeast(1:jpj) = MAXVAL(ztabramp(:,:), dim=1)
+         IF( lk_east ) THEN
+            ztabramp(:,:) = 0._wp
+            ind1 = jpiglo - nbghostcells - 1
+            DO ji = mi0(ind1), mi1(ind1)                 
+               ztabramp(ji,:) = ssumask(ji,:)
+            END DO
+            !
+            zmskeast(:) = 0._wp
+            zmskeast(1:jpj) = MAXVAL(ztabramp(:,:), dim=1)
+         ENDIF
 
          ! --- South --- !
-         ztabramp(:,:) = 0._wp
-         ind1 = 1+nbghostcells
-         DO jj = mj0(ind1), mj1(ind1)                 
-            ztabramp(:,jj) = ssvmask(:,jj)
-         END DO
-         !
-         zmsksouth(:) = 0._wp
-         zmsksouth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2)
+         IF( lk_south ) THEN
+            ztabramp(:,:) = 0._wp
+            ind1 = 1+nbghostcells
+            DO jj = mj0(ind1), mj1(ind1)                 
+               ztabramp(:,jj) = ssvmask(:,jj)
+            END DO
+            !
+            zmsksouth(:) = 0._wp
+            zmsksouth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2)
+         ENDIF
 
          ! --- North --- !
-         ztabramp(:,:) = 0._wp
-         ind1 = jpjglo - nbghostcells - 1
-         DO jj = mj0(ind1), mj1(ind1)                 
-            ztabramp(:,jj) = ssvmask(:,jj)
-         END DO
-         !
-         zmsknorth(:) = 0._wp
-         zmsknorth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2)
+         IF( lk_north) THEN
+            ztabramp(:,:) = 0._wp
+            ind1 = jpjglo - nbghostcells - 1
+            DO jj = mj0(ind1), mj1(ind1)                 
+               ztabramp(:,jj) = ssvmask(:,jj)
+            END DO
+            !
+            zmsknorth(:) = 0._wp
+            zmsknorth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2)
+         ENDIF
+
          ! JC: SPONGE MASKING TO BE SORTED OUT:
          zmskwest(:)  = 1._wp
@@ -192,75 +204,85 @@
 
          ! --- West --- !
-         ind1 = 1+nbghostcells
-         ind2 = 1+nbghostcells + ispongearea 
-         DO ji = mi0(ind1), mi1(ind2)   
-            DO jj = 1, jpj               
-               ztabramp(ji,jj) = REAL( ind2 - mig(ji) ) * z1_ispongearea * zmskwest(jj)
-            END DO
-         END DO
-
-         ! ghost cells:
-         ind1 = 1
-         ind2 = nbghostcells + 1
-         DO ji = mi0(ind1), mi1(ind2)   
-            DO jj = 1, jpj               
-               ztabramp(ji,jj) = zmskwest(jj)
-            END DO
-         END DO
+         IF(lk_west) THEN
+            ind1 = 1+nbghostcells
+            ind2 = 1+nbghostcells + ispongearea 
+            DO ji = mi0(ind1), mi1(ind2)   
+               DO jj = 1, jpj               
+                  ztabramp(ji,jj) = REAL( ind2 - mig(ji) ) * z1_ispongearea * zmskwest(jj)
+               END DO
+            END DO         
+
+            ! ghost cells:
+            ind1 = 1
+            ind2 = nbghostcells + 1
+            DO ji = mi0(ind1), mi1(ind2)   
+               DO jj = 1, jpj               
+                  ztabramp(ji,jj) = zmskwest(jj)
+               END DO
+            END DO
+         ENDIF
 
          ! --- East --- !
-         ind1 = jpiglo - nbghostcells - ispongearea
-         ind2 = jpiglo - nbghostcells
-         DO ji = mi0(ind1), mi1(ind2)
-            DO jj = 1, jpj
-               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mig(ji) - ind1 ) * z1_ispongearea) * zmskeast(jj)
-            ENDDO
-         END DO
-
-         ! ghost cells:
-         ind1 = jpiglo - nbghostcells
-         ind2 = jpiglo
-         DO ji = mi0(ind1), mi1(ind2)
-            DO jj = 1, jpj
-               ztabramp(ji,jj) = zmskeast(jj)
-            ENDDO
-         END DO
+         IF(lk_east) THEN
+            ind1 = jpiglo - nbghostcells - ispongearea
+            ind2 = jpiglo - nbghostcells
+            DO ji = mi0(ind1), mi1(ind2)
+
+               DO jj = 1, jpj
+                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mig(ji) - ind1 ) * z1_ispongearea) * zmskeast(jj)
+               ENDDO
+            END DO
+
+            ! ghost cells:
+            ind1 = jpiglo - nbghostcells
+            ind2 = jpiglo
+            DO ji = mi0(ind1), mi1(ind2)
+
+               DO jj = 1, jpj
+                  ztabramp(ji,jj) = zmskeast(jj)
+               ENDDO
+            END DO
+         ENDIF
 
          ! --- South --- !
-         ind1 = 1+nbghostcells
-         ind2 = 1+nbghostcells + jspongearea
-         DO jj = mj0(ind1), mj1(ind2) 
-            DO ji = 1, jpi
-               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - mjg(jj) ) * z1_jspongearea) * zmsksouth(ji)
-            END DO
-         END DO
-
-         ! ghost cells:
-         ind1 = 1
-         ind2 = nbghostcells + 1
-         DO jj = mj0(ind1), mj1(ind2) 
-            DO ji = 1, jpi
-               ztabramp(ji,jj) = zmsksouth(ji)
-            END DO
-         END DO
+         IF( lk_south ) THEN 
+            ind1 = 1+nbghostcells
+            ind2 = 1+nbghostcells + jspongearea
+            DO jj = mj0(ind1), mj1(ind2) 
+               DO ji = 1, jpi
+                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - mjg(jj) ) * z1_jspongearea) * zmsksouth(ji)
+               END DO
+            END DO
+
+            ! ghost cells:
+            ind1 = 1
+            ind2 = nbghostcells + 1
+            DO jj = mj0(ind1), mj1(ind2) 
+               DO ji = 1, jpi
+                  ztabramp(ji,jj) = zmsksouth(ji)
+               END DO
+            END DO
+         ENDIF
 
          ! --- North --- !
-         ind1 = jpjglo - nbghostcells - jspongearea
-         ind2 = jpjglo - nbghostcells
-         DO jj = mj0(ind1), mj1(ind2)
-            DO ji = 1, jpi
-               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mjg(jj) - ind1 ) * z1_jspongearea) * zmsknorth(ji)
-            END DO
-         END DO
-
-         ! ghost cells:
-         ind1 = jpjglo - nbghostcells
-         ind2 = jpjglo
-         DO jj = mj0(ind1), mj1(ind2)
-            DO ji = 1, jpi
-               ztabramp(ji,jj) = zmsknorth(ji)
-            END DO
-         END DO
-
+         IF( lk_north ) THEN  
+            ind1 = jpjglo - nbghostcells - jspongearea
+            ind2 = jpjglo - nbghostcells
+            DO jj = mj0(ind1), mj1(ind2)
+               DO ji = 1, jpi
+                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mjg(jj) - ind1 ) * z1_jspongearea) * zmsknorth(ji)
+               END DO
+            END DO
+
+            ! ghost cells:
+            ind1 = jpjglo - nbghostcells
+            ind2 = jpjglo
+            DO jj = mj0(ind1), mj1(ind2)
+               DO ji = 1, jpi
+                  ztabramp(ji,jj) = zmsknorth(ji)
+               END DO
+            END DO
+         ENDIF
+      
       ENDIF
 
@@ -334,5 +356,5 @@
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
       INTEGER  ::   iku, ikv
-      REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot, ztrelax
+      REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot
       REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv
       REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff
@@ -438,9 +460,4 @@
          ENDDO
 
-         !* set relaxation time scale
-         IF( l_1st_euler .AND. lk_agrif_fstep ) THEN   ;   ztrelax =   rn_trelax_tra  / (        rn_Dt )
-         ELSE                                          ;   ztrelax =   rn_trelax_tra  / (2._wp * rn_Dt )
-         ENDIF
-
          DO jn = 1, jpts            
             DO jk = 1, jpkm1
@@ -448,5 +465,5 @@
                DO jj = j1,j2
                   DO ji = i1,i2-1
-                     zabe1 = rn_sponge_tra * fspu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
+                     zabe1 = rn_sponge_tra * r1_Dt * fspu(ji,jj) * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
                      ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
                   END DO
@@ -455,5 +472,5 @@
                DO ji = i1,i2
                   DO jj = j1,j2-1
-                     zabe2 = rn_sponge_tra * fspv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
+                     zabe2 = rn_sponge_tra * r1_Dt * fspv(ji,jj) * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
                      ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
                   END DO
@@ -480,5 +497,5 @@
                         ! horizontal diffusive trends
                         ztsa = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) &
-                             &  - ztrelax * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn) 
+                             &  - rn_trelax_tra * r1_Dt * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn) 
                         ! add it to the general tracer trends
                         ts(ji,jj,jk,jn,Krhs_a) = ts(ji,jj,jk,jn,Krhs_a) + ztsa
@@ -507,5 +524,5 @@
 
       ! sponge parameters 
-      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax
+      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot
       REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff
       REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
@@ -595,8 +612,4 @@
          ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
 #endif
-         !* set relaxation time scale
-         IF( l_1st_euler .AND. lk_agrif_fstep ) THEN   ;   ztrelax =   rn_trelax_dyn  / (        rn_Dt )
-         ELSE                                          ;   ztrelax =   rn_trelax_dyn  / (2._wp * rn_Dt )
-         ENDIF
          !
          DO jk = 1, jpkm1                                 ! Horizontal slab
@@ -608,5 +621,5 @@
             DO jj = j1,j2
                DO ji = i1+1,i2   ! vector opt.
-                  zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb_a) * rn_sponge_dyn * fspt(ji,jj)
+                  zbtr = rn_sponge_dyn * r1_Dt * fspt(ji,jj) / e3t(ji,jj,jk,Kbb_a)
                   hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kbb_a) * ubdiff(ji  ,jj,jk) &
                                      &   -e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb_a) * ubdiff(ji-1,jj,jk) ) * zbtr
@@ -616,5 +629,5 @@
             DO jj = j1,j2-1
                DO ji = i1,i2   ! vector opt.
-                  zbtr = r1_e1e2f(ji,jj) * e3f(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj)
+                  zbtr = rn_sponge_dyn * r1_Dt * fspf(ji,jj) * e3f(ji,jj,jk) 
                   rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk)   &
                                     &   +e1u(ji,jj  ) * ubdiff(ji,jj  ,jk) ) * fmask(ji,jj,jk) * zbtr 
@@ -633,5 +646,5 @@
                      zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) )   &
                          & + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj) & 
-                         & - ztrelax  * fspu(ji,jj) * ubdiff(ji,jj,jk)
+                         & - rn_trelax_dyn * r1_Dt * fspu(ji,jj) * ubdiff(ji,jj,jk)
 
                      ! add it to the general momentum trends
@@ -646,4 +659,5 @@
 
          jmax = j2-1
+        ! IF (lk_north) jmax = MIN(jmax,nlcj-nbghostcells-2)   ! North
          IF ((nbondj == 1).OR.(nbondj == 2)) jmax = MIN(jmax,nlcj-nbghostcells-2)   ! North
 
@@ -684,5 +698,5 @@
       !
       INTEGER  ::   ji, jj, jk, imax
-      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax
+      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr, zhtot
       REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff
       REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
@@ -771,8 +785,4 @@
          vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
 # endif
-         !* set relaxation time scale
-         IF( l_1st_euler .AND. lk_agrif_fstep ) THEN   ;   ztrelax =   rn_trelax_dyn  / (        rn_Dt )
-         ELSE                                          ;   ztrelax =   rn_trelax_dyn  / (2._wp * rn_Dt )
-         ENDIF
          !
          DO jk = 1, jpkm1                                 ! Horizontal slab
@@ -784,5 +794,5 @@
             DO jj = j1+1,j2
                DO ji = i1,i2   ! vector opt.
-                  zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb_a) * rn_sponge_dyn * fspt(ji,jj)
+                  zbtr = rn_sponge_dyn * r1_Dt * fspt(ji,jj) / e3t(ji,jj,jk,Kbb_a)
                   hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kbb_a) * vbdiff(ji,jj  ,jk)  &
                                      &  -e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kbb_a) * vbdiff(ji,jj-1,jk)  ) * zbtr
@@ -791,5 +801,5 @@
             DO jj = j1,j2
                DO ji = i1,i2-1   ! vector opt.
-                  zbtr = r1_e1e2f(ji,jj) * e3f(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj)
+                  zbtr = rn_sponge_dyn * r1_Dt * fspf(ji,jj) * e3f(ji,jj,jk) 
                   rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) & 
                                     &  -e2v(ji  ,jj) * vbdiff(ji  ,jj,jk)  ) * fmask(ji,jj,jk) * zbtr
@@ -802,4 +812,5 @@
 
          imax = i2 - 1
+      !   IF(lk_east) imax = MIN(imax,nlci-nbghostcells-2)   ! East
          IF ((nbondi == 1).OR.(nbondi == 2))   imax = MIN(imax,nlci-nbghostcells-2)   ! East
 
@@ -808,5 +819,5 @@
                IF( .NOT. tabspongedone_u(ji,jj) ) THEN
                   DO jk = 1, jpkm1
-                     uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a)                                                               &
+                     uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a)                                                     &
                         & - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) )  &
                         & + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj  ,jk)) * r1_e1u(ji,jj)
@@ -822,8 +833,8 @@
                IF( .NOT. tabspongedone_v(ji,jj) ) THEN
                   DO jk = 1, jpkm1
-                     vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a)                                                                  &
+                     vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a)                                                        &
                         &  + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) )   &
-                        &  + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) * r1_e2v(ji,jj)                      &
-                        &  - ztrelax * fspv(ji,jj) * vbdiff(ji,jj,jk)
+                        &  + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) * r1_e2v(ji,jj)                          &
+                        &  - rn_trelax_dyn * r1_Dt * fspv(ji,jj) * vbdiff(ji,jj,jk)
                   END DO
                ENDIF
Index: NEMO/trunk/src/NST/agrif_oce_update.F90
===================================================================
--- NEMO/trunk/src/NST/agrif_oce_update.F90	(revision 13215)
+++ NEMO/trunk/src/NST/agrif_oce_update.F90	(revision 13216)
@@ -26,4 +26,5 @@
    USE domvvl         ! Need interpolation routines 
    USE vremap         ! Vertical remapping
+   USE lbclnk 
 
    IMPLICIT NONE
@@ -85,4 +86,8 @@
       Agrif_UseSpecialValueInUpdate = .FALSE.
       Agrif_SpecialValueFineGrid = 0.
+
+      use_sign_north = .TRUE.
+      sign_north = -1.
+
       !     
 # if ! defined DECAL_FEEDBACK
@@ -127,4 +132,6 @@
       END IF
       !
+      use_sign_north = .FALSE.
+      !
    END SUBROUTINE Agrif_Update_Dyn
 
@@ -148,4 +155,6 @@
 #  if defined VOL_REFLUX
       IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN
+         use_sign_north = .TRUE.
+         sign_north = -1.
          ! Refluxing on ssh:
 #  if defined DECAL_FEEDBACK_2D
@@ -156,4 +165,5 @@
          CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/ 0, 0/),locupdate2=(/-1,-1/),procname = reflux_sshv)
 #  endif
+         use_sign_north = .FALSE.
       END IF
 #  endif
@@ -826,5 +836,5 @@
    SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
       !!---------------------------------------------
-      !!           *** ROUTINE correct_u_bdy ***
+      !!           *** ROUTINE correct_v_bdy ***
       !!---------------------------------------------
       INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
Index: NEMO/trunk/src/NST/agrif_top_interp.F90
===================================================================
--- NEMO/trunk/src/NST/agrif_top_interp.F90	(revision 13215)
+++ NEMO/trunk/src/NST/agrif_top_interp.F90	(revision 13216)
@@ -119,5 +119,4 @@
             tr(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab_child(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
          END DO
-
       ENDIF
       !
Index: NEMO/trunk/src/NST/agrif_user.F90
===================================================================
--- NEMO/trunk/src/NST/agrif_user.F90	(revision 13215)
+++ NEMO/trunk/src/NST/agrif_user.F90	(revision 13216)
@@ -28,6 +28,4 @@
       !
       !                    !* Agrif initialization
-      CALL agrif_nemo_init
-      CALL Agrif_InitValues_cont_dom
       CALL Agrif_InitValues_cont
 # if defined key_top
@@ -40,41 +38,152 @@
    END SUBROUTINE Agrif_initvalues
 
-   SUBROUTINE Agrif_InitValues_cont_dom
-      !!----------------------------------------------------------------------
-      !!                 *** ROUTINE Agrif_InitValues_cont_dom ***
-      !!----------------------------------------------------------------------
-      !
-      CALL agrif_declare_var_dom
-      !
-   END SUBROUTINE Agrif_InitValues_cont_dom
-
-   SUBROUTINE agrif_declare_var_dom
-      !!----------------------------------------------------------------------
-      !!                 *** ROUTINE agrif_declare_var_dom ***
-      !!----------------------------------------------------------------------
-      USE par_oce, ONLY:  nbghostcells      
+   SUBROUTINE agrif_istate( Kbb, Kmm, Kaa )
+
+       USE domvvl
+       USE domain
+       USE par_oce
+       USE agrif_oce
+       USE agrif_oce_interp
+       USE oce
+       USE lib_mpp
+       USe lbclnk
+
+      INTEGER, INTENT(in)  :: Kbb, Kmm, Kaa
+      INTEGER :: jn
+
+      IF(lwp) WRITE(numout,*) ' '
+      IF(lwp) WRITE(numout,*) 'AGRIF: interp child initial state from parent'
+      IF(lwp) WRITE(numout,*) ' '
+
+      l_ini_child = .TRUE.
+      Agrif_SpecialValue    = 0._wp
+      Agrif_UseSpecialValue = .TRUE.
+      uu(:,:,:,:) = 0.  ;  vv(:,:,:,:) = 0.   ;  ts(:,:,:,:,:) = 0.
+       
+      Krhs_a = Kbb ; Kmm_a = Kbb
+
+      ! Brutal fix to pas 1x1 refinment. 
+  !    IF(Agrif_Irhox() == 1) THEN
+  !       CALL Agrif_Init_Variable(tsini_id, procname=agrif_initts)
+  !    ELSE
+      CALL Agrif_Init_Variable(tsini_id, procname=interptsn)
+
+  !    ENDIF
+! just for VORTEX because Parent velocities can actually be exactly zero
+!      Agrif_UseSpecialValue = .FALSE.
+      Agrif_UseSpecialValue = ln_spc_dyn
+      use_sign_north = .TRUE.
+      sign_north = -1.
+      CALL Agrif_Init_Variable(uini_id , procname=interpun )
+      CALL Agrif_Init_Variable(vini_id , procname=interpvn )
+      use_sign_north = .FALSE.
+
+      Agrif_UseSpecialValue = .FALSE.            !
+      l_ini_child = .FALSE.
+
+      Krhs_a = Kaa ; Kmm_a = Kmm
+
+      DO jn = 1, jpts
+         ts(:,:,:,jn,Kbb) = ts(:,:,:,jn,Kbb)*tmask(:,:,:)
+      END DO
+      uu(:,:,:,Kbb) =  uu(:,:,:,Kbb) * umask(:,:,:)     
+      vv(:,:,:,Kbb) =  vv(:,:,:,Kbb) * vmask(:,:,:) 
+
+
+      CALL lbc_lnk_multi( 'agrif_istate', uu(:,:,:,Kbb), 'U', -1. , vv(:,:,:,Kbb), 'V', -1. )
+      CALL lbc_lnk( 'agrif_istate', ts(:,:,:,:,Kbb), 'T', 1. )
+
+   END SUBROUTINE agrif_istate   
+
+   SUBROUTINE agrif_declare_var_ini
+      !!----------------------------------------------------------------------
+      !!                 *** ROUTINE agrif_declare_var ***
+      !!----------------------------------------------------------------------
+      USE agrif_util
+      USE agrif_oce
+      USE par_oce
+      USE zdf_oce 
+      USE oce
+      USE dom_oce
       !
       IMPLICIT NONE
       !
       INTEGER :: ind1, ind2, ind3
-      !!----------------------------------------------------------------------
+      External :: nemo_mapping
+      !!----------------------------------------------------------------------
+
+! In case of East-West periodicity, prevent AGRIF interpolation at east and west boundaries
+! The procnames will not be called at these boundaries
+      IF (jperio == 1) THEN
+         CALL Agrif_Set_NearCommonBorderX(.TRUE.)
+         CALL Agrif_Set_DistantCommonBorderX(.TRUE.)
+      ENDIF
+
+      IF ( .NOT. lk_south ) THEN
+         CALL Agrif_Set_NearCommonBorderY(.TRUE.)
+      ENDIF
 
       ! 1. Declaration of the type of variable which have to be interpolated
       !---------------------------------------------------------------------
       ind1 =     nbghostcells
-      ind2 = 1 + nbghostcells
-      ind3 = 2 + nbghostcells
-      CALL agrif_declare_variable((/1,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),e1u_id)
-      CALL agrif_declare_variable((/2,1/),(/ind3,ind2/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),e2v_id)
-
+      ind2 = 2 + nbghostcells_x
+      ind3 = 2 + nbghostcells_y_s
+
+      CALL agrif_declare_variable((/2,2,0/),(/ind2,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),e3t_id)
+      CALL agrif_declare_variable((/2,2/)  ,(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),mbkt_id)
+      CALL agrif_declare_variable((/2,2/)  ,(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),ht0_id)
+
+      CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),e1u_id)
+      CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),e2v_id)
+
+   
+      ! Initial or restart velues
+      
+      CALL agrif_declare_variable((/2,2,0,0/),(/ind2  ,ind3  ,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsini_id)
+      CALL agrif_declare_variable((/1,2,0,0/),(/ind2-1,ind3  ,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/)     ,uini_id ) 
+      CALL agrif_declare_variable((/2,1,0,0/),(/ind2  ,ind3-1,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/)     ,vini_id )
+      CALL agrif_declare_variable((/2,2/)    ,(/ind2,ind3/)        ,(/'x','y'/),(/1,1/),(/nlci,nlcj/),sshini_id)
+      ! 
+     
       ! 2. Type of interpolation
       !-------------------------
+      CALL Agrif_Set_bcinterp(e3t_id,interp=AGRIF_constant)
+
+      CALL Agrif_Set_bcinterp(mbkt_id,interp=AGRIF_constant)
+      CALL Agrif_Set_interp  (mbkt_id,interp=AGRIF_constant)
+      CALL Agrif_Set_bcinterp(ht0_id ,interp=AGRIF_constant)
+      CALL Agrif_Set_interp  (ht0_id ,interp=AGRIF_constant)
+
       CALL Agrif_Set_bcinterp( e1u_id, interp1=Agrif_linear, interp2=AGRIF_ppm    )
       CALL Agrif_Set_bcinterp( e2v_id, interp1=AGRIF_ppm   , interp2=Agrif_linear )
 
-      ! 3. Location of interpolation
+      ! Initial fields
+      CALL Agrif_Set_bcinterp(tsini_id ,interp=AGRIF_linear)
+      CALL Agrif_Set_interp  (tsini_id ,interp=AGRIF_linear)
+      CALL Agrif_Set_bcinterp(uini_id  ,interp=AGRIF_linear)
+      CALL Agrif_Set_interp  (uini_id  ,interp=AGRIF_linear)
+      CALL Agrif_Set_bcinterp(vini_id  ,interp=AGRIF_linear)
+      CALL Agrif_Set_interp  (vini_id  ,interp=AGRIF_linear)
+      CALL Agrif_Set_bcinterp(sshini_id,interp=AGRIF_linear)
+      CALL Agrif_Set_interp  (sshini_id,interp=AGRIF_linear)
+
+       ! 3. Location of interpolation
       !-----------------------------
+!      CALL Agrif_Set_bc(  e3t_id, (/-nn_sponge_len*Agrif_irhox(),ind1-1/) )  
+! JC: check near the boundary only until matching in sponge has been sorted out:
+      CALL Agrif_Set_bc(  e3t_id, (/0,ind1-1/) )  
+
+      ! extend the interpolation zone by 1 more point than necessary:
+      ! RB check here
+      CALL Agrif_Set_bc(  mbkt_id, (/-nn_sponge_len*Agrif_irhox()-2,ind1/) )
+      CALL Agrif_Set_bc(  ht0_id,  (/-nn_sponge_len*Agrif_irhox()-2,ind1/) )
+      
       CALL Agrif_Set_bc(e1u_id,(/0,ind1-1/))
-      CALL Agrif_Set_bc(e2v_id,(/0,ind1-1/))
+      CALL Agrif_Set_bc(e2v_id,(/0,ind1-1/))  
+
+      CALL Agrif_Set_bc( tsini_id , (/0,ind1-1/) ) ! if west,  rhox=3 and nbghost=3: columns 2 to 4
+      CALL Agrif_Set_bc( uini_id  , (/0,ind1-1/) ) 
+      CALL Agrif_Set_bc( vini_id  , (/0,ind1-1/) )
+      CALL Agrif_Set_bc( sshini_id, (/0,ind1-1/) )
 
       ! 4. Update type
@@ -87,37 +196,48 @@
       CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Average, update2=Agrif_Update_Copy)
 #endif
-
-   END SUBROUTINE agrif_declare_var_dom
-
-   SUBROUTINE Agrif_InitValues_cont
-      !!----------------------------------------------------------------------
-      !!                 *** ROUTINE Agrif_InitValues_cont ***
-      !!----------------------------------------------------------------------
-      USE agrif_oce
+      
+   !   CALL Agrif_Set_ExternalMapping(nemo_mapping)
+      !
+   END SUBROUTINE agrif_declare_var_ini
+
+
+   SUBROUTINE Agrif_Init_Domain( Kbb, Kmm, Kaa ) 
+      !!----------------------------------------------------------------------
+      !!                 *** ROUTINE Agrif_InitValues_cont_dom ***
+      !!----------------------------------------------------------------------
+  
+         !!----------------------------------------------------------------------
+         !!                 *** ROUTINE Agrif_InitValues_cont ***
+         !!
+         !! ** Purpose ::   Declaration of variables to be interpolated
+         !!----------------------------------------------------------------------
+      USE agrif_oce_update
       USE agrif_oce_interp
       USE agrif_oce_sponge
+      USE Agrif_Util
+      USE oce 
       USE dom_oce
-      USE oce
+      USE zdf_oce
+      USE nemogcm
+      USE agrif_oce
+      !
+      USE lbclnk
       USE lib_mpp
-      USE lbclnk
+      USE in_out_manager
       !
       IMPLICIT NONE
       !
-      INTEGER :: ji, jj
+      INTEGER, INTENT(in) ::  Kbb, Kmm, Kaa
+      !
       LOGICAL :: check_namelist
       CHARACTER(len=15) :: cl_check1, cl_check2, cl_check3, cl_check4 
-#if defined key_vertical
       REAL(wp), DIMENSION(jpi,jpj) ::   zk   ! workspace
-#endif
-      !!----------------------------------------------------------------------
-
-      ! 1. Declaration of the type of variable which have to be interpolated
-      !---------------------------------------------------------------------
-      CALL agrif_declare_var
-
-      ! 2. First interpolations of potentially non zero fields
-      !-------------------------------------------------------
-
-#if defined key_vertical
+      INTEGER :: ji, jj, jk
+      !!----------------------------------------------------------------------
+    
+     ! CALL Agrif_Declare_Var_ini
+
+      IF( agrif_oce_alloc()  > 0 )   CALL ctl_warn('agrif agrif_oce_alloc: allocation of arrays failed')
+
       ! Build consistent parent bathymetry and number of levels
       ! on the child grid 
@@ -126,6 +246,8 @@
       mbkt_parent(:,:) = 0
       !
-      CALL Agrif_Bc_variable(ht0_id ,calledweight=1.,procname=interpht0 )
-      CALL Agrif_Bc_variable(mbkt_id,calledweight=1.,procname=interpmbkt)
+  !    CALL Agrif_Bc_variable(ht0_id ,calledweight=1.,procname=interpht0 )
+  !    CALL Agrif_Bc_variable(mbkt_id,calledweight=1.,procname=interpmbkt)
+      CALL Agrif_Init_Variable(ht0_id , procname=interpht0 )
+      CALL Agrif_Init_Variable(mbkt_id, procname=interpmbkt)
       !
       ! Assume step wise change of bathymetry near interface
@@ -149,12 +271,98 @@
       ENDIF
       !
-      CALL lbc_lnk( 'Agrif_InitValues_cont', hu0_parent, 'U', 1. )
-      CALL lbc_lnk( 'Agrif_InitValues_cont', hv0_parent, 'V', 1. )
+      CALL lbc_lnk( 'Agrif_Init_Domain', hu0_parent, 'U', 1. )
+      CALL lbc_lnk( 'Agrif_Init_Domain', hv0_parent, 'V', 1. )
       zk(:,:) = REAL( mbku_parent(:,:), wp )   ;   CALL lbc_lnk( 'Agrif_InitValues_cont', zk, 'U', 1. )
-      mbku_parent(:,:) = MAX( NINT( zk(:,:) ), 1 )
+      mbku_parent(:,:) = MAX( NINT( zk(:,:) ), 1 ) ;
       zk(:,:) = REAL( mbkv_parent(:,:), wp )   ;   CALL lbc_lnk( 'Agrif_InitValues_cont', zk, 'V', 1. )
       mbkv_parent(:,:) = MAX( NINT( zk(:,:) ), 1 )   
-#endif
-
+
+      IF ( ln_init_chfrpar ) THEN 
+         CALL Agrif_Init_Variable(sshini_id, procname=agrif_initssh)
+         CALL lbc_lnk( 'Agrif_Init_Domain', ssh(:,:,Kbb), 'T', 1. )
+         DO jk = 1, jpk
+               e3t(:,:,jk,Kbb) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb)  ) &
+                        &             / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
+                        &              + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
+         END DO
+      ENDIF
+
+      ! check if masks and bathymetries match
+      IF(ln_chk_bathy) THEN
+         Agrif_UseSpecialValue = .FALSE.
+         !
+         IF(lwp) WRITE(numout,*) ' '
+         IF(lwp) WRITE(numout,*) 'AGRIF: Check Bathymetry and masks near bdys. Level: ', Agrif_Level()
+         !
+         kindic_agr = 0
+         IF( .NOT. l_vremap ) THEN
+            !
+            ! check if tmask and vertical scale factors agree with parent in sponge area:
+            CALL Agrif_Bc_variable(e3t_id,calledweight=1.,procname=interpe3t)
+            !
+         ELSE
+            !
+            ! In case of vertical interpolation, check only that total depths agree between child and parent:
+            DO ji = 1, jpi
+               DO jj = 1, jpj
+                  IF ((mbkt_parent(ji,jj)/=0).AND.(ABS(ht0_parent(ji,jj)-ht_0(ji,jj))>1.e-3)) kindic_agr = kindic_agr + 1
+                  IF ((mbku_parent(ji,jj)/=0).AND.(ABS(hu0_parent(ji,jj)-hu_0(ji,jj))>1.e-3)) kindic_agr = kindic_agr + 1
+                  IF ((mbkv_parent(ji,jj)/=0).AND.(ABS(hv0_parent(ji,jj)-hv_0(ji,jj))>1.e-3)) kindic_agr = kindic_agr + 1
+               END DO
+            END DO
+
+            CALL mpp_sum( 'agrif_user', kindic_agr )
+            IF( kindic_agr /= 0 ) THEN
+               CALL ctl_stop('==> Child Bathymetry is NOT correct near boundaries.')
+            ELSE
+               IF(lwp) WRITE(numout,*) '==> Child Bathymetry is ok near boundaries.'
+               IF(lwp) WRITE(numout,*) ' '
+            ENDIF  
+         ENDIF
+      ENDIF
+
+      IF( l_vremap ) THEN
+      ! Additional constrain that should be removed someday:
+         IF ( Agrif_Parent(jpk).GT.jpk ) THEN
+            CALL ctl_stop( ' With l_vremap, child grids must have jpk greater or equal to the parent value' )
+         ENDIF
+      ENDIF
+      !
+   END SUBROUTINE Agrif_Init_Domain
+
+
+   SUBROUTINE Agrif_InitValues_cont
+         !!----------------------------------------------------------------------
+         !!                 *** ROUTINE Agrif_InitValues_cont ***
+         !!
+         !! ** Purpose ::   Declaration of variables to be interpolated
+         !!----------------------------------------------------------------------
+      USE agrif_oce_update
+      USE agrif_oce_interp
+      USE agrif_oce_sponge
+      USE Agrif_Util
+      USE oce 
+      USE dom_oce
+      USE zdf_oce
+      USE nemogcm
+      USE agrif_oce
+      !
+      USE lbclnk
+      USE lib_mpp
+      USE in_out_manager
+      !
+      IMPLICIT NONE
+      !
+      LOGICAL :: check_namelist
+      CHARACTER(len=15) :: cl_check1, cl_check2, cl_check3, cl_check4 
+      REAL(wp), DIMENSION(jpi,jpj) ::   zk   ! workspace
+      INTEGER :: ji, jj
+
+      ! 1. Declaration of the type of variable which have to be interpolated
+      !---------------------------------------------------------------------
+      CALL agrif_declare_var
+
+      ! 2. First interpolations of potentially non zero fields
+      !-------------------------------------------------------
       Agrif_SpecialValue    = 0._wp
       Agrif_UseSpecialValue = .TRUE.
@@ -163,8 +371,10 @@
       tabspongedone_tsn = .FALSE.
       CALL Agrif_Bc_variable(tsn_sponge_id,calledweight=1.,procname=interptsn_sponge)
-      ! reset ts(:,:,:,:,Krhs_a) to zero
+      ! reset tsa to zero
       ts(:,:,:,:,Krhs_a) = 0._wp
 
       Agrif_UseSpecialValue = ln_spc_dyn
+      use_sign_north = .TRUE.
+      sign_north = -1.
       CALL Agrif_Bc_variable(un_interp_id,calledweight=1.,procname=interpun)
       CALL Agrif_Bc_variable(vn_interp_id,calledweight=1.,procname=interpvn)
@@ -175,4 +385,5 @@
       tabspongedone_v = .FALSE.
       CALL Agrif_Bc_variable(vn_sponge_id,calledweight=1.,procname=interpvn_sponge)
+      use_sign_north = .FALSE.
       uu(:,:,:,Krhs_a) = 0._wp
       vv(:,:,:,Krhs_a) = 0._wp
@@ -185,43 +396,20 @@
       IF ( ln_dynspg_ts ) THEN
          Agrif_UseSpecialValue = ln_spc_dyn
+         use_sign_north = .TRUE.
+         sign_north = -1.
          CALL Agrif_Bc_variable(unb_id,calledweight=1.,procname=interpunb)
          CALL Agrif_Bc_variable(vnb_id,calledweight=1.,procname=interpvnb)
          CALL Agrif_Bc_variable(ub2b_interp_id,calledweight=1.,procname=interpub2b)
          CALL Agrif_Bc_variable(vb2b_interp_id,calledweight=1.,procname=interpvb2b)
+         use_sign_north = .FALSE.
          ubdy(:,:) = 0._wp
          vbdy(:,:) = 0._wp
       ENDIF
-
-      Agrif_UseSpecialValue = .FALSE.
-
-      ! 3. Some controls
+      Agrif_UseSpecialValue = .FALSE. 
+
       !-----------------
       check_namelist = .TRUE.
 
       IF( check_namelist ) THEN 
-
-         ! Check time steps           
-         IF( NINT(Agrif_Rhot()) * NINT(rn_Dt) .NE. Agrif_Parent(rn_Dt) ) THEN
-            WRITE(cl_check1,*)  NINT(Agrif_Parent(rn_Dt))
-            WRITE(cl_check2,*)  NINT(rn_Dt)
-            WRITE(cl_check3,*)  NINT(Agrif_Parent(rn_Dt)/Agrif_Rhot())
-            CALL ctl_stop( 'Incompatible time step between ocean grids',   &
-                  &               'parent grid value : '//cl_check1    ,   & 
-                  &               'child  grid value : '//cl_check2    ,   & 
-                  &               'value on child grid should be changed to : '//cl_check3 )
-         ENDIF
-
-         ! Check run length
-         IF( Agrif_IRhot() * (Agrif_Parent(nitend)- &
-               Agrif_Parent(nit000)+1) .NE. (nitend-nit000+1) ) THEN
-            WRITE(cl_check1,*)  (Agrif_Parent(nit000)-1)*Agrif_IRhot() + 1
-            WRITE(cl_check2,*)   Agrif_Parent(nitend)   *Agrif_IRhot()
-            CALL ctl_warn( 'Incompatible run length between grids'                      ,   &
-                  &               'nit000 on fine grid will be changed to : '//cl_check1,   &
-                  &               'nitend on fine grid will be changed to : '//cl_check2    )
-            nit000 = (Agrif_Parent(nit000)-1)*Agrif_IRhot() + 1
-            nitend =  Agrif_Parent(nitend)   *Agrif_IRhot()
-         ENDIF
-
          ! Check free surface scheme
          IF ( ( Agrif_Parent(ln_dynspg_ts ).AND.ln_dynspg_exp ).OR.&
@@ -251,48 +439,6 @@
             STOP
          ENDIF
-
-      ENDIF
-
-      ! check if masks and bathymetries match
-      IF(ln_chk_bathy) THEN
-         Agrif_UseSpecialValue = .FALSE.
-         !
-         IF(lwp) WRITE(numout,*) ' '
-         IF(lwp) WRITE(numout,*) 'AGRIF: Check Bathymetry and masks near bdys. Level: ', Agrif_Level()
-         !
-         kindic_agr = 0
-# if ! defined key_vertical
-         !
-         ! check if tmask and vertical scale factors agree with parent in sponge area:
-         CALL Agrif_Bc_variable(e3t_id,calledweight=1.,procname=interpe3t)
-         !
-# else
-         !
-         ! In case of vertical interpolation, check only that total depths agree between child and parent:
-         DO ji = 1, jpi
-            DO jj = 1, jpj
-               IF ((mbkt_parent(ji,jj)/=0).AND.(ABS(ht0_parent(ji,jj)-ht_0(ji,jj))>1.e-3)) kindic_agr = kindic_agr + 1
-               IF ((mbku_parent(ji,jj)/=0).AND.(ABS(hu0_parent(ji,jj)-hu_0(ji,jj))>1.e-3)) kindic_agr = kindic_agr + 1
-               IF ((mbkv_parent(ji,jj)/=0).AND.(ABS(hv0_parent(ji,jj)-hv_0(ji,jj))>1.e-3)) kindic_agr = kindic_agr + 1
-            END DO
-         END DO
-# endif
-         CALL mpp_sum( 'agrif_user', kindic_agr )
-         IF( kindic_agr /= 0 ) THEN
-            CALL ctl_stop('==> Child Bathymetry is NOT correct near boundaries.')
-         ELSE
-            IF(lwp) WRITE(numout,*) '==> Child Bathymetry is ok near boundaries.'
-            IF(lwp) WRITE(numout,*) ' '
-         END IF  
-         !    
-      ENDIF
-
-# if defined key_vertical
-      ! Additional constrain that should be removed someday:
-      IF ( Agrif_Parent(jpk).GT.jpk ) THEN
-	 CALL ctl_stop( ' With key_vertical, child grids must have jpk greater or equal to the parent value' )
-      ENDIF
-# endif
-      ! 
+      ENDIF
+
    END SUBROUTINE Agrif_InitValues_cont
 
@@ -314,65 +460,57 @@
       ! 1. Declaration of the type of variable which have to be interpolated
       !---------------------------------------------------------------------
+
       ind1 =     nbghostcells
-      ind2 = 1 + nbghostcells
-      ind3 = 2 + nbghostcells
+      ind2 = 2 + nbghostcells_x
+      ind3 = 2 + nbghostcells_y_s
+
 # if defined key_vertical
-      CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_id)
-      CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_sponge_id)
-
-      CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_interp_id)
-      CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_interp_id)
-      CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_update_id)
-      CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_update_id)
-      CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_sponge_id)
-      CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_sponge_id)
+      CALL agrif_declare_variable((/2,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_id)
+      CALL agrif_declare_variable((/2,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_sponge_id)
+      CALL agrif_declare_variable((/1,2,0,0/),(/ind2-1,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_interp_id)
+      CALL agrif_declare_variable((/2,1,0,0/),(/ind2,ind3-1,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_interp_id)
+      CALL agrif_declare_variable((/1,2,0,0/),(/ind2-1,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_update_id)
+      CALL agrif_declare_variable((/2,1,0,0/),(/ind2,ind3-1,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_update_id)
+      CALL agrif_declare_variable((/1,2,0,0/),(/ind2-1,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_sponge_id)
+      CALL agrif_declare_variable((/2,1,0,0/),(/ind2,ind3-1,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_sponge_id)
 # else
-      CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_id)
-      CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_sponge_id)
-
-      CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_interp_id)
-      CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_interp_id)
-      CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_update_id)
-      CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_update_id)
-      CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_sponge_id)
-      CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_sponge_id)
+      CALL agrif_declare_variable((/2,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_id)
+      CALL agrif_declare_variable((/2,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_sponge_id)
+      CALL agrif_declare_variable((/1,2,0,0/),(/ind2-1,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_interp_id)
+      CALL agrif_declare_variable((/2,1,0,0/),(/ind2,ind3-1,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_interp_id)
+      CALL agrif_declare_variable((/1,2,0,0/),(/ind2-1,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_update_id)
+      CALL agrif_declare_variable((/2,1,0,0/),(/ind2,ind3-1,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_update_id)
+      CALL agrif_declare_variable((/1,2,0,0/),(/ind2-1,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_sponge_id)
+      CALL agrif_declare_variable((/2,1,0,0/),(/ind2,ind3-1,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_sponge_id)
 # endif
-
-      CALL agrif_declare_variable((/2,2,0/),(/ind3,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),e3t_id)
-
-# if defined key_vertical
-      CALL agrif_declare_variable((/2,2/),(/ind3,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),mbkt_id)
-      CALL agrif_declare_variable((/2,2/),(/ind3,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),ht0_id)
-# endif
-
-      CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,3/),scales_t_id)
-
-      CALL agrif_declare_variable((/1,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),unb_id)
-      CALL agrif_declare_variable((/2,1/),(/ind3,ind2/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),vnb_id)
-      CALL agrif_declare_variable((/1,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),ub2b_interp_id)
-      CALL agrif_declare_variable((/2,1/),(/ind3,ind2/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),vb2b_interp_id)
-      CALL agrif_declare_variable((/1,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),ub2b_update_id)
-      CALL agrif_declare_variable((/2,1/),(/ind3,ind2/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),vb2b_update_id)
-
-      CALL agrif_declare_variable((/2,2/),(/ind3,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),sshn_id)
-
-      IF( ln_zdftke.OR.ln_zdfgls ) THEN
+      CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),unb_id)
+      CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),vnb_id)
+      CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),ub2b_interp_id)
+      CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),vb2b_interp_id)
+      CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),ub2b_update_id)
+      CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),vb2b_update_id)
+
+      CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),sshn_id)
+
+
+      IF( ln_zdftke.OR.ln_zdfgls ) THEN  ! logical not known at this point
 !         CALL agrif_declare_variable((/2,2,0/),(/ind3,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/), en_id)
 !         CALL agrif_declare_variable((/2,2,0/),(/ind3,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),avt_id)
 # if defined key_vertical
-         CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),avm_id)
+         CALL agrif_declare_variable((/2,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),avm_id)
 # else
-         CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),avm_id)
+         CALL agrif_declare_variable((/2,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),avm_id)
 # endif
       ENDIF
-
+     
       ! 2. Type of interpolation
       !-------------------------
       CALL Agrif_Set_bcinterp(tsn_id,interp=AGRIF_linear)
-
       CALL Agrif_Set_bcinterp(un_interp_id,interp1=Agrif_linear,interp2=AGRIF_ppm)
       CALL Agrif_Set_bcinterp(vn_interp_id,interp1=AGRIF_ppm,interp2=Agrif_linear)
 
       CALL Agrif_Set_bcinterp(tsn_sponge_id,interp=AGRIF_linear)
+      CALL Agrif_Set_bcinterp(un_sponge_id,interp1=Agrif_linear,interp2=AGRIF_ppm)
+      CALL Agrif_Set_bcinterp(vn_sponge_id,interp1=AGRIF_ppm,interp2=Agrif_linear)
 
       CALL Agrif_Set_bcinterp(sshn_id,interp=AGRIF_linear)
@@ -390,17 +528,8 @@
 !<
 
-      CALL Agrif_Set_bcinterp(un_sponge_id,interp1=Agrif_linear,interp2=AGRIF_ppm)
-      CALL Agrif_Set_bcinterp(vn_sponge_id,interp1=AGRIF_ppm,interp2=Agrif_linear)
-
-      CALL Agrif_Set_bcinterp(e3t_id,interp=AGRIF_constant)
-
-# if defined key_vertical
-      CALL Agrif_Set_bcinterp(mbkt_id,interp=AGRIF_constant)
-      CALL Agrif_Set_bcinterp(ht0_id ,interp=AGRIF_constant)
-# endif
-
-      IF( ln_zdftke.OR.ln_zdfgls )   CALL Agrif_Set_bcinterp( avm_id, interp=AGRIF_linear )
-
-      ! 3. Location of interpolation
+      IF( ln_zdftke.OR.ln_zdfgls )  CALL Agrif_Set_bcinterp( avm_id, interp=AGRIF_linear )
+    
+
+       ! 3. Location of interpolation
       !-----------------------------
       CALL Agrif_Set_bc(       tsn_id, (/0,ind1-1/) ) ! if west,  rhox=3 and nbghost=3: columns 2 to 4
@@ -418,19 +547,8 @@
       CALL Agrif_Set_bc( vb2b_interp_id, (/0,ind1-1/) )
 
-!      CALL Agrif_Set_bc(  e3t_id, (/-nn_sponge_len*Agrif_irhox(),ind1-1/) )  
-! JC: check near the boundary only until matching in sponge has been sorted out:
-      CALL Agrif_Set_bc(  e3t_id, (/0,ind1-1/) )  
-
-# if defined key_vertical 
-      ! extend the interpolation zone by 1 more point than necessary:
-      CALL Agrif_Set_bc(  mbkt_id, (/-nn_sponge_len*Agrif_irhox()-2,ind1/) )
-      CALL Agrif_Set_bc(  ht0_id,  (/-nn_sponge_len*Agrif_irhox()-2,ind1/) )
-# endif
-
-      IF( ln_zdftke.OR.ln_zdfgls )   CALL Agrif_Set_bc( avm_id, (/0,ind1/) )
+      IF( ln_zdftke.OR.ln_zdfgls ) CALL Agrif_Set_bc( avm_id, (/0,ind1/) )
 
       ! 4. Update type
       !--------------- 
-      CALL Agrif_Set_Updatetype(scales_t_id, update = AGRIF_Update_Average)
 
 # if defined UPD_HIGH
@@ -444,9 +562,9 @@
       CALL Agrif_Set_Updatetype(e3t_id, update = Agrif_Update_Full_Weighting)
 
-      IF( ln_zdftke.OR.ln_zdfgls ) THEN
+  !    IF( ln_zdftke.OR.ln_zdfgls ) THEN
 !         CALL Agrif_Set_Updatetype( en_id, update = AGRIF_Update_Full_Weighting)
 !         CALL Agrif_Set_Updatetype(avt_id, update = AGRIF_Update_Full_Weighting)
 !         CALL Agrif_Set_Updatetype(avm_id, update = AGRIF_Update_Full_Weighting)
-      ENDIF
+   !   ENDIF
 
 #else
@@ -460,9 +578,9 @@
       CALL Agrif_Set_Updatetype(e3t_id, update = AGRIF_Update_Average)
 
-      IF( ln_zdftke.OR.ln_zdfgls ) THEN
+ !     IF( ln_zdftke.OR.ln_zdfgls ) THEN
 !         CALL Agrif_Set_Updatetype( en_id, update = AGRIF_Update_Average)
 !         CALL Agrif_Set_Updatetype(avt_id, update = AGRIF_Update_Average)
 !         CALL Agrif_Set_Updatetype(avm_id, update = AGRIF_Update_Average)
-      ENDIF
+ !     ENDIF
 
 #endif
@@ -472,7 +590,4 @@
 #if defined key_si3
 SUBROUTINE Agrif_InitValues_cont_ice
-      !!----------------------------------------------------------------------
-      !!                 *** ROUTINE Agrif_InitValues_cont_ice ***
-      !!----------------------------------------------------------------------
       USE Agrif_Util
       USE sbc_oce, ONLY : nn_fsbc  ! clem: necessary otherwise Agrif_Parent(nn_fsbc) = nn_fsbc
@@ -482,11 +597,7 @@
       USE agrif_ice_interp
       USE lib_mpp
-      !
-      IMPLICIT NONE
-      !!----------------------------------------------------------------------
-      !
-      ! Declaration of the type of variable which have to be interpolated (parent=>child)
-      !----------------------------------------------------------------------------------
-      CALL agrif_declare_var_ice
+      !!----------------------------------------------------------------------
+      !!                 *** ROUTINE Agrif_InitValues_cont_ice ***
+      !!----------------------------------------------------------------------
 
       ! Controls
@@ -495,5 +606,5 @@
       !          the run must satisfy CFL=Uice/(dx/dt) < 0.6/nn_fsbc(child)
       !          therefore, if nn_fsbc(child)>1 one must reduce the time-step in proportion to nn_fsbc(child), which is not acceptable
-      !       If a solution is found, the following stop could be removed because the rest of the code take nn_fsbc(child) into account
+      !       If a solution is found, the following stop could be removed because the rest of the code take nn_fsbc(child) into account     
       IF( nn_fsbc > 1 )  CALL ctl_stop('nn_fsbc(child) must be set to 1 otherwise agrif and sea-ice may not work properly')
 
@@ -516,12 +627,13 @@
       !!                 *** ROUTINE agrif_declare_var_ice ***
       !!----------------------------------------------------------------------
+
       USE Agrif_Util
       USE ice
-      USE par_oce, ONLY : nbghostcells
+      USE par_oce, ONLY : nbghostcells, nbghostcells_x, nbghostcells_y_s
       !
       IMPLICIT NONE
       !
       INTEGER :: ind1, ind2, ind3
-      !!----------------------------------------------------------------------
+         !!----------------------------------------------------------------------
       !
       ! 1. Declaration of the type of variable which have to be interpolated (parent=>child)
@@ -532,10 +644,15 @@
       !                            2,2 = two ghost lines
       !-------------------------------------------------------------------------------------
+
       ind1 =     nbghostcells
-      ind2 = 1 + nbghostcells
-      ind3 = 2 + nbghostcells
-      CALL agrif_declare_variable((/2,2,0/),(/ind3,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpl*(8+nlay_s+nlay_i)/),tra_ice_id)
-      CALL agrif_declare_variable((/1,2/)  ,(/ind2,ind3/)  ,(/'x','y'/)    ,(/1,1/)  ,(/nlci,nlcj/)                      ,u_ice_id  )
-      CALL agrif_declare_variable((/2,1/)  ,(/ind3,ind2/)  ,(/'x','y'/)    ,(/1,1/)  ,(/nlci,nlcj/)                      ,v_ice_id  )
+      ind2 = 2 + nbghostcells_x
+      ind3 = 2 + nbghostcells_y_s
+      CALL agrif_declare_variable((/2,2,0/),(/ind2,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpl*(8+nlay_s+nlay_i)/),tra_ice_id)
+      CALL agrif_declare_variable((/1,2/)  ,(/ind2-1,ind3/)  ,(/'x','y'/)    ,(/1,1/)  ,(/nlci,nlcj/)                      ,u_ice_id  )
+      CALL agrif_declare_variable((/2,1/)  ,(/ind2,ind3-1/)  ,(/'x','y'/)    ,(/1,1/)  ,(/nlci,nlcj/)                      ,v_ice_id  )
+
+      CALL agrif_declare_variable((/2,2,0/),(/ind3,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpl*(8+nlay_s+nlay_i)/),tra_iceini_id)
+      CALL agrif_declare_variable((/1,2/)  ,(/ind2-1,ind3/)  ,(/'x','y'/)    ,(/1,1/)  ,(/nlci,nlcj/)                      ,u_iceini_id  )
+      CALL agrif_declare_variable((/2,1/)  ,(/ind2,ind3-1/)  ,(/'x','y'/)    ,(/1,1/)  ,(/nlci,nlcj/)                      ,v_iceini_id  )
 
       ! 2. Set interpolations (normal & tangent to the grid cell for velocities)
@@ -545,4 +662,11 @@
       CALL Agrif_Set_bcinterp(v_ice_id  , interp1 = AGRIF_ppm   ,interp2 = Agrif_linear)
 
+      CALL Agrif_Set_bcinterp(tra_iceini_id, interp  = AGRIF_linear)
+      CALL Agrif_Set_interp  (tra_iceini_id, interp  = AGRIF_linear)
+      CALL Agrif_Set_bcinterp(u_iceini_id  , interp  = AGRIF_linear  )
+      CALL Agrif_Set_interp  (u_iceini_id  , interp  = AGRIF_linear   )
+      CALL Agrif_Set_bcinterp(v_iceini_id  , interp  = AGRIF_linear)
+      CALL Agrif_Set_interp  (v_iceini_id  , interp  = AGRIF_linear)
+
       ! 3. Set location of interpolations
       !----------------------------------
@@ -550,4 +674,8 @@
       CALL Agrif_Set_bc(u_ice_id  ,(/0,ind1/))
       CALL Agrif_Set_bc(v_ice_id  ,(/0,ind1/))
+
+      CALL Agrif_Set_bc(tra_iceini_id,(/0,ind1/))
+      CALL Agrif_Set_bc(u_iceini_id  ,(/0,ind1/))
+      CALL Agrif_Set_bc(v_iceini_id  ,(/0,ind1/))
 
       ! 4. Set update type in case 2 ways (child=>parent) (normal & tangent to the grid cell for velocities)
@@ -557,9 +685,9 @@
       CALL Agrif_Set_Updatetype(u_ice_id  , update1 = Agrif_Update_Average       , update2 = Agrif_Update_Full_Weighting)
       CALL Agrif_Set_Updatetype(v_ice_id  , update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average       )
-#else
+# else
       CALL Agrif_Set_Updatetype(tra_ice_id, update  = AGRIF_Update_Average)
       CALL Agrif_Set_Updatetype(u_ice_id  , update1 = Agrif_Update_Copy   , update2 = Agrif_Update_Average)
       CALL Agrif_Set_Updatetype(v_ice_id  , update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy   )
-#endif
+# endif
 
    END SUBROUTINE agrif_declare_var_ice
@@ -585,36 +713,39 @@
       USE agrif_top_sponge
       !!
-      IMPLICIT NONE
-      !
-      CHARACTER(len=10) :: cl_check1, cl_check2, cl_check3
-      LOGICAL :: check_namelist
-      !!----------------------------------------------------------------------
-
-      ! 1. Declaration of the type of variable which have to be interpolated
-      !---------------------------------------------------------------------
-      CALL agrif_declare_var_top
-
-      ! 2. First interpolations of potentially non zero fields
-      !-------------------------------------------------------
-      Agrif_SpecialValue=0._wp
-      Agrif_UseSpecialValue = .TRUE.
-      CALL Agrif_Bc_variable(trn_id,calledweight=1.,procname=interptrn)
-      Agrif_UseSpecialValue = .FALSE.
-      CALL Agrif_Sponge
-      tabspongedone_trn = .FALSE.
-      CALL Agrif_Bc_variable(trn_sponge_id,calledweight=1.,procname=interptrn_sponge)
-      ! reset ts(:,:,:,:,Krhs_a) to zero
-      tr(:,:,:,:,Krhs_a) = 0._wp
-
-      ! 3. Some controls
-      !-----------------
-      check_namelist = .TRUE.
-
-      IF( check_namelist ) THEN
-         ! Check time steps
-      IF( NINT(Agrif_Rhot()) * NINT(rn_Dt) .NE. Agrif_Parent(rn_Dt) ) THEN
-         WRITE(cl_check1,*)  Agrif_Parent(rn_Dt)
-         WRITE(cl_check2,*)  rn_Dt
-         WRITE(cl_check3,*)  rn_Dt*Agrif_Rhot()
+  
+   !!
+   IMPLICIT NONE
+   !
+   CHARACTER(len=10) :: cl_check1, cl_check2, cl_check3
+   LOGICAL :: check_namelist
+      !!----------------------------------------------------------------------
+
+
+   ! 1. Declaration of the type of variable which have to be interpolated
+   !---------------------------------------------------------------------
+   CALL agrif_declare_var_top
+
+   ! 2. First interpolations of potentially non zero fields
+   !-------------------------------------------------------
+   Agrif_SpecialValue=0.
+   Agrif_UseSpecialValue = .TRUE.
+   CALL Agrif_Bc_variable(trn_id,calledweight=1.,procname=interptrn)
+   Agrif_UseSpecialValue = .FALSE.
+   CALL Agrif_Sponge
+   tabspongedone_trn = .FALSE.
+   CALL Agrif_Bc_variable(trn_sponge_id,calledweight=1.,procname=interptrn_sponge)
+   ! reset tsa to zero
+   tra(:,:,:,:) = 0.
+
+   ! 3. Some controls
+   !-----------------
+   check_namelist = .TRUE.
+
+   IF( check_namelist ) THEN
+      ! Check time steps
+      IF( NINT(Agrif_Rhot()) * NINT(rdt) .NE. Agrif_Parent(rdt) ) THEN
+         WRITE(cl_check1,*)  Agrif_Parent(rdt)
+         WRITE(cl_check2,*)  rdt
+         WRITE(cl_check3,*)  rdt*Agrif_Rhot()
          CALL ctl_stop( 'incompatible time step between grids',   &
                &               'parent grid value : '//cl_check1    ,   & 
@@ -635,5 +766,4 @@
          nitend =  Agrif_Parent(nitend)   *Agrif_IRhot()
       ENDIF
-
    ENDIF
    !
@@ -655,15 +785,19 @@
       !!----------------------------------------------------------------------
 
+
+
+!RB_CMEMS : declare here init for top      
       ! 1. Declaration of the type of variable which have to be interpolated
       !---------------------------------------------------------------------
       ind1 =     nbghostcells
-      ind2 = 1 + nbghostcells
-      ind3 = 2 + nbghostcells
+      ind2 = 2 + nbghostcells_x
+      ind3 = 2 + nbghostcells_y_s
 # if defined key_vertical
-      CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra+1/),trn_id)
-      CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra+1/),trn_sponge_id)
+      CALL agrif_declare_variable((/2,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra+1/),trn_id)
+      CALL agrif_declare_variable((/2,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra+1/),trn_sponge_id)
 # else
-      CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra/),trn_id)
-      CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra/),trn_sponge_id)
+! LAURENT: STRANGE why (3,3) here ?
+      CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra/),trn_id)
+      CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra/),trn_sponge_id)
 # endif
 
@@ -705,13 +839,14 @@
       !!                     *** ROUTINE agrif_init ***
       !!----------------------------------------------------------------------
-      USE agrif_oce 
-      USE agrif_ice
-      USE in_out_manager
-      USE lib_mpp
+   USE agrif_oce 
+   USE agrif_ice
+   USE dom_oce
+   USE in_out_manager
+   USE lib_mpp
       !!
       IMPLICIT NONE
       !
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namagrif/ ln_agrif_2way, rn_sponge_tra, rn_sponge_dyn, rn_trelax_tra, rn_trelax_dyn, &
+      NAMELIST/namagrif/ ln_agrif_2way, ln_init_chfrpar, rn_sponge_tra, rn_sponge_dyn, rn_trelax_tra, rn_trelax_dyn, &
                        & ln_spc_dyn, ln_chk_bathy
       !!--------------------------------------------------------------------------------------
@@ -729,19 +864,37 @@
          WRITE(numout,*) '   Namelist namagrif : set AGRIF parameters'
          WRITE(numout,*) '      Two way nesting activated ln_agrif_2way         = ', ln_agrif_2way
-         WRITE(numout,*) '      sponge coefficient for tracers    rn_sponge_tra = ', rn_sponge_tra, ' m^2/s'
-         WRITE(numout,*) '      sponge coefficient for dynamics   rn_sponge_tra = ', rn_sponge_dyn, ' m^2/s'
-         WRITE(numout,*) '      time relaxation for tracers       rn_trelax_tra = ', rn_trelax_tra, ' ad.'
-         WRITE(numout,*) '      time relaxation for dynamics      rn_trelax_dyn = ', rn_trelax_dyn, ' ad.'
+         WRITE(numout,*) '      child initial state from parent ln_init_chfrpar = ', ln_init_chfrpar
+         WRITE(numout,*) '      ad. sponge coeft for tracers      rn_sponge_tra = ', rn_sponge_tra
+         WRITE(numout,*) '      ad. sponge coeft for dynamics     rn_sponge_tra = ', rn_sponge_dyn
+         WRITE(numout,*) '      ad. time relaxation for tracers   rn_trelax_tra = ', rn_trelax_tra
+         WRITE(numout,*) '      ad. time relaxation for dynamics  rn_trelax_dyn = ', rn_trelax_dyn
          WRITE(numout,*) '      use special values for dynamics   ln_spc_dyn    = ', ln_spc_dyn
          WRITE(numout,*) '      check bathymetry                  ln_chk_bathy  = ', ln_chk_bathy
       ENDIF
-      !
-      !
-      IF( agrif_oce_alloc()  > 0 )   CALL ctl_warn('agrif agrif_oce_alloc: allocation of arrays failed')
+
+      lk_west  = .NOT. ( Agrif_Ix() == 1 )
+      lk_east  = .NOT. ( Agrif_Ix() + nbcellsx/AGRIF_Irhox() == Agrif_Parent(jpiglo) -1 )
+      lk_south = .NOT. ( Agrif_Iy() == 1 )
+      lk_north = .NOT. ( Agrif_Iy() + nbcellsy/AGRIF_Irhoy() == Agrif_Parent(jpjglo) -1 )
+
+      !
+      ! Set the number of ghost cells according to periodicity
+      nbghostcells_x = nbghostcells
+      nbghostcells_y_s = nbghostcells
+      nbghostcells_y_n = nbghostcells
+      !
+      IF ( jperio == 1 ) nbghostcells_x = 0
+      IF ( .NOT. lk_south ) nbghostcells_y_s = 0
+
+      ! Some checks
+      IF( jpiglo /= nbcellsx + 2 + 2*nbghostcells_x )   &
+          CALL ctl_stop( 'STOP', 'agrif_nemo_init: Agrif children requires jpiglo == nbcellsx + 2 + 2*nbghostcells_x' )
+      IF( jpjglo /= nbcellsy + 2 + nbghostcells_y_s + nbghostcells_y_n )   &
+          CALL ctl_stop( 'STOP', 'agrif_nemo_init: Agrif children requires jpjglo == nbcellsy + 2 + nbghostcells_y_s + nbghostcells_y_n' )
+      IF( ln_use_jattr )   CALL ctl_stop( 'STOP', 'agrif_nemo_init:Agrif children requires ln_use_jattr = .false. ' )
       !
    END SUBROUTINE agrif_nemo_init
 
 # if defined key_mpp_mpi
-
    SUBROUTINE Agrif_InvLoc( indloc, nprocloc, i, indglob )
       !!----------------------------------------------------------------------
@@ -803,4 +956,231 @@
 # endif
 
+   SUBROUTINE nemo_mapping(ndim,ptx,pty,bounds,bounds_chunks,correction_required,nb_chunks)
+      !!----------------------------------------------------------------------
+      !!                   *** ROUTINE Nemo_mapping ***
+      !!----------------------------------------------------------------------
+      USE dom_oce
+      !!
+      IMPLICIT NONE
+      !
+      INTEGER :: ndim
+      INTEGER :: ptx, pty
+      INTEGER, DIMENSION(ndim,2,2) :: bounds
+      INTEGER, DIMENSION(:,:,:,:), ALLOCATABLE :: bounds_chunks
+      LOGICAL, DIMENSION(:), ALLOCATABLE :: correction_required
+      INTEGER :: nb_chunks
+      !
+      INTEGER :: i
+
+      IF (agrif_debug_interp) THEN
+         DO i=1,ndim
+            WRITE(*,*) 'direction = ',i,bounds(i,1,2),bounds(i,2,2)
+         ENDDO
+      ENDIF
+
+      IF( bounds(2,2,2) > jpjglo) THEN
+         IF( bounds(2,1,2) <=jpjglo) THEN
+            nb_chunks = 2
+            ALLOCATE(bounds_chunks(nb_chunks,ndim,2,2))
+            ALLOCATE(correction_required(nb_chunks))
+            DO i = 1,nb_chunks
+               bounds_chunks(i,:,:,:) = bounds
+            END DO
+        
+      ! FIRST CHUNCK (for j<=jpjglo)   
+
+      ! Original indices
+            bounds_chunks(1,1,1,1) = bounds(1,1,2)
+            bounds_chunks(1,1,2,1) = bounds(1,2,2)
+            bounds_chunks(1,2,1,1) = bounds(2,1,2)
+            bounds_chunks(1,2,2,1) = jpjglo
+
+            bounds_chunks(1,1,1,2) = bounds(1,1,2)
+            bounds_chunks(1,1,2,2) = bounds(1,2,2)
+            bounds_chunks(1,2,1,2) = bounds(2,1,2)
+            bounds_chunks(1,2,2,2) = jpjglo
+
+      ! Correction required or not
+            correction_required(1)=.FALSE.
+       
+      ! SECOND CHUNCK (for j>jpjglo)
+
+      ! Original indices
+            bounds_chunks(2,1,1,1) = bounds(1,1,2)
+            bounds_chunks(2,1,2,1) = bounds(1,2,2)
+            bounds_chunks(2,2,1,1) = jpjglo-2
+            bounds_chunks(2,2,2,1) = bounds(2,2,2)
+
+      ! Where to find them
+      ! We use the relation TAB(ji,jj)=TAB(jpiglo-ji+2,jpjglo-2-(jj-jpjglo))
+
+            IF( ptx == 2) THEN ! T, V points
+               bounds_chunks(2,1,1,2) = jpiglo-bounds(1,2,2)+2
+               bounds_chunks(2,1,2,2) = jpiglo-bounds(1,1,2)+2
+            ELSE ! U, F points
+               bounds_chunks(2,1,1,2) = jpiglo-bounds(1,2,2)+1
+               bounds_chunks(2,1,2,2) = jpiglo-bounds(1,1,2)+1       
+            ENDIF
+
+            IF( pty == 2) THEN ! T, U points
+               bounds_chunks(2,2,1,2) = jpjglo-2-(bounds(2,2,2) -jpjglo)
+               bounds_chunks(2,2,2,2) = jpjglo-2-(jpjglo-2      -jpjglo)
+            ELSE ! V, F points
+               bounds_chunks(2,2,1,2) = jpjglo-3-(bounds(2,2,2) -jpjglo)
+               bounds_chunks(2,2,2,2) = jpjglo-3-(jpjglo-2      -jpjglo)
+            ENDIF
+      ! Correction required or not
+            correction_required(2)=.TRUE.
+
+         ELSE
+            nb_chunks = 1
+            ALLOCATE(bounds_chunks(nb_chunks,ndim,2,2))
+            ALLOCATE(correction_required(nb_chunks))
+            DO i=1,nb_chunks
+               bounds_chunks(i,:,:,:) = bounds
+            END DO
+
+            bounds_chunks(1,1,1,1) = bounds(1,1,2)
+            bounds_chunks(1,1,2,1) = bounds(1,2,2)
+            bounds_chunks(1,2,1,1) = bounds(2,1,2)
+            bounds_chunks(1,2,2,1) = bounds(2,2,2)
+
+            bounds_chunks(1,1,1,2) = jpiglo-bounds(1,2,2)+2
+            bounds_chunks(1,1,2,2) = jpiglo-bounds(1,1,2)+2
+
+            bounds_chunks(1,2,1,2) = jpjglo-2-(bounds(2,2,2)-jpjglo)
+            bounds_chunks(1,2,2,2) = jpjglo-2-(bounds(2,1,2)-jpjglo)
+
+            IF( ptx == 2) THEN ! T, V points
+               bounds_chunks(1,1,1,2) = jpiglo-bounds(1,2,2)+2
+               bounds_chunks(1,1,2,2) = jpiglo-bounds(1,1,2)+2
+            ELSE ! U, F points
+               bounds_chunks(1,1,1,2) = jpiglo-bounds(1,2,2)+1
+               bounds_chunks(1,1,2,2) = jpiglo-bounds(1,1,2)+1    	
+            ENDIF
+
+            IF (pty == 2) THEN ! T, U points
+               bounds_chunks(1,2,1,2) = jpjglo-2-(bounds(2,2,2) -jpjglo)
+               bounds_chunks(1,2,2,2) = jpjglo-2-(bounds(2,1,2) -jpjglo)
+            ELSE ! V, F points
+               bounds_chunks(1,2,1,2) = jpjglo-3-(bounds(2,2,2) -jpjglo)
+               bounds_chunks(1,2,2,2) = jpjglo-3-(bounds(2,1,2) -jpjglo)
+            ENDIF
+
+            correction_required(1)=.TRUE.          
+         ENDIF
+
+      ELSE IF (bounds(1,1,2) < 1) THEN
+         IF (bounds(1,2,2) > 0) THEN
+            nb_chunks = 2
+            ALLOCATE(correction_required(nb_chunks))
+            correction_required=.FALSE.
+            ALLOCATE(bounds_chunks(nb_chunks,ndim,2,2))
+            DO i=1,nb_chunks
+               bounds_chunks(i,:,:,:) = bounds
+            END DO
+              
+            bounds_chunks(1,1,1,2) = bounds(1,1,2)+jpiglo-2
+            bounds_chunks(1,1,2,2) = 1+jpiglo-2
+          
+            bounds_chunks(1,1,1,1) = bounds(1,1,2)
+            bounds_chunks(1,1,2,1) = 1
+       
+            bounds_chunks(2,1,1,2) = 2
+            bounds_chunks(2,1,2,2) = bounds(1,2,2)
+          
+            bounds_chunks(2,1,1,1) = 2
+            bounds_chunks(2,1,2,1) = bounds(1,2,2)
+
+         ELSE
+            nb_chunks = 1
+            ALLOCATE(correction_required(nb_chunks))
+            correction_required=.FALSE.
+            ALLOCATE(bounds_chunks(nb_chunks,ndim,2,2))
+            DO i=1,nb_chunks
+               bounds_chunks(i,:,:,:) = bounds
+            END DO    
+            bounds_chunks(1,1,1,2) = bounds(1,1,2)+jpiglo-2
+            bounds_chunks(1,1,2,2) = bounds(1,2,2)+jpiglo-2
+          
+            bounds_chunks(1,1,1,1) = bounds(1,1,2)
+           bounds_chunks(1,1,2,1) = bounds(1,2,2)
+         ENDIF
+      ELSE
+         nb_chunks=1  
+         ALLOCATE(correction_required(nb_chunks))
+         correction_required=.FALSE.
+         ALLOCATE(bounds_chunks(nb_chunks,ndim,2,2))
+         DO i=1,nb_chunks
+            bounds_chunks(i,:,:,:) = bounds
+         END DO
+         bounds_chunks(1,1,1,2) = bounds(1,1,2)
+         bounds_chunks(1,1,2,2) = bounds(1,2,2)
+         bounds_chunks(1,2,1,2) = bounds(2,1,2)
+         bounds_chunks(1,2,2,2) = bounds(2,2,2)
+          
+         bounds_chunks(1,1,1,1) = bounds(1,1,2)
+         bounds_chunks(1,1,2,1) = bounds(1,2,2)
+         bounds_chunks(1,2,1,1) = bounds(2,1,2)
+         bounds_chunks(1,2,2,1) = bounds(2,2,2)              
+      ENDIF
+        
+   END SUBROUTINE nemo_mapping
+
+   FUNCTION agrif_external_switch_index(ptx,pty,i1,isens)
+
+   USE dom_oce
+
+   INTEGER :: ptx, pty, i1, isens
+   INTEGER :: agrif_external_switch_index
+
+   IF( isens == 1 ) THEN
+      IF( ptx == 2 ) THEN ! T, V points
+         agrif_external_switch_index = jpiglo-i1+2
+      ELSE ! U, F points
+         agrif_external_switch_index = jpiglo-i1+1      
+      ENDIF
+   ELSE IF( isens ==2 ) THEN
+      IF ( pty == 2 ) THEN ! T, U points
+         agrif_external_switch_index = jpjglo-2-(i1 -jpjglo)
+      ELSE ! V, F points
+         agrif_external_switch_index = jpjglo-3-(i1 -jpjglo)
+      ENDIF
+   ENDIF
+
+   END FUNCTION agrif_external_switch_index
+
+   SUBROUTINE Correct_field(tab2d,i1,i2,j1,j2)
+      !!----------------------------------------------------------------------
+      !!                   *** ROUTINE Correct_field ***
+      !!----------------------------------------------------------------------
+   
+   USE dom_oce
+   USE agrif_oce
+
+   INTEGER :: i1,i2,j1,j2
+   REAL(wp), DIMENSION(i1:i2,j1:j2) :: tab2d
+
+   INTEGER :: i,j
+   REAL(wp), DIMENSION(i1:i2,j1:j2) :: tab2dtemp
+
+   tab2dtemp = tab2d
+
+   IF( .NOT. use_sign_north ) THEN
+      DO j=j1,j2
+         DO i=i1,i2
+            tab2d(i,j)=tab2dtemp(i2-(i-i1),j2-(j-j1))
+         END DO
+      END DO
+   ELSE
+      DO j=j1,j2
+         DO i=i1,i2
+            tab2d(i,j)=sign_north * tab2dtemp(i2-(i-i1),j2-(j-j1))
+         END DO
+      END DO
+   ENDIF
+
+   END SUBROUTINE Correct_field
+
 #else
    SUBROUTINE Subcalledbyagrif
Index: NEMO/trunk/src/OCE/DOM/dom_oce.F90
===================================================================
--- NEMO/trunk/src/OCE/DOM/dom_oce.F90	(revision 13215)
+++ NEMO/trunk/src/OCE/DOM/dom_oce.F90	(revision 13216)
@@ -215,4 +215,5 @@
 #if defined key_agrif
    LOGICAL, PUBLIC, PARAMETER ::   lk_agrif = .TRUE.    !: agrif flag
+   LOGICAL, PUBLIC            ::   lk_south, lk_north, lk_west, lk_east !: Child grid boundaries (interpolation or not)
 #else
    LOGICAL, PUBLIC, PARAMETER ::   lk_agrif = .FALSE.   !: agrif flag
Index: NEMO/trunk/src/OCE/DOM/domain.F90
===================================================================
--- NEMO/trunk/src/OCE/DOM/domain.F90	(revision 13215)
+++ NEMO/trunk/src/OCE/DOM/domain.F90	(revision 13216)
@@ -187,6 +187,11 @@
       ENDIF
       !
+
       IF( lk_c1d         )   CALL cor_c1d       ! 1D configuration: Coriolis set at T-point
       !
+
+#if defined key_agrif
+      IF( .NOT. Agrif_Root() ) CALL Agrif_Init_Domain( Kbb, Kmm, Kaa )
+#endif
       IF( ln_meshmask    )   CALL dom_wri       ! Create a domain file
       IF( .NOT.ln_rstart )   CALL dom_ctl       ! Domain control
@@ -307,4 +312,11 @@
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namrun in configuration namelist' )
       IF(lwm) WRITE ( numond, namrun )
+
+#if defined key_agrif
+      IF( .NOT. Agrif_Root() ) THEN
+            nn_it000 = (Agrif_Parent(nn_it000)-1)*Agrif_IRhot() + 1
+            nn_itend =  Agrif_Parent(nn_itend)   *Agrif_IRhot()
+      ENDIF
+#endif
       !
       IF(lwp) THEN                  ! control print
@@ -403,4 +415,10 @@
       IF(lwm) WRITE( numond, namdom )
       !
+#if defined key_agrif
+      IF( .NOT. Agrif_Root() ) THEN
+            rn_Dt = Agrif_Parent(rn_Dt) / Agrif_Rhot()
+      ENDIF
+#endif
+      !
       IF(lwp) THEN
          WRITE(numout,*)
Index: NEMO/trunk/src/OCE/DOM/istate.F90
===================================================================
--- NEMO/trunk/src/OCE/DOM/istate.F90	(revision 13215)
+++ NEMO/trunk/src/OCE/DOM/istate.F90	(revision 13216)
@@ -34,4 +34,8 @@
    USE lib_mpp         ! MPP library
    USE restart         ! restart
+#if defined key_agrif
+   USE agrif_oce_interp
+   USE agrif_oce
+#endif   
 
    IMPLICIT NONE
@@ -69,5 +73,5 @@
 !!gm  Why not include in the first call of dta_tsd ?  
 !!gm  probably associated with the use of internal damping...
-                     CALL dta_tsd_init        ! Initialisation of T & S input data
+       CALL dta_tsd_init        ! Initialisation of T & S input data
 !!gm to be moved in usrdef of C1D case
 !      IF( lk_c1d )   CALL dta_uvd_init        ! Initialization of U & V input data
@@ -83,4 +87,17 @@
 #endif
 
+#if defined key_agrif
+      IF ( (.NOT.Agrif_root()).AND.ln_init_chfrpar ) THEN
+         numror = 0                           ! define numror = 0 -> no restart file to read
+         ln_1st_euler = .true.                ! Set time-step indicator at nit000 (euler forward)
+         CALL day_init 
+         CALL agrif_istate( Kbb, Kmm, Kaa )   ! Interp from parent
+         !
+         ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) 
+         ssh (:,:,Kmm)     = ssh(:,:,Kbb)
+         uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
+         vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
+      ELSE
+#endif
       IF( ln_rstart ) THEN                    ! Restart from a file
          !                                    ! -------------------
@@ -99,4 +116,7 @@
             !
             ssh(:,:,Kbb)   = 0._wp               ! set the ocean at rest
+            uu  (:,:,:,Kbb) = 0._wp
+            vv  (:,:,:,Kbb) = 0._wp  
+            !
             IF( ll_wd ) THEN
                ssh(:,:,Kbb) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD 
@@ -110,7 +130,5 @@
                END_2D
             ENDIF 
-            uu  (:,:,:,Kbb) = 0._wp
-            vv  (:,:,:,Kbb) = 0._wp  
-            !
+             !
          ELSE                                 ! user defined initial T and S
             CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )         
@@ -147,4 +165,7 @@
          ! 
       ENDIF 
+#if defined key_agrif
+      ENDIF
+#endif
       ! 
       ! Initialize "now" and "before" barotropic velocities:
Index: NEMO/trunk/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/dynspg_ts.F90	(revision 13215)
+++ NEMO/trunk/src/OCE/DYN/dynspg_ts.F90	(revision 13216)
@@ -513,4 +513,10 @@
          END_2D
          !
+         ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
+         IF( ln_bdy )   CALL bdy_ssh( ssha_e )
+#if defined key_agrif
+         IF( .NOT.Agrif_Root() )   CALL agrif_ssh_ts( jn )
+#endif
+
          CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp,  zhU, 'U', -1._wp,  zhV, 'V', -1._wp )
          !
@@ -525,9 +531,4 @@
          END IF
          !
-         ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
-         IF( ln_bdy )   CALL bdy_ssh( ssha_e )
-#if defined key_agrif
-         IF( .NOT.Agrif_Root() )   CALL agrif_ssh_ts( jn )
-#endif
          !  
          ! Sea Surface Height at u-,v-points (vvl case only)
@@ -643,9 +644,17 @@
          ENDIF
        
-         IF( .NOT.ln_linssh ) THEN   !* Update ocean depth (variable volume case only)
+         IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
             hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
             hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) )
             hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)
             hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) )
+         ENDIF
+         !                                                 ! open boundaries
+         IF( ln_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
+#if defined key_agrif                                                           
+         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
+#endif
+         !
+         IF( .NOT.ln_linssh ) THEN   !* Update ocean depth (variable volume case only)
             CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  &
                  &                         , hu_e , 'U',  1._wp, hv_e , 'V',  1._wp  &
@@ -655,10 +664,4 @@
          ENDIF
          !
-         !
-         !                                                 ! open boundaries
-         IF( ln_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
-#if defined key_agrif                                                           
-         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
-#endif
          !                                             !* Swap
          !                                             !  ----
Index: NEMO/trunk/src/OCE/DYN/sshwzv.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/sshwzv.F90	(revision 13215)
+++ NEMO/trunk/src/OCE/DYN/sshwzv.F90	(revision 13216)
@@ -204,31 +204,40 @@
          ! Mask vertical velocity at first/last columns/row 
          ! inside computational domain (cosmetic) 
-         ! --- West --- !
-         DO ji = mi0(2), mi1(2)
-            DO jj = 1, jpj
-               pww(ji,jj,:) = 0._wp 
+         ! --- West --- !         
+         IF( lk_west) THEN
+            DO ji = mi0(2), mi1(2)
+               DO jj = 1, jpj
+                  pww(ji,jj,:) = 0._wp 
+               ENDDO
             ENDDO
-         ENDDO
+         ENDIF
          !
          ! --- East --- !
-         DO ji = mi0(jpiglo-1), mi1(jpiglo-1)
-            DO jj = 1, jpj
-               pww(ji,jj,:) = 0._wp
+         IF( lk_east) THEN
+            DO ji = mi0(jpiglo-1), mi1(jpiglo-1)
+               DO jj = 1, jpj
+                  pww(ji,jj,:) = 0._wp
+               ENDDO
             ENDDO
-         ENDDO
+         ENDIF
          !
          ! --- South --- !
-         DO jj = mj0(2), mj1(2)
-            DO ji = 1, jpi
-               pww(ji,jj,:) = 0._wp
+         IF( lk_south) THEN
+            DO jj = mj0(2), mj1(2)
+               DO ji = 1, jpi
+                  pww(ji,jj,:) = 0._wp
+               ENDDO
             ENDDO
-         ENDDO
+         ENDIF
          !
          ! --- North --- !
-         DO jj = mj0(jpjglo-1), mj1(jpjglo-1)
-            DO ji = 1, jpi
-               pww(ji,jj,:) = 0._wp
+         IF( lk_north) THEN
+            DO jj = mj0(jpjglo-1), mj1(jpjglo-1)
+               DO ji = 1, jpi
+                  pww(ji,jj,:) = 0._wp
+               ENDDO
             ENDDO
-         ENDDO
+         ENDIF
+         !
       ENDIF 
 #endif 
Index: NEMO/trunk/src/OCE/FLO/floblk.F90
===================================================================
--- NEMO/trunk/src/OCE/FLO/floblk.F90	(revision 13215)
+++ NEMO/trunk/src/OCE/FLO/floblk.F90	(revision 13216)
@@ -41,4 +41,8 @@
       INTEGER, INTENT( in  ) ::   Kbb, Kmm ! ocean time level indices
       !!
+#ifndef key_agrif
+
+!RB super quick fix to compile with agrif
+
       INTEGER :: jfl              ! dummy loop arguments
       INTEGER :: ind, ifin, iloop
@@ -364,4 +368,5 @@
          GO TO 222
       ENDIF
+#endif
       !
       !
Index: NEMO/trunk/src/OCE/LBC/lib_mpp.F90
===================================================================
--- NEMO/trunk/src/OCE/LBC/lib_mpp.F90	(revision 13215)
+++ NEMO/trunk/src/OCE/LBC/lib_mpp.F90	(revision 13216)
@@ -137,7 +137,7 @@
 
    ! Communications summary report
-   CHARACTER(len=128), DIMENSION(:), ALLOCATABLE ::   crname_lbc                   !: names of lbc_lnk calling routines
-   CHARACTER(len=128), DIMENSION(:), ALLOCATABLE ::   crname_glb                   !: names of global comm calling routines
-   CHARACTER(len=128), DIMENSION(:), ALLOCATABLE ::   crname_dlg                   !: names of delayed global comm calling routines
+   CHARACTER(len=lca), DIMENSION(:), ALLOCATABLE ::   crname_lbc                   !: names of lbc_lnk calling routines
+   CHARACTER(len=lca), DIMENSION(:), ALLOCATABLE ::   crname_glb                   !: names of global comm calling routines
+   CHARACTER(len=lca), DIMENSION(:), ALLOCATABLE ::   crname_dlg                   !: names of delayed global comm calling routines
    INTEGER, PUBLIC                               ::   ncom_stp = 0                 !: copy of time step # istp
    INTEGER, PUBLIC                               ::   ncom_fsbc = 1                !: copy of sbc time step # nn_fsbc
Index: NEMO/trunk/src/OCE/LBC/mppini.F90
===================================================================
--- NEMO/trunk/src/OCE/LBC/mppini.F90	(revision 13215)
+++ NEMO/trunk/src/OCE/LBC/mppini.F90	(revision 13216)
@@ -102,4 +102,9 @@
             &           'the domain is lay out for distributed memory computing!' )
          !
+#if defined key_agrif
+    IF (.NOT.agrif_root()) THEN
+      call agrif_nemo_init()
+    ENDIF
+#endif
    END SUBROUTINE mpp_init
 
@@ -333,9 +338,5 @@
 #if defined key_agrif
       IF( .NOT. Agrif_Root() ) THEN       ! AGRIF children: specific setting (cf. agrif_user.F90)
-         IF( jpiglo /= nbcellsx + 2 + 2*nbghostcells )   &
-            CALL ctl_stop( 'STOP', 'mpp_init: Agrif children requires jpiglo == nbcellsx + 2 + 2*nbghostcells' )
-         IF( jpjglo /= nbcellsy + 2 + 2*nbghostcells )   &
-            CALL ctl_stop( 'STOP', 'mpp_init: Agrif children requires jpjglo == nbcellsy + 2 + 2*nbghostcells' )
-         IF( ln_use_jattr )   CALL ctl_stop( 'STOP', 'mpp_init:Agrif children requires ln_use_jattr = .false. ' )
+         CALL agrif_nemo_init()
       ENDIF
 #endif
Index: NEMO/trunk/src/OCE/OBS/diaobs.F90
===================================================================
--- NEMO/trunk/src/OCE/OBS/diaobs.F90	(revision 13215)
+++ NEMO/trunk/src/OCE/OBS/diaobs.F90	(revision 13216)
@@ -94,5 +94,5 @@
    TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) ::   profdataqc   !: Profile data after quality control
 
-   CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE ::   cobstypesprof, cobstypessurf   !: Profile & surface obs types
+   CHARACTER(len=lca), PUBLIC, DIMENSION(:), ALLOCATABLE ::   cobstypesprof, cobstypessurf   !: Profile & surface obs types
 
    !!----------------------------------------------------------------------
Index: NEMO/trunk/src/OCE/SBC/sbcmod.F90
===================================================================
--- NEMO/trunk/src/OCE/SBC/sbcmod.F90	(revision 13215)
+++ NEMO/trunk/src/OCE/SBC/sbcmod.F90	(revision 13216)
@@ -121,4 +121,5 @@
 #endif
       !
+      !
       IF(lwp) THEN                  !* Control print
          WRITE(numout,*) '   Namelist namsbc (partly overwritten with CPP key setting)'
Index: NEMO/trunk/src/OCE/STO/stopar.F90
===================================================================
--- NEMO/trunk/src/OCE/STO/stopar.F90	(revision 13215)
+++ NEMO/trunk/src/OCE/STO/stopar.F90	(revision 13216)
@@ -56,6 +56,6 @@
    INTEGER,          DIMENSION(:),       ALLOCATABLE :: sto3d_ord  ! order of autoregressive process
 
-   CHARACTER(len=1), DIMENSION(:),       ALLOCATABLE :: sto2d_typ  ! nature of grid point (T, U, V, W, F, I)
-   CHARACTER(len=1), DIMENSION(:),       ALLOCATABLE :: sto3d_typ  ! nature of grid point (T, U, V, W, F, I)
+   CHARACTER(len=lca), DIMENSION(:),       ALLOCATABLE :: sto2d_typ  ! nature of grid point (T, U, V, W, F, I)
+   CHARACTER(len=lca), DIMENSION(:),       ALLOCATABLE :: sto3d_typ  ! nature of grid point (T, U, V, W, F, I)
    REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto2d_sgn  ! control of the sign accross the north fold
    REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto3d_sgn  ! control of the sign accross the north fold
Index: NEMO/trunk/src/OCE/USR/usrdef_hgr.F90
===================================================================
--- NEMO/trunk/src/OCE/USR/usrdef_hgr.F90	(revision 13215)
+++ NEMO/trunk/src/OCE/USR/usrdef_hgr.F90	(revision 13216)
@@ -95,4 +95,5 @@
 #if defined key_agrif
       ! ! Upper left longitude and latitude from parent:
+      ! Laurent: Should be modify in case of an east-west cyclic parent grid
       IF (.NOT.Agrif_root()) THEN
          zlam0 = zlam1 + Agrif_irhox() * REAL(Agrif_Parent(jpjglo)-2 , wp) * ze1deg * zcos_alpha  &
Index: NEMO/trunk/src/OCE/USR/usrdef_nam.F90
===================================================================
--- NEMO/trunk/src/OCE/USR/usrdef_nam.F90	(revision 13215)
+++ NEMO/trunk/src/OCE/USR/usrdef_nam.F90	(revision 13216)
@@ -74,6 +74,6 @@
 #if defined key_agrif
       IF( .NOT. Agrif_Root() ) THEN
-         kpi  = nbcellsx + 2 + 2*nbghostcells
-         kpj  = nbcellsy + 2 + 2*nbghostcells
+         kpi  = nbcellsx + 2 + 2*nbghostcells_x
+         kpj  = nbcellsy + 2 + 2*nbghostcells_y_s
       ENDIF
 #endif
Index: NEMO/trunk/src/OCE/nemogcm.F90
===================================================================
--- NEMO/trunk/src/OCE/nemogcm.F90	(revision 13215)
+++ NEMO/trunk/src/OCE/nemogcm.F90	(revision 13216)
@@ -143,11 +143,7 @@
 #if defined key_agrif
       Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
-      CALL Agrif_Declare_Var_dom   ! AGRIF: set the meshes for DOM
-      CALL Agrif_Declare_Var       !  "      "   "   "      "  DYN/TRA 
+      CALL Agrif_Declare_Var       !  "      "   "   "      "  DYN/TRA
 # if defined key_top
       CALL Agrif_Declare_Var_top   !  "      "   "   "      "  TOP
-# endif
-# if defined key_si3
-      CALL Agrif_Declare_Var_ice   !  "      "   "   "      "  Sea ice
 # endif
 #endif
@@ -400,5 +396,7 @@
       ! Initialise time level indices
       Nbb = 1; Nnn = 2; Naa = 3; Nrhs = Naa
-
+#if defined key_agrif
+      Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
+#endif 
       !                             !-------------------------------!
       !                             !  NEMO general initialization  !
@@ -415,5 +413,12 @@
       IF( lk_c1d       )   CALL     c1d_init        ! 1D column configuration
                            CALL     wad_init        ! Wetting and drying options
+
+#if defined key_agrif
+     CALL Agrif_Declare_Var_ini   !  "      "   "   "      "  DOM
+#endif
                            CALL     dom_init( Nbb, Nnn, Naa, "OPA") ! Domain
+
+
+
       IF( ln_crs       )   CALL     crs_init(      Nnn )       ! coarsened grid: domain initialization 
       IF( sn_cfctl%l_prtctl )   &
@@ -436,5 +441,5 @@
       ENDIF
       !
-      
+
                            CALL  istate_init( Nbb, Nnn, Naa )    ! ocean initial state (Dynamics and tracers)
 
Index: NEMO/trunk/src/OCE/par_kind.F90
===================================================================
--- NEMO/trunk/src/OCE/par_kind.F90	(revision 13215)
+++ NEMO/trunk/src/OCE/par_kind.F90	(revision 13216)
@@ -31,5 +31,6 @@
    
    !                                                                !!** Integer **
-   INTEGER, PUBLIC, PARAMETER ::   lc = 256                          !: Lenght of Character strings
+   INTEGER, PUBLIC, PARAMETER ::   lc  = 256                          !: Lenght of Character strings
+   INTEGER, PUBLIC, PARAMETER ::   lca = 400                          !: Lenght of Character arrays
 
    !!----------------------------------------------------------------------
Index: NEMO/trunk/src/OCE/par_oce.F90
===================================================================
--- NEMO/trunk/src/OCE/par_oce.F90	(revision 13215)
+++ NEMO/trunk/src/OCE/par_oce.F90	(revision 13216)
@@ -47,7 +47,10 @@
    ! global domain size for AGRIF     !!! * total AGRIF computational domain *
    INTEGER, PUBLIC            ::   nbug_in_agrif_conv_do_not_remove_or_modify = 1 - 1
-   INTEGER, PUBLIC, PARAMETER ::   nbghostcells = 3                             !: number of ghost cells
-   INTEGER, PUBLIC            ::   nbcellsx   ! = jpiglo - 2 - 2*nbghostcells   !: number of cells in i-direction
-   INTEGER, PUBLIC            ::   nbcellsy   ! = jpjglo - 2 - 2*nbghostcells   !: number of cells in j-direction
+   INTEGER, PUBLIC, PARAMETER ::   nbghostcells = 3 !: number of ghost cells: default value
+   INTEGER, PUBLIC            ::   nbghostcells_x   !: number of ghost cells in i-direction
+   INTEGER, PUBLIC            ::   nbghostcells_y_s   !: number of ghost cells in j-direction at south
+   INTEGER, PUBLIC            ::   nbghostcells_y_n   !: number of ghost cells in j-direction at north                       !: number of ghost cells
+   INTEGER, PUBLIC            ::   nbcellsx   ! = jpiglo - 2 - 2*nbghostcells_x   !: number of cells in i-direction
+   INTEGER, PUBLIC            ::   nbcellsy   ! = jpjglo - 2 - 2*nbghostcells-y   !: number of cells in j-direction
 
    ! local domain size                !!! * local computational domain *
Index: NEMO/trunk/src/OCE/stpctl.F90
===================================================================
--- NEMO/trunk/src/OCE/stpctl.F90	(revision 13215)
+++ NEMO/trunk/src/OCE/stpctl.F90	(revision 13216)
@@ -119,37 +119,27 @@
       !                                   !==            test of local extrema           ==!
       !                                   !==  done by all processes at every time step  ==!
-      !
-      ! define zmax default value. needed for land processors
-      IF( ll_colruns ) THEN    ! default value: must not be kept when calling mpp_max -> must be as small as possible
-         zmax(:) = -HUGE(1._wp)
-      ELSE                     ! default value: must not give true for any of the tests bellow (-> avoid manipulating HUGE...)
-         zmax(:) =  0._wp
-         zmax(3) = -1._wp      ! avoid salinity minimum at 0.
-      ENDIF
-      !
       llmsk(:,:,1) = ssmask(:,:) == 1._wp
-      IF( COUNT( llmsk(:,:,1) ) > 0 ) THEN   ! avoid huge values sent back for land processors...
-         IF( ll_wd ) THEN
-            zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) )   ! ssh max
+      IF( ll_wd ) THEN
+         zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) )   ! ssh max
+      ELSE
+         zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm)           ), mask = llmsk(:,:,1) )   ! ssh max
+      ENDIF
+      llmsk(:,:,:) = umask(:,:,:) == 1._wp
+      zmax(2) = MAXVAL(  ABS( uu(:,:,:,Kmm) ), mask = llmsk )                     ! velocity max (zonal only)
+      llmsk(:,:,:) = tmask(:,:,:) == 1._wp
+      zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     ! minus salinity max
+      zmax(4) = MAXVAL(  ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     !       salinity max
+      IF( ll_colruns .OR. jpnij == 1 ) THEN     ! following variables are used only in the netcdf file
+         zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  ! minus temperature max
+         zmax(6) = MAXVAL(  ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  !       temperature max
+         IF( ln_zad_Aimp ) THEN
+            zmax(7) = MAXVAL(   Cu_adv(:,:,:)   , mask = llmsk )                  ! partitioning coeff. max
+            llmsk(:,:,:) = wmask(:,:,:) == 1._wp
+            zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = llmsk )                  ! implicit vertical vel. max
          ELSE
-            zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm)           ), mask = llmsk(:,:,1) )   ! ssh max
-         ENDIF
-      ENDIF
-      zmax(2) = MAXVAL( ABS( uu(:,:,:,Kmm) ) )                                       ! velocity max (zonal only)
-      llmsk(:,:,:) = tmask(:,:,:) == 1._wp
-      IF( COUNT( llmsk(:,:,:) ) > 0 ) THEN   ! avoid huge values sent back for land processors...
-         zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     ! minus salinity max
-         zmax(4) = MAXVAL(  ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     !       salinity max
-         IF( ll_colruns .OR. jpnij == 1 ) THEN     ! following variables are used only in the netcdf file
-            zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  ! minus temperature max
-            zmax(6) = MAXVAL(  ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  !       temperature max
-            IF( ln_zad_Aimp ) THEN
-               zmax(7) = MAXVAL(   Cu_adv(:,:,:)   , mask = llmsk )                  ! partitioning coeff. max
-               llmsk(:,:,:) = wmask(:,:,:) == 1._wp
-               IF( COUNT( llmsk(:,:,:) ) > 0 ) THEN   ! avoid huge values sent back for land processors...
-                  zmax(8) = MAXVAL(ABS( wi(:,:,:) ), mask = llmsk )                  ! implicit vertical vel. max
-               ENDIF
-            ENDIF
-         ENDIF
+            zmax(7:8) = 0._wp
+         ENDIF
+      ELSE
+         zmax(5:8) = 0._wp
       ENDIF
       zmax(9) = REAL( nstop, wp )                                              ! stop indicator
Index: NEMO/trunk/src/SAS/nemogcm.F90
===================================================================
--- NEMO/trunk/src/SAS/nemogcm.F90	(revision 13215)
+++ NEMO/trunk/src/SAS/nemogcm.F90	(revision 13216)
@@ -91,11 +91,7 @@
 #if defined key_agrif
       Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
-      CALL Agrif_Declare_Var_dom   ! AGRIF: set the meshes for DOM
       CALL Agrif_Declare_Var       !  "      "   "   "      "  DYN/TRA 
 # if defined key_top
       CALL Agrif_Declare_Var_top   !  "      "   "   "      "  TOP
-# endif
-# if defined key_si3
-      CALL Agrif_Declare_Var_ice   !  "      "   "   "      "  Sea ice
 # endif
 #endif
@@ -345,4 +341,7 @@
       ! Initialise time level indices
       Nbb = 1; Nnn = 2; Naa = 3; Nrhs = Naa
+#if defined key_agrif
+      Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
+#endif 
 
       !                             !-------------------------------!
@@ -358,4 +357,7 @@
                            CALL phy_cst         ! Physical constants
                            CALL eos_init        ! Equation of seawater
+#if defined key_agrif
+     CALL Agrif_Declare_Var_ini   !  "      "   "   "      "  DOM
+#endif
                            CALL dom_init( Nbb, Nnn, Naa, 'SAS') ! Domain
       IF( sn_cfctl%l_prtctl )   &
Index: NEMO/trunk/tests/CANAL/MY_SRC/stpctl.F90
===================================================================
--- NEMO/trunk/tests/CANAL/MY_SRC/stpctl.F90	(revision 13215)
+++ NEMO/trunk/tests/CANAL/MY_SRC/stpctl.F90	(revision 13216)
@@ -19,12 +19,12 @@
    USE dom_oce         ! ocean space and time domain variables 
    USE c1d             ! 1D vertical configuration
+   USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables
+   USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy
+   !  
    USE diawri          ! Standard run outputs       (dia_wri_state routine)
-   !
    USE in_out_manager  ! I/O manager
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp         ! distributed memory computing
-   USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables
-   USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy
-
+   !
    USE netcdf          ! NetCDF library
    IMPLICIT NONE
@@ -33,6 +33,6 @@
    PUBLIC stp_ctl           ! routine called by step.F90
 
-   INTEGER  ::   idrun, idtime, idssh, idu, ids1, ids2, idt1, idt2, idc1, idw1, istatus
-   LOGICAL  ::   lsomeoce
+   INTEGER                ::   nrunid   ! netcdf file id
+   INTEGER, DIMENSION(8)  ::   nvarid   ! netcdf variable id
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -42,5 +42,5 @@
 CONTAINS
 
-   SUBROUTINE stp_ctl( kt, Kbb, Kmm, kindic )
+   SUBROUTINE stp_ctl( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE stp_ctl  ***
@@ -50,5 +50,5 @@
       !! ** Method  : - Save the time step in numstp
       !!              - Print it each 50 time steps
-      !!              - Stop the run IF problem encountered by setting indic=-3
+      !!              - Stop the run IF problem encountered by setting nstop > 0
       !!                Problems checked: |ssh| maximum larger than 10 m
       !!                                  |U|   maximum larger than 10 m/s 
@@ -57,143 +57,241 @@
       !! ** Actions :   "time.step" file = last ocean time-step
       !!                "run.stat"  file = run statistics
-      !!                nstop indicator sheared among all local domain (lk_mpp=T)
+      !!                 nstop indicator sheared among all local domain
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
-      INTEGER, INTENT(in   ) ::   Kbb, Kmm      ! ocean time level index
-      INTEGER, INTENT(inout) ::   kindic   ! error indicator
-      !!
-      INTEGER                ::   ji, jj, jk          ! dummy loop indices
-      INTEGER, DIMENSION(2)  ::   ih                  ! min/max loc indices
-      INTEGER, DIMENSION(3)  ::   iu, is1, is2        ! min/max loc indices
-      REAL(wp)               ::   zzz                 ! local real 
-      REAL(wp), DIMENSION(9) ::   zmax
-      LOGICAL                ::   ll_wrtstp, ll_colruns, ll_wrtruns
-      CHARACTER(len=20) :: clname
-      !!----------------------------------------------------------------------
-      !
-      ll_wrtstp  = ( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend )
-      ll_colruns = ll_wrtstp .AND. ( sn_cfctl%l_runstat )
-      ll_wrtruns = ll_colruns .AND. lwm
-      IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'stp_ctl : time-stepping control'
-         WRITE(numout,*) '~~~~~~~'
-         !                                ! open time.step file
-         IF( lwm ) CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
-         !                                ! open run.stat file(s) at start whatever
-         !                                ! the value of sn_cfctl%ptimincr
-         IF( lwm .AND. ( sn_cfctl%l_runstat ) ) THEN
+      INTEGER, INTENT(in   ) ::   Kmm      ! ocean time level index
+      !!
+      INTEGER                         ::   ji                                    ! dummy loop indices
+      INTEGER                         ::   idtime, istatus
+      INTEGER , DIMENSION(9)          ::   iareasum, iareamin, iareamax
+      INTEGER , DIMENSION(3,4)        ::   iloc                                  ! min/max loc indices
+      REAL(wp)                        ::   zzz                                   ! local real 
+      REAL(wp), DIMENSION(9)          ::   zmax, zmaxlocal
+      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns
+      LOGICAL, DIMENSION(jpi,jpj,jpk) ::   llmsk
+      CHARACTER(len=20)               ::   clname
+      !!----------------------------------------------------------------------
+      IF( nstop > 0 .AND. ngrdstop > -1 )   RETURN   !   stpctl was already called by a child grid
+      !
+      ll_wrtstp  = ( MOD( kt-nit000, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend )
+      ll_colruns = ll_wrtstp .AND. sn_cfctl%l_runstat .AND. jpnij > 1 
+      ll_wrtruns = ( ll_colruns .OR. jpnij == 1 ) .AND. lwm
+      !
+      IF( kt == nit000 ) THEN
+         !
+         IF( lwp ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'stp_ctl : time-stepping control'
+            WRITE(numout,*) '~~~~~~~'
+         ENDIF
+         !                                ! open time.step    ascii file, done only by 1st subdomain
+         IF( lwm )   CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+         !
+         IF( ll_wrtruns ) THEN
+            !                             ! open run.stat     ascii file, done only by 1st subdomain
             CALL ctl_opn( numrun, 'run.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+            !                             ! open run.stat.nc netcdf file, done only by 1st subdomain
             clname = 'run.stat.nc'
             IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
-            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, idrun )
-            istatus = NF90_DEF_DIM( idrun, 'time', NF90_UNLIMITED, idtime )
-            istatus = NF90_DEF_VAR( idrun, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), idssh )
-            istatus = NF90_DEF_VAR( idrun,   'abs_u_max', NF90_DOUBLE, (/ idtime /), idu   )
-            istatus = NF90_DEF_VAR( idrun,       's_min', NF90_DOUBLE, (/ idtime /), ids1  )
-            istatus = NF90_DEF_VAR( idrun,       's_max', NF90_DOUBLE, (/ idtime /), ids2  )
-            istatus = NF90_DEF_VAR( idrun,       't_min', NF90_DOUBLE, (/ idtime /), idt1  )
-            istatus = NF90_DEF_VAR( idrun,       't_max', NF90_DOUBLE, (/ idtime /), idt2  )
+            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, nrunid )
+            istatus = NF90_DEF_DIM( nrunid, 'time', NF90_UNLIMITED, idtime )
+            istatus = NF90_DEF_VAR( nrunid, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), nvarid(1) )
+            istatus = NF90_DEF_VAR( nrunid,   'abs_u_max', NF90_DOUBLE, (/ idtime /), nvarid(2) )
+            istatus = NF90_DEF_VAR( nrunid,       's_min', NF90_DOUBLE, (/ idtime /), nvarid(3) )
+            istatus = NF90_DEF_VAR( nrunid,       's_max', NF90_DOUBLE, (/ idtime /), nvarid(4) )
+            istatus = NF90_DEF_VAR( nrunid,       't_min', NF90_DOUBLE, (/ idtime /), nvarid(5) )
+            istatus = NF90_DEF_VAR( nrunid,       't_max', NF90_DOUBLE, (/ idtime /), nvarid(6) )
             IF( ln_zad_Aimp ) THEN
-               istatus = NF90_DEF_VAR( idrun,   'abs_wi_max', NF90_DOUBLE, (/ idtime /), idw1  )
-               istatus = NF90_DEF_VAR( idrun,       'Cf_max', NF90_DOUBLE, (/ idtime /), idc1  )
+               istatus = NF90_DEF_VAR( nrunid,   'Cf_max', NF90_DOUBLE, (/ idtime /), nvarid(7) )
+               istatus = NF90_DEF_VAR( nrunid,'abs_wi_max',NF90_DOUBLE, (/ idtime /), nvarid(8) )
             ENDIF
-            istatus = NF90_ENDDEF(idrun)
-            zmax(8:9) = 0._wp    ! initialise to zero in case ln_zad_Aimp option is not in use
-         ENDIF
-      ENDIF
-      IF( kt == nit000 )   lsomeoce = COUNT( ssmask(:,:) == 1._wp ) > 0
-      !
-      IF(lwm .AND. ll_wrtstp) THEN        !==  current time step  ==!   ("time.step" file)
+            istatus = NF90_ENDDEF(nrunid)
+         ENDIF
+         !    
+      ENDIF
+      !
+      !                                   !==              write current time step              ==!
+      !                                   !==  done only by 1st subdomain at writting timestep  ==!
+      IF( lwm .AND. ll_wrtstp ) THEN
          WRITE ( numstp, '(1x, i8)' )   kt
          REWIND( numstp )
       ENDIF
-      !
-      !                                   !==  test of extrema  ==!
+      !                                   !==            test of local extrema           ==!
+      !                                   !==  done by all processes at every time step  ==!
+      llmsk(:,:,1) = ssmask(:,:) == 1._wp
       IF( ll_wd ) THEN
-         zmax(1) = MAXVAL(  ABS( ssh(:,:,Kmm) + ssh_ref*tmask(:,:,1) )  )        ! ssh max 
+         zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) )   ! ssh max
       ELSE
-         zmax(1) = MAXVAL(  ABS( ssh(:,:,Kmm) )  )                               ! ssh max
-      ENDIF
-      zmax(2) = MAXVAL(  ABS( uu(:,:,:,Kmm) )  )                                  ! velocity max (zonal only)
-      zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp )   ! minus salinity max
-      zmax(4) = MAXVAL(  ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp )   !       salinity max
-      zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp )   ! minus temperature max
-      zmax(6) = MAXVAL(  ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp )   !       temperature max
-      zmax(7) = REAL( nstop , wp )                                            ! stop indicator
-      IF( ln_zad_Aimp ) THEN
-         zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = wmask(:,:,:) == 1._wp ) ! implicit vertical vel. max
-         zmax(9) = MAXVAL(   Cu_adv(:,:,:)   , mask = tmask(:,:,:) == 1._wp ) ! partitioning coeff. max
-      ENDIF
-      !
+         zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm)           ), mask = llmsk(:,:,1) )   ! ssh max
+      ENDIF
+      llmsk(:,:,:) = umask(:,:,:) == 1._wp
+      zmax(2) = MAXVAL(  ABS( uu(:,:,:,Kmm) ), mask = llmsk )                     ! velocity max (zonal only)
+      llmsk(:,:,:) = tmask(:,:,:) == 1._wp
+      zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     ! minus salinity max
+      zmax(4) = MAXVAL(  ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     !       salinity max
+      IF( ll_colruns .OR. jpnij == 1 ) THEN     ! following variables are used only in the netcdf file
+         zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  ! minus temperature max
+         zmax(6) = MAXVAL(  ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  !       temperature max
+         IF( ln_zad_Aimp ) THEN
+            zmax(7) = MAXVAL(   Cu_adv(:,:,:)   , mask = llmsk )                  ! partitioning coeff. max
+            llmsk(:,:,:) = wmask(:,:,:) == 1._wp
+            zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = llmsk )                  ! implicit vertical vel. max
+         ELSE
+            zmax(7:8) = 0._wp
+         ENDIF
+      ELSE
+         zmax(5:8) = 0._wp
+      ENDIF
+      zmax(9) = REAL( nstop, wp )                                              ! stop indicator
+      !                                   !==               get global extrema             ==!
+      !                                   !==  done by all processes if writting run.stat  ==!
       IF( ll_colruns ) THEN
+         zmaxlocal(:) = zmax(:)
          CALL mpp_max( "stpctl", zmax )          ! max over the global domain
-         nstop = NINT( zmax(7) )                 ! nstop indicator sheared among all local domains
-      ENDIF
-      !                                   !==  run statistics  ==!   ("run.stat" files)
+         nstop = NINT( zmax(9) )                 ! update nstop indicator (now sheared among all local domains)
+      ENDIF
+      !                                   !==              write "run.stat" files              ==!
+      !                                   !==  done only by 1st subdomain at writting timestep  ==!
       IF( ll_wrtruns ) THEN
          WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3), zmax(4)
-         istatus = NF90_PUT_VAR( idrun, idssh, (/ zmax(1)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( idrun,   idu, (/ zmax(2)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( idrun,  ids1, (/-zmax(3)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( idrun,  ids2, (/ zmax(4)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( idrun,  idt1, (/-zmax(5)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( idrun,  idt2, (/ zmax(6)/), (/kt/), (/1/) )
+         istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) )
+         istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) )
+         istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/-zmax(3)/), (/kt/), (/1/) )
+         istatus = NF90_PUT_VAR( nrunid, nvarid(4), (/ zmax(4)/), (/kt/), (/1/) )
+         istatus = NF90_PUT_VAR( nrunid, nvarid(5), (/-zmax(5)/), (/kt/), (/1/) )
+         istatus = NF90_PUT_VAR( nrunid, nvarid(6), (/ zmax(6)/), (/kt/), (/1/) )
          IF( ln_zad_Aimp ) THEN
-            istatus = NF90_PUT_VAR( idrun,  idw1, (/ zmax(8)/), (/kt/), (/1/) )
-            istatus = NF90_PUT_VAR( idrun,  idc1, (/ zmax(9)/), (/kt/), (/1/) )
-         ENDIF
-         IF( MOD( kt , 100 ) == 0 ) istatus = NF90_SYNC(idrun)
-         IF( kt == nitend         ) istatus = NF90_CLOSE(idrun)
+            istatus = NF90_PUT_VAR( nrunid, nvarid(7), (/ zmax(7)/), (/kt/), (/1/) )
+            istatus = NF90_PUT_VAR( nrunid, nvarid(8), (/ zmax(8)/), (/kt/), (/1/) )
+         ENDIF
+         IF( kt == nitend )   istatus = NF90_CLOSE(nrunid)
       END IF
-      !                                   !==  error handling  ==!
-      IF( ( sn_cfctl%l_glochk .OR. lsomeoce ) .AND. (   &  ! domain contains some ocean points, check for sensible ranges
-         &  zmax(1) >   20._wp .OR.   &                    ! too large sea surface height ( > 20 m )
-         &  zmax(2) >   10._wp .OR.   &                    ! too large velocity ( > 10 m/s)
-!!$         &  zmax(3) >=   0._wp .OR.   &                    ! negative or zero sea surface salinity
-!!$         &  zmax(4) >= 100._wp .OR.   &                    ! too large sea surface salinity ( > 100 )
-!!$         &  zmax(4) <    0._wp .OR.   &                    ! too large sea surface salinity (keep this line for sea-ice)
-         &  ISNAN( zmax(1) + zmax(2) + zmax(3) ) ) ) THEN   ! NaN encounter in the tests
-         IF( lk_mpp .AND. sn_cfctl%l_glochk ) THEN
-            ! have use mpp_max (because sn_cfctl%l_glochk=.T. and distributed)
-            CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,Kmm))        , ssmask(:,:)  , zzz, ih  )
-            CALL mpp_maxloc( 'stpctl', ABS(uu(:,:,:,Kmm))          , umask (:,:,:), zzz, iu  )
-            CALL mpp_minloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is1 )
-            CALL mpp_maxloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is2 )
+      !                                   !==               error handling               ==!
+      !                                   !==  done by all processes at every time step  ==!
+      !
+      IF(   zmax(1) >   20._wp .OR.   &                   ! too large sea surface height ( > 20 m )
+         &  zmax(2) >   10._wp .OR.   &                   ! too large velocity ( > 10 m/s)
+!!$         &  zmax(3) >=   0._wp .OR.   &                   ! negative or zero sea surface salinity
+!!$         &  zmax(4) >= 100._wp .OR.   &                   ! too large sea surface salinity ( > 100 )
+!!$         &  zmax(4) <    0._wp .OR.   &                   ! too large sea surface salinity (keep this line for sea-ice)
+         &  ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR.   &               ! NaN encounter in the tests
+         &  ABS(   zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN   ! Infinity encounter in the tests
+         !
+         iloc(:,:) = 0
+         IF( ll_colruns ) THEN   ! zmax is global, so it is the same on all subdomains -> no dead lock with mpp_maxloc
+            ! first: close the netcdf file, so we can read it
+            IF( lwm .AND. kt /= nitend )   istatus = NF90_CLOSE(nrunid)
+            ! get global loc on the min/max
+            CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,         Kmm)), ssmask(:,:  ), zzz, iloc(1:2,1) )   ! mpp_maxloc ok if mask = F 
+            CALL mpp_maxloc( 'stpctl', ABS( uu(:,:,:,       Kmm)),  umask(:,:,:), zzz, iloc(1:3,2) )
+            CALL mpp_minloc( 'stpctl',      ts(:,:,:,jp_sal,Kmm) ,  tmask(:,:,:), zzz, iloc(1:3,3) )
+            CALL mpp_maxloc( 'stpctl',      ts(:,:,:,jp_sal,Kmm) ,  tmask(:,:,:), zzz, iloc(1:3,4) )
+            ! find which subdomain has the max.
+            iareamin(:) = jpnij+1   ;   iareamax(:) = 0   ;   iareasum(:) = 0
+            DO ji = 1, 9
+               IF( zmaxlocal(ji) == zmax(ji) ) THEN
+                  iareamin(ji) = narea   ;   iareamax(ji) = narea   ;   iareasum(ji) = 1
+               ENDIF
+            END DO
+            CALL mpp_min( "stpctl", iareamin )         ! min over the global domain
+            CALL mpp_max( "stpctl", iareamax )         ! max over the global domain
+            CALL mpp_sum( "stpctl", iareasum )         ! sum over the global domain
+         ELSE                    ! find local min and max locations:
+            ! if we are here, this means that the subdomain contains some oce points -> no need to test the mask used in maxloc
+            iloc(1:2,1) = MAXLOC( ABS( ssh(:,:,         Kmm)), mask = ssmask(:,:  ) == 1._wp ) + (/ nimpp - 1, njmpp - 1    /)
+            iloc(1:3,2) = MAXLOC( ABS(  uu(:,:,:,       Kmm)), mask =  umask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)
+            iloc(1:3,3) = MINLOC(       ts(:,:,:,jp_sal,Kmm) , mask =  tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)
+            iloc(1:3,4) = MAXLOC(       ts(:,:,:,jp_sal,Kmm) , mask =  tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)
+            iareamin(:) = narea   ;   iareamax(:) = narea   ;   iareasum(:) = 1         ! this is local information
+         ENDIF
+         !
+         WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests'
+         CALL wrt_line( ctmp2, kt, '|ssh| max',  zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) )
+         CALL wrt_line( ctmp3, kt, '|U|   max',  zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) )
+         CALL wrt_line( ctmp4, kt, 'Sal   min', -zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) )
+         CALL wrt_line( ctmp5, kt, 'Sal   max',  zmax(4), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) )
+         IF( Agrif_Root() ) THEN
+            WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort* files'
          ELSE
-            ! find local min and max locations
-            ih(:)  = MAXLOC( ABS( ssh(:,:,Kmm)   )                              ) + (/ nimpp - 1, njmpp - 1    /)
-            iu(:)  = MAXLOC( ABS( uu  (:,:,:,Kmm) )                              ) + (/ nimpp - 1, njmpp - 1, 0 /)
-            is1(:) = MINLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)
-            is2(:) = MAXLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)
-         ENDIF
-         
-         WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests'
-         WRITE(ctmp2,9100) kt,   zmax(1), ih(1) , ih(2)
-         WRITE(ctmp3,9200) kt,   zmax(2), iu(1) , iu(2) , iu(3)
-         WRITE(ctmp4,9300) kt, - zmax(3), is1(1), is1(2), is1(3)
-         WRITE(ctmp5,9400) kt,   zmax(4), is2(1), is2(2), is2(3)
-         WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort.nc file'
-         
+            WRITE(ctmp6,*) '      ===> output of last computed fields in '//TRIM(Agrif_CFixed())//'_output.abort* files'
+         ENDIF
+         !
          CALL dia_wri_state( Kmm, 'output.abort' )     ! create an output.abort file
-         
-         IF( .NOT. sn_cfctl%l_glochk ) THEN
-            WRITE(ctmp8,*) 'E R R O R message from sub-domain: ', narea
-            CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp8, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ctmp6 )
-         ELSE
-            CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6, ' ' )
-         ENDIF
-
-         kindic = -3
-         !
-      ENDIF
-      !
-9100  FORMAT (' kt=',i8,'   |ssh| max: ',1pg11.4,', at  i j  : ',2i5)
-9200  FORMAT (' kt=',i8,'   |U|   max: ',1pg11.4,', at  i j k: ',3i5)
-9300  FORMAT (' kt=',i8,'   S     min: ',1pg11.4,', at  i j k: ',3i5)
-9400  FORMAT (' kt=',i8,'   S     max: ',1pg11.4,', at  i j k: ',3i5)
+         !
+         IF( ll_colruns .or. jpnij == 1 ) THEN   ! all processes synchronized -> use lwp to print in opened ocean.output files
+            IF(lwp) THEN   ;   CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 )
+            ELSE           ;   nstop = MAX(1, nstop)   ! make sure nstop > 0 (automatically done when calling ctl_stop)
+            ENDIF
+         ELSE                                    ! only mpi subdomains with errors are here -> STOP now
+            CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 )
+         ENDIF
+         !
+      ENDIF
+      !
+      IF( nstop > 0 ) THEN                                                  ! an error was detected and we did not abort yet...
+         ngrdstop = Agrif_Fixed()                                           ! store which grid got this error
+         IF( .NOT. ll_colruns .AND. jpnij > 1 )   CALL ctl_stop( 'STOP' )   ! we must abort here to avoid MPI deadlock
+      ENDIF
+      !
 9500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16)
       !
    END SUBROUTINE stp_ctl
+
+
+   SUBROUTINE wrt_line( cdline, kt, cdprefix, pval, kloc, ksum, kmin, kmax )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE wrt_line  ***
+      !!
+      !! ** Purpose :   write information line
+      !!
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*),      INTENT(  out) ::   cdline
+      CHARACTER(len=*),      INTENT(in   ) ::   cdprefix
+      REAL(wp),              INTENT(in   ) ::   pval
+      INTEGER, DIMENSION(3), INTENT(in   ) ::   kloc
+      INTEGER,               INTENT(in   ) ::   kt, ksum, kmin, kmax
+      !
+      CHARACTER(len=80) ::   clsuff
+      CHARACTER(len=9 ) ::   clkt, clsum, clmin, clmax
+      CHARACTER(len=9 ) ::   cli, clj, clk
+      CHARACTER(len=1 ) ::   clfmt
+      CHARACTER(len=4 ) ::   cl4   ! needed to be able to compile with Agrif, I don't know why
+      INTEGER           ::   ifmtk
+      !!----------------------------------------------------------------------
+      WRITE(clkt , '(i9)') kt
+      
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpnij  ,wp))) + 1     ! how many digits to we need to write ? (we decide max = 9)
+      !!! WRITE(clsum, '(i'//clfmt//')') ksum                   ! this is creating a compilation error with AGRIF
+      cl4 = '(i'//clfmt//')'   ;   WRITE(clsum, cl4) ksum
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1    ! how many digits to we need to write ? (we decide max = 9)
+      cl4 = '(i'//clfmt//')'   ;   WRITE(clmin, cl4) kmin-1
+                                   WRITE(clmax, cl4) kmax-1
+      !
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpiglo,wp))) + 1      ! how many digits to we need to write jpiglo? (we decide max = 9)
+      cl4 = '(i'//clfmt//')'   ;   WRITE(cli, cl4) kloc(1)      ! this is ok with AGRIF
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpjglo,wp))) + 1      ! how many digits to we need to write jpjglo? (we decide max = 9)
+      cl4 = '(i'//clfmt//')'   ;   WRITE(clj, cl4) kloc(2)      ! this is ok with AGRIF
+      !
+      IF( ksum == 1 ) THEN   ;   WRITE(clsuff,9100) TRIM(clmin)
+      ELSE                   ;   WRITE(clsuff,9200) TRIM(clsum), TRIM(clmin), TRIM(clmax)
+      ENDIF
+      IF(kloc(3) == 0) THEN
+         ifmtk = INT(LOG10(REAL(jpk,wp))) + 1                   ! how many digits to we need to write jpk? (we decide max = 9)
+         clk = REPEAT(' ', ifmtk)                               ! create the equivalent in blank string
+         WRITE(cdline,9300) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj), clk(1:ifmtk), TRIM(clsuff)
+      ELSE
+         WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpk,wp))) + 1      ! how many digits to we need to write jpk? (we decide max = 9)
+         !!! WRITE(clk, '(i'//clfmt//')') kloc(3)               ! this is creating a compilation error with AGRIF
+         cl4 = '(i'//clfmt//')'   ;   WRITE(clk, cl4) kloc(3)   ! this is ok with AGRIF
+         WRITE(cdline,9400) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj),    TRIM(clk), TRIM(clsuff)
+      ENDIF
+      !
+9100  FORMAT('MPI rank ', a)
+9200  FORMAT('found in ', a, ' MPI tasks, spread out among ranks ', a, ' to ', a)
+9300  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j   ', a, ' ', a, ' ', a, ' ', a)
+9400  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j k ', a, ' ', a, ' ', a, ' ', a)
+      !
+   END SUBROUTINE wrt_line
+
 
    !!======================================================================
Index: NEMO/trunk/tests/VORTEX/EXPREF/1_namelist_cfg
===================================================================
--- NEMO/trunk/tests/VORTEX/EXPREF/1_namelist_cfg	(revision 13215)
+++ NEMO/trunk/tests/VORTEX/EXPREF/1_namelist_cfg	(revision 13216)
@@ -98,8 +98,6 @@
 &namagrif      !  AGRIF zoom                                            ("key_agrif")
 !-----------------------------------------------------------------------
-   ln_spc_dyn    = .true.  !  use 0 as special value for dynamics
-   rn_sponge_tra =  800.   !  coefficient for tracer   sponge layer [m2/s]
-   rn_sponge_dyn =  800.   !  coefficient for dynamics sponge layer [m2/s]
-   ln_chk_bathy  = .FALSE. !
+   rn_sponge_tra =  0.00768   !  coefficient for tracer   sponge layer []
+   rn_sponge_dyn =  0.00768   !  coefficient for dynamics sponge layer []
 /
 !!======================================================================
