diff --git a/.gitlab-ci-mercator.yml b/.gitlab-ci-mercator.yml
index ed644ce4e4354b378efd8c66cfbcec48e2120008..73135bce09ef8ff7106553dee258c6dc2d1fd2e8 100644
--- a/.gitlab-ci-mercator.yml
+++ b/.gitlab-ci-mercator.yml
@@ -3,8 +3,8 @@
   parallel:
     matrix:
       - HPC: [aa,belenos]
-        #CONFIG: [ORCA2_ICE_PISCES,ORCA2_OFF_PISCES,AMM12,AGRIF_DEMO,WED025,GYRE_PISCES,ORCA2_SAS_ICE,ORCA2_ICE_OBS,C1D_PAPA]
-        CONFIG: [ORCA2_ICE_PISCES,ORCA2_OFF_PISCES,AMM12,WED025,GYRE_PISCES,ORCA2_SAS_ICE,ORCA2_ICE_OBS,C1D_PAPA]
+        CONFIG: [ORCA2_ICE_PISCES,ORCA2_OFF_PISCES,AMM12,AGRIF_DEMO,WED025,GYRE_PISCES,ORCA2_SAS_ICE,ORCA2_ICE_OBS,C1D_PAPA]
+        #CONFIG: [ORCA2_ICE_PISCES,ORCA2_OFF_PISCES,AMM12,WED025,GYRE_PISCES,ORCA2_SAS_ICE,ORCA2_ICE_OBS,C1D_PAPA]
 
 # HPC & testcases lists to be tested with SETTE
 .parallel_HPC_TST:
diff --git a/cfgs/AGRIF_DEMO/cpp_AGRIF_DEMO.fcm b/cfgs/AGRIF_DEMO/cpp_AGRIF_DEMO.fcm
index 0e63a0c88c5a15ff518eba03d436c3b64111a310..238173ed2dd7c244a4f50ea46f30d2b6fc9decb4 100644
--- a/cfgs/AGRIF_DEMO/cpp_AGRIF_DEMO.fcm
+++ b/cfgs/AGRIF_DEMO/cpp_AGRIF_DEMO.fcm
@@ -1 +1 @@
-bld::tool::fppkeys   key_si3 key_top key_xios key_agrif key_qco key_vco_1d3d
+bld::tool::fppkeys   key_si3 key_top key_xios key_agrif key_qco  key_vco_3d
diff --git a/src/NST/agrif_oce_interp.F90 b/src/NST/agrif_oce_interp.F90
index 19b62117bdade68832dcd509f8c4aa3811a9a347..1675ffee3a347bfed6c70e1f6cea04dd5c0735b7 100644
--- a/src/NST/agrif_oce_interp.F90
+++ b/src/NST/agrif_oce_interp.F90
@@ -1720,7 +1720,7 @@ CONTAINS
       !!----------------------------------------------------------------------  
       !    
       IF( before ) THEN
-         IF ( lk_vco_1d3d ) THEN
+         IF ( l_zps ) THEN
             DO jk = k1, k2
                DO jj = j1, j2
                   DO ji = i1, i2
diff --git a/src/OCE/BDY/bdyini.F90 b/src/OCE/BDY/bdyini.F90
index 12fe7e728d4fb4a91b95153228337405a7beabaa..989ca8162a8189b3d049ab88b0c4d21e967330e7 100644
--- a/src/OCE/BDY/bdyini.F90
+++ b/src/OCE/BDY/bdyini.F90
@@ -1594,12 +1594,10 @@ CONTAINS
       DO ib = 1, nbdysegw
          ! get mask at boundary extremities:
          ztestmask(1:2)=0.
-         DO ji = 1, jpi
-            DO jj = 1, jpj             
-              IF( mig(ji,0) == jpiwob(ib) .AND. mjg(jj,0) == jpjwdt(ib) )   ztestmask(1) = tmask(ji,jj,1)
-              IF( mig(ji,0) == jpiwob(ib) .AND. mjg(jj,0) == jpjwft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
-            END DO
-         END DO
+         DO_2D (0, 0, 0, 0) 
+            IF( mig(ji,0) == jpiwob(ib) .AND. mjg(jj,0) == jpjwdt(ib) )   ztestmask(1) = tmask(ji,jj,1)
+            IF( mig(ji,0) == jpiwob(ib) .AND. mjg(jj,0) == jpjwft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
+         END_2D
          CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
 
          IF (ztestmask(1)==1) THEN 
@@ -1630,12 +1628,10 @@ CONTAINS
       DO ib = 1, nbdysege
          ! get mask at boundary extremities:
          ztestmask(1:2)=0.
-         DO ji = 1, jpi
-            DO jj = 1, jpj             
-              IF( mig(ji,0) == jpieob(ib)+1 .AND. mjg(jj,0) == jpjedt(ib) )   ztestmask(1) = tmask(ji,jj,1)
-              IF( mig(ji,0) == jpieob(ib)+1 .AND. mjg(jj,0) == jpjeft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
-            END DO
-         END DO
+         DO_2D (0, 0, 0, 0) 
+            IF( mig(ji,0) == jpieob(ib)+1 .AND. mjg(jj,0) == jpjedt(ib) )   ztestmask(1) = tmask(ji,jj,1)
+            IF( mig(ji,0) == jpieob(ib)+1 .AND. mjg(jj,0) == jpjeft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
+         END_2D 
          CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
 
          IF (ztestmask(1)==1) THEN
@@ -1666,12 +1662,10 @@ CONTAINS
       DO ib = 1, nbdysegs
          ! get mask at boundary extremities:
          ztestmask(1:2)=0.
-         DO ji = 1, jpi
-            DO jj = 1, jpj             
-              IF( mjg(jj,0) == jpjsob(ib) .AND. mig(ji,0) == jpisdt(ib) )   ztestmask(1) = tmask(ji,jj,1)
-              IF( mjg(jj,0) == jpjsob(ib) .AND. mig(ji,0) == jpisft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
-            END DO
-         END DO
+         DO_2D (0, 0, 0, 0) 
+            IF( mjg(jj,0) == jpjsob(ib) .AND. mig(ji,0) == jpisdt(ib) )   ztestmask(1) = tmask(ji,jj,1)
+            IF( mjg(jj,0) == jpjsob(ib) .AND. mig(ji,0) == jpisft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
+         END_2D 
          CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
 
          IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN
@@ -1688,12 +1682,10 @@ CONTAINS
       DO ib = 1, nbdysegn
          ! get mask at boundary extremities:
          ztestmask(1:2)=0.
-         DO ji = 1, jpi
-            DO jj = 1, jpj             
-               IF( mjg(jj,0) == jpjnob(ib)+1 .AND. mig(ji,0) == jpindt(ib) )   ztestmask(1) = tmask(ji,jj,1)
-               IF( mjg(jj,0) == jpjnob(ib)+1 .AND. mig(ji,0) == jpinft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
-            END DO
-         END DO
+         DO_2D (0, 0, 0, 0) 
+            IF( mjg(jj,0) == jpjnob(ib)+1 .AND. mig(ji,0) == jpindt(ib) )   ztestmask(1) = tmask(ji,jj,1)
+            IF( mjg(jj,0) == jpjnob(ib)+1 .AND. mig(ji,0) == jpinft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
+         END_2D
          CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
 
          IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN
diff --git a/src/OCE/DOM/domzgr.F90 b/src/OCE/DOM/domzgr.F90
index ca1139c20dfde85ae434874f9bd5b2cbe6289de5..7106646015e5d6a1048bd5783ae324d9516ee2d0 100644
--- a/src/OCE/DOM/domzgr.F90
+++ b/src/OCE/DOM/domzgr.F90
@@ -22,7 +22,6 @@ MODULE domzgr
 
    !!----------------------------------------------------------------------
    !!   dom_zgr       : read or set the ocean vertical coordinate system
-   !!   zgr_read      : read the vertical information in the domain configuration file
    !!   zgr_top_bot   : ocean top and bottom level for t-, u, and v-points with 1 as minimum value
    !!---------------------------------------------------------------------
    USE oce            ! ocean variables
@@ -85,8 +84,15 @@ CONTAINS
          WRITE(numout,*)
          WRITE(numout,*) 'dom_zgr : vertical coordinate'
          WRITE(numout,*) '~~~~~~~'
-         IF( ln_linssh ) WRITE(numout,*) '          linear free surface: the vertical mesh does not change in time'
+         IF( lk_linssh ) WRITE(numout,*) '          linear free surface: the vertical mesh does not change in time'
       ENDIF
+      !                                ! Control keys
+      IF(.NOT.( lk_linssh.OR.lk_qco ) .OR. ( lk_linssh.AND.lk_qco ))   &
+         &              CALL ctl_stop( 'STOP','domzgr: either key_linssh or key_qco MUST be set up in cpp_* file !' )
+      !
+      IF(.NOT.lk_vco_1d .AND. .NOT.lk_vco_1d3d .AND. .NOT.lk_vco_3d )   &
+         &              CALL ctl_stop( 'STOP','domzgr: either key_vco_1d, key_vco_1d3d or key_vco_3d MUST be set up in cpp_* file !' )
+      !
       !                             !==============================!
       IF( ln_read_cfg ) THEN        !==  read in domcfg.nc file  ==!
          !                          !==============================!
@@ -136,19 +142,30 @@ CONTAINS
          CALL iom_get( inum, jpdom_global, 'bottom_level' , z2d(:,:)   )   ! last wet T-points
          k_bot(:,:) = NINT( z2d(:,:) )
          !
+         IF( iom_varid( inum, 'mbku', ldstop = .FALSE. ) > 0 ) THEN
+            IF(lwp) WRITE(numout,*) '          mbku, mbkv & mbkf read in ', TRIM(cn_domcfg), ' file'
+            CALL iom_get( inum, jpdom_global, 'mbku', z2d, cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
+            mbku(:,:) = NINT( z2d(:,:) )
+            CALL iom_get( inum, jpdom_global, 'mbkv', z2d, cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
+            mbkv(:,:) = NINT( z2d(:,:) )
+            CALL iom_get( inum, jpdom_global, 'mbkf', z2d, cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
+            mbkf(:,:) = NINT( z2d(:,:) )
+            is_mbkuvf = 1
+         ELSE
+            is_mbkuvf = 0
+         ENDIF
+         !
          !                                !--------------------!
          IF(     lk_vco_1d ) THEN         !--  z-coordinate  --!   use only 1D arrays for all gdep and e3 fields
             !                             !--------------------!
-            l_zco = .TRUE.                     ! old logical   ==> to be removed
-            l_zps = .FALSE.
-            l_sco = .FALSE.
+            !
+            IF(.NOT. l_zco )   CALL ctl_stop( 'STOP','domzgr: only l_zco is compatible with key_vco_1d. Fix domcfg or cpp_* file !' )
             !
             !                             !-----------------------!
          ELSEIF( lk_vco_1d3d ) THEN       !--  z-partial cells  --!   use 3D t-level e3
             !                             !-----------------------!
-            l_zco = .TRUE.
-            l_zps = .FALSE.                    ! old logical   ==> to be removed
-            l_sco = .FALSE.
+            !
+            IF( l_sco )   CALL ctl_stop( 'STOP','domzgr: l_sco is NOT compatible with key_vco_1d3d. Fix domcfg or cpp_* file !' )
             !
             !                                   ! t-level: 3D reference   (include partial cell)
             CALL iom_get( inum, jpdom_global, 'e3t_0'  , e3t_3d, cd_type = 'T', psgn = 1._wp, kfill = jpfillcopy )
@@ -159,9 +176,6 @@ CONTAINS
             !                             !--------------------!               
          ELSEIF( lk_vco_3d ) THEN         !--  s-coordinate  --!   use 3D for all gdep and e3 fields
             !                             !--------------------!
-            l_zco = .FALSE.
-            l_zps = .FALSE.
-            l_sco = .TRUE.                     ! old logical   ==> to be removed
             !
             !                                   !* depth         : 3D reference   (without partial cell)
 !!st no gdep._0 in ORCA_domcfg.nc from sette
@@ -179,16 +193,25 @@ CONTAINS
             !
 !!st PATCH for ORCA see above no gdep._0 in ORCA_domcfg.nc
             IF(  iom_varid( inum, 'gdept_0' , ldstop = .FALSE. ) > 0  .AND.  &
-               & iom_varid( inum, 'gdepw_0' , ldstop = .FALSE. ) > 0    ) THEN
+                 & iom_varid( inum, 'gdepw_0' , ldstop = .FALSE. ) > 0    ) THEN
                CALL ctl_warn( 'zgr_read : old definition of depths and scale factors used ', & 
-                  &           '           depths at t- and w-points read in the domain configuration file')
+                    &           '           depths at t- and w-points read in the domain configuration file')
                CALL iom_get( inum, jpdom_global , 'gdept_0' , gdept_3d, kfill = jpfillcopy )
                CALL iom_get( inum, jpdom_global , 'gdepw_0' , gdepw_3d, kfill = jpfillcopy )
                !
             ELSE                                !- depths computed from e3. scale factors
                CALL e3_to_depth( e3t_3d   , e3w_3d   , gdept_3d   , gdepw_3d    )    ! 3D depths
             ENDIF
-!!st
+
+!!st            IF(l_zps) THEN
+!!st               DO jk=1, jpk
+!!st                   e3w_3d(:,:,jk) = e3w_1d(jk)
+!!st                  e3uw_3d(:,:,jk) = e3w_1d(jk)
+!!st                  e3vw_3d(:,:,jk) = e3w_1d(jk)
+!!st                  gdept_3d(:,:,jk) = gdept_1d(jk)
+!!st                  gdepw_3d(:,:,jk) = gdepw_1d(jk)
+!!st               END DO  
+!!st            ENDIF
             !                                   !* reference depth for negative bathy (wetting and drying only)
             IF( ll_wd )   CALL iom_get( inum,  'rn_wd_ref_depth' , ssh_ref   )
             !
@@ -224,50 +247,49 @@ CONTAINS
          IF(lwp) WRITE(numout,*) '          User defined vertical mesh (usr_def_zgr)'
          is_mbkuvf = 0
          !
-         IF( lk_linssh .OR. lk_qco ) THEN
-            !                                !--------------------!
-            IF( lk_vco_1d ) THEN             !--  z-coordinate  --!   use only 1D arrays for all gdep and e3 fields
-               !                             !--------------------!
-               !
-               
-               CALL usr_def_zgr( l_zco  , l_zps  , l_sco, ln_isfcav,   &
-                  &              k_top   , k_bot                      ,   &    ! 1st & last ocean level
-                  &              gdept_1d, gdepw_1d, e3t_1d, e3w_1d   )        ! 1D gridpoints depth
-               !
-               IF( l_sco )   CALL ctl_stop( 'STOP','domzgr: key_vco_1d and l_sco=T are incompatible. Fix usrdef_zgr !' )
-               IF( l_zps )   CALL ctl_stop( 'STOP','domzgr: key_vco_1d and l_zps=T are incompatible. Fix usrdef_zgr !' )
-               !
-               !                             !-----------------------!
-            ELSEIF( lk_vco_1d3d ) THEN       !--  z-partial cells  --!   use 3D t-level e3
-               !                             !-----------------------!
-               !
-               CALL usr_def_zgr( l_zco  , l_zps  , l_sco, ln_isfcav   ,   &
-                  &              k_top   , k_bot                      ,   &    ! 1st & last ocean level
-                  &              gdept_1d, gdepw_1d, e3t_1d, e3w_1d   ,   &    ! 1D gridpoints depth
-                  &              e3t_3d  , e3u_3d  , e3v_3d, e3f_3d       )    ! vertical scale factors
-               !
-               IF( l_sco )   CALL ctl_stop( 'STOP','domzgr: key_vco_1d3d and l_sco=T are incompatible. Fix usrdef_zgr !' )
-               !
-               ! make sure that periodicities are properly applied 
-               CALL lbc_lnk( 'dom_zgr', e3t_3d, 'T', 1._wp,   e3u_3d, 'U', 1._wp,  e3v_3d, 'V', 1._wp, e3f_3d, 'F', 1._wp,   &
-                    &                     kfillmode = jpfillcopy )   ! do not put 0 over closed boundaries
-               !
-               !                             !--------------------!               
-            ELSEIF( lk_vco_3d ) THEN         !--  s-coordinate  --!   use 3D for all gdep and e3 fields
-               !                             !--------------------!
-               !
-               CALL usr_def_zgr( l_zco  , l_zps  , l_sco, ln_isfcav   ,   &
-                  &              k_top   , k_bot                      ,   &    ! 1st & last ocean level
-                  &              gdept_1d, gdepw_1d, e3t_1d, e3w_1d   ,   &    ! 1D gridpoints depth
-                  &              e3t_3d  , e3u_3d  , e3v_3d, e3f_3d   ,   &    ! vertical scale factors
-                  &              gdept_3d, gdepw_3d                   ,   &    ! gridpoints depth 
-                  &              e3w_3d  , e3uw_3d , e3vw_3d               )   ! vertical scale factors
-               CALL lbc_lnk( 'dom_zgr', gdept_3d, 'T', 1._wp, gdepw_3d, 'W', 1._wp,                                            &
-                  &                       e3t_3d, 'T', 1._wp,   e3u_3d, 'U', 1._wp,  e3v_3d, 'V', 1._wp, e3f_3d, 'F', 1._wp,   &
-                  &                       e3w_3d, 'W', 1._wp,  e3uw_3d, 'U', 1._wp, e3vw_3d, 'V', 1._wp,                       &   
-                  &                     kfillmode = jpfillcopy )   ! do not put 0 over closed boundaries
-            ENDIF
+         !                                !--------------------!
+         IF( lk_vco_1d ) THEN             !--  z-coordinate  --!   use only 1D arrays for all gdep and e3 fields
+            !                             !--------------------!
+            !
+            
+            CALL usr_def_zgr( l_zco  , l_zps  , l_sco, ln_isfcav,   &
+                 &              k_top   , k_bot                      ,   &    ! 1st & last ocean level
+                 &              gdept_1d, gdepw_1d, e3t_1d, e3w_1d   )        ! 1D gridpoints depth
+            !
+            IF( l_sco )   CALL ctl_stop( 'STOP','domzgr: key_vco_1d and l_sco=T are incompatible. Fix usrdef_zgr !' )
+            IF( l_zps )   CALL ctl_stop( 'STOP','domzgr: key_vco_1d and l_zps=T are incompatible. Fix usrdef_zgr !' )
+            !
+            !                             !-----------------------!
+         ELSEIF( lk_vco_1d3d ) THEN       !--  z-partial cells  --!   use 3D t-level e3
+            !                             !-----------------------!
+            !
+            CALL usr_def_zgr( l_zco  , l_zps  , l_sco, ln_isfcav   ,   &
+                 &              k_top   , k_bot                      ,   &    ! 1st & last ocean level
+                 &              gdept_1d, gdepw_1d, e3t_1d, e3w_1d   ,   &    ! 1D gridpoints depth
+                 &              e3t_3d  , e3u_3d  , e3v_3d, e3f_3d       )    ! vertical scale factors
+            !
+            IF( l_sco )   CALL ctl_stop( 'STOP','domzgr: key_vco_1d3d and l_sco=T are incompatible. Fix usrdef_zgr !' )
+            !
+            ! make sure that periodicities are properly applied 
+            CALL lbc_lnk( 'dom_zgr', e3t_3d, 'T', 1._wp,   e3u_3d, 'U', 1._wp,  e3v_3d, 'V', 1._wp, e3f_3d, 'F', 1._wp,   &
+                 &                     kfillmode = jpfillcopy )   ! do not put 0 over closed boundaries
+            !
+            !                             !--------------------!               
+         ELSEIF( lk_vco_3d ) THEN         !--  s-coordinate  --!   use 3D for all gdep and e3 fields
+            !                             !--------------------!
+            !
+            CALL usr_def_zgr( l_zco  , l_zps  , l_sco, ln_isfcav   ,   &
+                 &              k_top   , k_bot                      ,   &    ! 1st & last ocean level
+                 &              gdept_1d, gdepw_1d, e3t_1d, e3w_1d   ,   &    ! 1D gridpoints depth
+                 &              e3t_3d  , e3u_3d  , e3v_3d, e3f_3d   ,   &    ! vertical scale factors
+                 &              gdept_3d, gdepw_3d                   ,   &    ! gridpoints depth 
+                 &              e3w_3d  , e3uw_3d , e3vw_3d               )   ! vertical scale factors
+            CALL lbc_lnk( 'dom_zgr', gdept_3d, 'T', 1._wp, gdepw_3d, 'W', 1._wp,                                            &
+                 &                       e3t_3d, 'T', 1._wp,   e3u_3d, 'U', 1._wp,  e3v_3d, 'V', 1._wp, e3f_3d, 'F', 1._wp,   &
+                 &                       e3w_3d, 'W', 1._wp,  e3uw_3d, 'U', 1._wp, e3vw_3d, 'V', 1._wp,                       &   
+                 &                     kfillmode = jpfillcopy )   ! do not put 0 over closed boundaries
          ENDIF
+         !
          ztopbot(:,:,1) = REAL(k_top, wp)
          ztopbot(:,:,2) = REAL(k_bot, wp)
          CALL lbc_lnk( 'dom_zgr', ztopbot, 'T', 1._wp, kfillmode = jpfillcopy )   ! do not put 0 over closed boundaries
@@ -377,135 +399,135 @@ CONTAINS
    END SUBROUTINE dom_zgr
 
 
-   SUBROUTINE zgr_read( ld_zco  , ld_zps  , ld_sco  , ld_isfcav,   &   ! type of vertical coordinate
-      &                 pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d  ,   &   ! 1D reference vertical coordinate
-      &                 pdept , pdepw ,                            &   ! 3D t & w-points depth
-      &                 pe3t  , pe3u  , pe3v   , pe3f ,            &   ! vertical scale factors
-      &                 pe3w  , pe3uw , pe3vw         ,            &   !     -      -      -
-      &                 k_top  , k_bot  ,                          &   ! top & bottom ocean level
-      &                 k_mbkuvf  , k_bot_u  , k_bot_v  , k_bot_f  )   ! U/V/F points bottom levels
-      !!---------------------------------------------------------------------
-      !!              ***  ROUTINE zgr_read  ***
-      !!
-      !! ** Purpose :   Read the vertical information in the domain configuration file
-      !!
-      !!----------------------------------------------------------------------
-      LOGICAL                   , INTENT(out) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags
-      LOGICAL                   , INTENT(out) ::   ld_isfcav                   ! under iceshelf cavity flag
-      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth       [m]
-      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d           ! 1D vertical scale factors [m]
-      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pdept, pdepw                ! grid-point depth          [m]
-      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors    [m]
-      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3w , pe3uw, pe3vw         !    -       -      -
-      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top , k_bot               ! first & last ocean level
-      INTEGER                   , INTENT(out) ::   k_mbkuvf                    ! ==1 if mbku, mbkv, mbkf are in file
-      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_bot_u , k_bot_v, k_bot_f  ! bottom levels at U/V/F points
-      !
-      INTEGER  ::   ji,jj,jk     ! dummy loop index
-      INTEGER  ::   inum, iatt
-      REAL(WP) ::   z_zco, z_zps, z_sco, z_cav
-      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D workspace
-      CHARACTER(len=7) ::   catt   ! 'zco', 'zps, 'sco' or 'UNKNOWN'
-      !!----------------------------------------------------------------------
-      !
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) '   zgr_read : read the vertical coordinates in ', TRIM( cn_domcfg ), ' file'
-         WRITE(numout,*) '   ~~~~~~~~'
-      ENDIF
-      !
-      CALL iom_open( cn_domcfg, inum )
-      !
-      !                          !* type of vertical coordinate
-      CALL iom_getatt( inum, 'VertCoord', catt )   ! returns 'UNKNOWN' if not found
-      ld_zco = catt == 'zco'          ! default = .false.
-      ld_zps = catt == 'zps'          ! default = .false.
-      ld_sco = catt == 'sco'          ! default = .false.
-      !                          !* ocean cavities under iceshelves
-      CALL iom_getatt( inum,    'IsfCav', iatt )   ! returns -999 if not found
-      ld_isfcav = iatt == 1           ! default = .false.
-      !
-      ! ------- keep compatibility with OLD VERSION... start -------
-      IF( catt == 'UNKNOWN' ) THEN
-         CALL iom_get( inum,    'ln_zco', z_zco )   ;   ld_zco = z_zco /= 0._wp
-         CALL iom_get( inum,    'ln_zps', z_zps )   ;   ld_zps = z_zps /= 0._wp
-         CALL iom_get( inum,    'ln_sco', z_sco )   ;   ld_sco = z_sco /= 0._wp
-      ENDIF
-      IF( iatt == -999 ) THEN
-         CALL iom_get( inum, 'ln_isfcav', z_cav )   ;   ld_isfcav = z_cav /= 0._wp
-      ENDIF
-      ! ------- keep compatibility with OLD VERSION... end -------
-      !
-      !                          !* ocean top and bottom level
-      CALL iom_get( inum, jpdom_global, 'top_level'    , z2d   )   ! 1st wet T-points (ISF)
-      k_top(:,:) = NINT( z2d(:,:) )
-      CALL iom_get( inum, jpdom_global, 'bottom_level' , z2d   )   ! last wet T-points
-      k_bot(:,:) = NINT( z2d(:,:) )
-      !
-      !                          !* vertical scale factors
-      CALL iom_get( inum, jpdom_unknown, 'e3t_1d'  , pe3t_1d  )                     ! 1D reference coordinate
-      CALL iom_get( inum, jpdom_unknown, 'e3w_1d'  , pe3w_1d  )
-      !
-      CALL iom_get( inum, jpdom_global, 'e3t_0'  , pe3t , cd_type = 'T', psgn = 1._wp, kfill = jpfillcopy )    ! 3D coordinate
-      CALL iom_get( inum, jpdom_global, 'e3u_0'  , pe3u , cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
-      CALL iom_get( inum, jpdom_global, 'e3v_0'  , pe3v , cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
-      CALL iom_get( inum, jpdom_global, 'e3f_0'  , pe3f , cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
-      CALL iom_get( inum, jpdom_global, 'e3w_0'  , pe3w , cd_type = 'W', psgn = 1._wp, kfill = jpfillcopy )
-      CALL iom_get( inum, jpdom_global, 'e3uw_0' , pe3uw, cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
-      CALL iom_get( inum, jpdom_global, 'e3vw_0' , pe3vw, cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
-      !
-      !                          !* depths
-      !                                   !- old depth definition (obsolescent feature)
-      IF(  iom_varid( inum, 'gdept_1d', ldstop = .FALSE. ) > 0  .AND.  &
-         & iom_varid( inum, 'gdepw_1d', ldstop = .FALSE. ) > 0  .AND.  &
-         & iom_varid( inum, 'gdept_0' , ldstop = .FALSE. ) > 0  .AND.  &
-         & iom_varid( inum, 'gdepw_0' , ldstop = .FALSE. ) > 0    ) THEN
-         CALL ctl_warn( 'zgr_read : old definition of depths and scale factors used ', & 
-            &           '           depths at t- and w-points read in the domain configuration file')
-         CALL iom_get( inum, jpdom_unknown, 'gdept_1d', pdept_1d )   
-         CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', pdepw_1d )
-         CALL iom_get( inum, jpdom_global , 'gdept_0' , pdept, kfill = jpfillcopy )
-         CALL iom_get( inum, jpdom_global , 'gdepw_0' , pdepw, kfill = jpfillcopy )
-         !
-      ELSE                                !- depths computed from e3. scale factors
-         CALL e3_to_depth( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d )    ! 1D reference depth
-         CALL e3_to_depth( pe3t   , pe3w   , pdept   , pdepw    )    ! 3D depths
-#if defined key_qco && key_isf
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpk )        ! vertical sum at partial cell xxxx other level  
-            IF( jk == k_top(ji,jj) ) THEN                               ! first ocean point : partial cell
-               pdept(ji,jj,jk) = pdepw(ji,jj,jk  ) + 0.5_wp * pe3w(ji,jj,jk)   ! = risfdep + 1/2 e3w_0(mikt)
-            ELSE                                                        !  other levels
-               pdept(ji,jj,jk) = pdept(ji,jj,jk-1) +          pe3w(ji,jj,jk) 
-            ENDIF
-         END_3D
-#endif
-         IF(lwp) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
-            WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
-            WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
-         ENDIF
-      ENDIF
-      !
-      IF( iom_varid( inum, 'mbku', ldstop = .FALSE. ) > 0 ) THEN
-         IF(lwp) WRITE(numout,*) '          mbku, mbkv & mbkf read in ', TRIM(cn_domcfg), ' file'
-         CALL iom_get( inum, jpdom_global, 'mbku', z2d, cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
-         k_bot_u(:,:) = NINT( z2d(:,:) )
-         CALL iom_get( inum, jpdom_global, 'mbkv', z2d, cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
-         k_bot_v(:,:) = NINT( z2d(:,:) )
-         CALL iom_get( inum, jpdom_global, 'mbkf', z2d, cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
-         k_bot_f(:,:) = NINT( z2d(:,:) )
-         k_mbkuvf = 1
-      ELSE
-         k_mbkuvf = 0
-      ENDIF
-      !
-      ! reference depth for negative bathy (wetting and drying only)
-      IF( ll_wd )  CALL iom_get( inum,  'rn_wd_ref_depth' , ssh_ref   )
-      !
-      CALL iom_close( inum )
-      !
-   END SUBROUTINE zgr_read
+!!st   SUBROUTINE zgr_read( ld_zco  , ld_zps  , ld_sco  , ld_isfcav,   &   ! type of vertical coordinate
+!!st      &                 pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d  ,   &   ! 1D reference vertical coordinate
+!!st      &                 pdept , pdepw ,                            &   ! 3D t & w-points depth
+!!st      &                 pe3t  , pe3u  , pe3v   , pe3f ,            &   ! vertical scale factors
+!!st      &                 pe3w  , pe3uw , pe3vw         ,            &   !     -      -      -
+!!st      &                 k_top  , k_bot  ,                          &   ! top & bottom ocean level
+!!st      &                 k_mbkuvf  , k_bot_u  , k_bot_v  , k_bot_f  )   ! U/V/F points bottom levels
+!!st      !!---------------------------------------------------------------------
+!!st      !!              ***  ROUTINE zgr_read  ***
+!!st      !!
+!!st      !! ** Purpose :   Read the vertical information in the domain configuration file
+!!st      !!
+!!st      !!----------------------------------------------------------------------
+!!st      LOGICAL                   , INTENT(out) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags
+!!st      LOGICAL                   , INTENT(out) ::   ld_isfcav                   ! under iceshelf cavity flag
+!!st      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth       [m]
+!!st      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d           ! 1D vertical scale factors [m]
+!!st      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pdept, pdepw                ! grid-point depth          [m]
+!!st      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors    [m]
+!!st      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3w , pe3uw, pe3vw         !    -       -      -
+!!st      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top , k_bot               ! first & last ocean level
+!!st      INTEGER                   , INTENT(out) ::   k_mbkuvf                    ! ==1 if mbku, mbkv, mbkf are in file
+!!st      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_bot_u , k_bot_v, k_bot_f  ! bottom levels at U/V/F points
+!!st      !
+!!st      INTEGER  ::   ji,jj,jk     ! dummy loop index
+!!st      INTEGER  ::   inum, iatt
+!!st      REAL(WP) ::   z_zco, z_zps, z_sco, z_cav
+!!st      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D workspace
+!!st      CHARACTER(len=7) ::   catt   ! 'zco', 'zps, 'sco' or 'UNKNOWN'
+!!st      !!----------------------------------------------------------------------
+!!st      !
+!!st      IF(lwp) THEN
+!!st         WRITE(numout,*)
+!!st         WRITE(numout,*) '   zgr_read : read the vertical coordinates in ', TRIM( cn_domcfg ), ' file'
+!!st         WRITE(numout,*) '   ~~~~~~~~'
+!!st      ENDIF
+!!st      !
+!!st      CALL iom_open( cn_domcfg, inum )
+!!st      !
+!!st      !                          !* type of vertical coordinate
+!!st      CALL iom_getatt( inum, 'VertCoord', catt )   ! returns 'UNKNOWN' if not found
+!!st      ld_zco = catt == 'zco'          ! default = .false.
+!!st      ld_zps = catt == 'zps'          ! default = .false.
+!!st      ld_sco = catt == 'sco'          ! default = .false.
+!!st      !                          !* ocean cavities under iceshelves
+!!st      CALL iom_getatt( inum,    'IsfCav', iatt )   ! returns -999 if not found
+!!st      ld_isfcav = iatt == 1           ! default = .false.
+!!st      !
+!!st      ! ------- keep compatibility with OLD VERSION... start -------
+!!st      IF( catt == 'UNKNOWN' ) THEN
+!!st         CALL iom_get( inum,    'ln_zco', z_zco )   ;   ld_zco = z_zco /= 0._wp
+!!st         CALL iom_get( inum,    'ln_zps', z_zps )   ;   ld_zps = z_zps /= 0._wp
+!!st         CALL iom_get( inum,    'ln_sco', z_sco )   ;   ld_sco = z_sco /= 0._wp
+!!st      ENDIF
+!!st      IF( iatt == -999 ) THEN
+!!st         CALL iom_get( inum, 'ln_isfcav', z_cav )   ;   ld_isfcav = z_cav /= 0._wp
+!!st      ENDIF
+!!st      ! ------- keep compatibility with OLD VERSION... end -------
+!!st      !
+!!st      !                          !* ocean top and bottom level
+!!st      CALL iom_get( inum, jpdom_global, 'top_level'    , z2d   )   ! 1st wet T-points (ISF)
+!!st      k_top(:,:) = NINT( z2d(:,:) )
+!!st      CALL iom_get( inum, jpdom_global, 'bottom_level' , z2d   )   ! last wet T-points
+!!st      k_bot(:,:) = NINT( z2d(:,:) )
+!!st      !
+!!st      !                          !* vertical scale factors
+!!st      CALL iom_get( inum, jpdom_unknown, 'e3t_1d'  , pe3t_1d  )                     ! 1D reference coordinate
+!!st      CALL iom_get( inum, jpdom_unknown, 'e3w_1d'  , pe3w_1d  )
+!!st      !
+!!st      CALL iom_get( inum, jpdom_global, 'e3t_0'  , pe3t , cd_type = 'T', psgn = 1._wp, kfill = jpfillcopy )    ! 3D coordinate
+!!st      CALL iom_get( inum, jpdom_global, 'e3u_0'  , pe3u , cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
+!!st      CALL iom_get( inum, jpdom_global, 'e3v_0'  , pe3v , cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
+!!st      CALL iom_get( inum, jpdom_global, 'e3f_0'  , pe3f , cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
+!!st      CALL iom_get( inum, jpdom_global, 'e3w_0'  , pe3w , cd_type = 'W', psgn = 1._wp, kfill = jpfillcopy )
+!!st      CALL iom_get( inum, jpdom_global, 'e3uw_0' , pe3uw, cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
+!!st      CALL iom_get( inum, jpdom_global, 'e3vw_0' , pe3vw, cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
+!!st      !
+!!st      !                          !* depths
+!!st      !                                   !- old depth definition (obsolescent feature)
+!!st      IF(  iom_varid( inum, 'gdept_1d', ldstop = .FALSE. ) > 0  .AND.  &
+!!st         & iom_varid( inum, 'gdepw_1d', ldstop = .FALSE. ) > 0  .AND.  &
+!!st         & iom_varid( inum, 'gdept_0' , ldstop = .FALSE. ) > 0  .AND.  &
+!!st         & iom_varid( inum, 'gdepw_0' , ldstop = .FALSE. ) > 0    ) THEN
+!!st         CALL ctl_warn( 'zgr_read : old definition of depths and scale factors used ', & 
+!!st            &           '           depths at t- and w-points read in the domain configuration file')
+!!st         CALL iom_get( inum, jpdom_unknown, 'gdept_1d', pdept_1d )   
+!!st         CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', pdepw_1d )
+!!st         CALL iom_get( inum, jpdom_global , 'gdept_0' , pdept, kfill = jpfillcopy )
+!!st         CALL iom_get( inum, jpdom_global , 'gdepw_0' , pdepw, kfill = jpfillcopy )
+!!st         !
+!!st      ELSE                                !- depths computed from e3. scale factors
+!!st         CALL e3_to_depth( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d )    ! 1D reference depth
+!!st         CALL e3_to_depth( pe3t   , pe3w   , pdept   , pdepw    )    ! 3D depths
+!!st#if defined key_qco && key_isf
+!!st         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpk )        ! vertical sum at partial cell xxxx other level  
+!!st            IF( jk == k_top(ji,jj) ) THEN                               ! first ocean point : partial cell
+!!st               pdept(ji,jj,jk) = pdepw(ji,jj,jk  ) + 0.5_wp * pe3w(ji,jj,jk)   ! = risfdep + 1/2 e3w_0(mikt)
+!!st            ELSE                                                        !  other levels
+!!st               pdept(ji,jj,jk) = pdept(ji,jj,jk-1) +          pe3w(ji,jj,jk) 
+!!st            ENDIF
+!!st         END_3D
+!!st#endif
+!!st         IF(lwp) THEN
+!!st            WRITE(numout,*)
+!!st            WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
+!!st            WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
+!!st            WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
+!!st         ENDIF
+!!st      ENDIF
+!!st      !
+!!st      IF( iom_varid( inum, 'mbku', ldstop = .FALSE. ) > 0 ) THEN
+!!st         IF(lwp) WRITE(numout,*) '          mbku, mbkv & mbkf read in ', TRIM(cn_domcfg), ' file'
+!!st         CALL iom_get( inum, jpdom_global, 'mbku', z2d, cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
+!!st         k_bot_u(:,:) = NINT( z2d(:,:) )
+!!st         CALL iom_get( inum, jpdom_global, 'mbkv', z2d, cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
+!!st         k_bot_v(:,:) = NINT( z2d(:,:) )
+!!st         CALL iom_get( inum, jpdom_global, 'mbkf', z2d, cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
+!!st         k_bot_f(:,:) = NINT( z2d(:,:) )
+!!st         k_mbkuvf = 1
+!!st      ELSE
+!!st         k_mbkuvf = 0
+!!st      ENDIF
+!!st      !
+!!st      ! reference depth for negative bathy (wetting and drying only)
+!!st      IF( ll_wd )  CALL iom_get( inum,  'rn_wd_ref_depth' , ssh_ref   )
+!!st      !
+!!st      CALL iom_close( inum )
+!!st      !
+!!st   END SUBROUTINE zgr_read
 
 
    SUBROUTINE zgr_top_bot( k_top, k_bot, k_mbkuvf )
diff --git a/src/OCE/SBC/sbcfwb.F90 b/src/OCE/SBC/sbcfwb.F90
index 902a2061a3f04c96f84ab8f63a9443d36bedaf5e..5d8eb980a4cd0b3262069df5142f2fe6a1d84024 100644
--- a/src/OCE/SBC/sbcfwb.F90
+++ b/src/OCE/SBC/sbcfwb.F90
@@ -113,6 +113,7 @@ CONTAINS
             nn_hvolg_mth   = Agrif_parent(nn_hvolg_mth)
          ENDIF
 #endif
+         IF ( nn_ice/=2 ) nn_fwb_voltype = 2 ! Enforce liquid volume conservation with no sea-ice 
 
          IF(lwp) THEN
             WRITE(numout,*)
@@ -198,11 +199,6 @@ CONTAINS
          ! isf cavities are excluded because it can feedback to the melting with generation of inhibition of plumes
          ! and in case of no melt, it can generate HSSW.
          !
-         IF( nn_ice == 0 ) THEN
-            snwice_mass_b(:,:) = 0.e0               ! no sea-ice model is being used : no snow+ice mass
-            snwice_fmass (:,:) = 0.e0
-         ENDIF
-         !
       ENDIF
 
 
@@ -279,10 +275,11 @@ CONTAINS
                   SELECT CASE (nn_fwb_voltype)
                   CASE( 1 )
                      agrif_tmp(1) = glob_sum( 'sbcfwb', e1e2t(A2D(0)) * tmask_agrif(A2D(0)) * ( ssh(A2D(0),Kmm) + snwice_mass_b(A2D(0)) * r1_rho0 ))
+                     CALL Agrif_step_child_adj(glob_sum_volume_ice_oce_agrif) 
                   CASE( 2 )
                      agrif_tmp(1) = glob_sum( 'sbcfwb', e1e2t(A2D(0)) * tmask_agrif(A2D(0)) * ssh(A2D(0),Kmm))
+                     CALL Agrif_step_child_adj(glob_sum_volume_oce_agrif) 
                   END SELECT
-                  CALL Agrif_step_child_adj(glob_sum_volume_agrif) ! Get value over child grids
                   CALL mpp_min('sbcfwb', agrif_tmp(:)) ! Required with // sisters to populate the value of each grid on each processor
                   a_fwb = SUM(agrif_tmp) * rho0 / area ! Sum over all grids
                   IF ( a_fwb_b == 999._wp ) a_fwb_b = a_fwb
@@ -358,10 +355,11 @@ CONTAINS
                SELECT CASE (nn_fwb_voltype)
                CASE( 1 )
                   agrif_tmp(1) = glob_sum( 'sbcfwb', e1e2t(A2D(0)) * tmask_agrif(A2D(0)) * ( ssh(A2D(0),Kmm) + snwice_mass_b(A2D(0)) * r1_rho0 ))
+                  CALL Agrif_step_child_adj(glob_sum_volume_ice_oce_agrif) ! Get value over child grids
                CASE( 2 )
                   agrif_tmp(1) = glob_sum( 'sbcfwb', e1e2t(A2D(0)) * tmask_agrif(A2D(0)) * ssh(A2D(0),Kmm) )
+                  CALL Agrif_step_child_adj(glob_sum_volume_oce_agrif)     ! Get value over child grids
                END SELECT
-               CALL Agrif_step_child_adj(glob_sum_volume_agrif) ! Get value over child grids
                CALL mpp_min('sbcfwb', agrif_tmp(:)) ! Required with // sisters to populate the value of each grid on each processor
                a_fwb = SUM(agrif_tmp) ! Sum over all grids
 #else          
@@ -481,9 +479,16 @@ CONTAINS
       CASE ( 4 )                             !==  global mean fwf set to zero (ISOMIP case) ==!
          !
          IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
-            !                                                  ! fwf global mean (excluding ocean to ice/snow exchanges) 
-            emp_corr = glob_sum( 'sbcfwb', e1e2t(A2D(0)) * ( emp(A2D(0)) - rnf(A2D(0)) - fwfisf_cav(A2D(0)) - fwfisf_par(A2D(0)) &
-               &                                                      - snwice_fmass(A2D(0)) ) ) / area
+            !                                                  ! fwf global mean (excluding ocean to ice/snow exchanges)
+            SELECT CASE (nn_fwb_voltype)
+            CASE( 1 )
+               emp_corr = glob_sum( 'sbcfwb', e1e2t(A2D(0)) * ( emp(A2D(0)) - rnf(A2D(0)) - fwfisf_cav(A2D(0)) - fwfisf_par(A2D(0)) &
+                  &                                                      - snwice_fmass(A2D(0)) ) ) / area
+            CASE( 2 )
+               emp_corr = glob_sum( 'sbcfwb', e1e2t(A2D(0)) * ( emp(A2D(0)) - rnf(A2D(0)) - fwfisf_cav(A2D(0)) - fwfisf_par(A2D(0)) &
+                  &                                                                             ) ) / area
+            END SELECT
+!!st check with clem
             ! clem: use y_fwfnow instead to improve performance? 
             !y_fwfnow(1) = local_sum( e1e2t(A2D(0)) * ( emp(A2D(0)) - rnf(A2D(0)) - fwfisf_cav(A2D(0)) - fwfisf_par(A2D(0)) &
             !   &                                                   - snwice_fmass(A2D(0)) ) )
@@ -558,26 +563,38 @@ CONTAINS
 
    END SUBROUTINE glob_sum_area_agrif   
 
-   SUBROUTINE glob_sum_volume_agrif()
+   SUBROUTINE glob_sum_volume_ice_oce_agrif()
       !!---------------------------------------------------------------------
-      !!           ***  Compute volume with embedded zooms  ***
+      !!     ***  Compute volume with embedded zooms (ice + liquid)  ***
       !!----------------------------------------------------------------------
       INTEGER :: igrid
       !
       IF (Agrif_root()) RETURN
 
       igrid = agrif_fixed() + 1
-      IF (( nn_ice==2 ).OR.(nn_fwb_voltype==1)) THEN
-         ! NB: we use "now" value for snwice_mass over child grids since it has not been updated yet at the time
-         !     this call is made (e.g. when starting a new step over the parent grid)
+      IF ( nn_ice==2 ) THEN
+         ! NB1: nn_ice is known on child grids at this stage
+         ! NB2: we use "now" value for snwice_mass over child grids since it has not been updated yet at the time
+         !      this call is made (e.g. when starting a new step over the parent grid)
          agrif_tmp(igrid) = glob_sum( 'sbcfwb', e1e2t(:,:) * tmask_agrif(:,:) * ( ssh(:,:,Kmm_a) + snwice_mass(:,:) * r1_rho0 ))
       ELSE
          agrif_tmp(igrid) = glob_sum( 'sbcfwb', e1e2t(:,:) * tmask_agrif(:,:) * ssh(:,:,Kmm_a) )
       ENDIF
 
-   END SUBROUTINE glob_sum_volume_agrif 
+   END SUBROUTINE glob_sum_volume_ice_oce_agrif 
 
+   SUBROUTINE glob_sum_volume_oce_agrif()
+      !!---------------------------------------------------------------------
+      !!       ***  Compute volume with embedded zooms (liquid only)  ***
+      !!---------------------------------------------------------------------
+      INTEGER :: igrid
+      !
+      IF (Agrif_root()) RETURN
+
+      igrid = agrif_fixed() + 1
+      agrif_tmp(igrid) = glob_sum( 'sbcfwb', e1e2t(:,:) * tmask_agrif(:,:) * ssh(:,:,Kmm_a) )
 
+   END SUBROUTINE glob_sum_volume_oce_agrif 
 
 #endif
 
diff --git a/src/OCE/nemogcm.F90 b/src/OCE/nemogcm.F90
index 6ac75ede19ca8a4fc37a59e8c4092a3864f48713..34b228312106779beea66f074fbad9e3dcac1551 100644
--- a/src/OCE/nemogcm.F90
+++ b/src/OCE/nemogcm.F90
@@ -64,15 +64,11 @@ MODULE nemogcm
 #if defined key_nemocice_decomp
    USE ice_domain_size, only: nx_global, ny_global
 #endif
-#if defined key_qco   ||   defined key_linssh
 # if defined key_RK3
    USE stprk3
 # else
    USE stpmlf         ! NEMO time-stepping               (stp_MLF   routine)
 # endif
-#else
-   USE step           ! NEMO time-stepping                 (stp     routine)
-#endif
    !
    USE lib_mpp        ! distributed memory computing
    USE mppini         ! shared/distributed memory setting (mpp_init routine)
@@ -168,15 +164,11 @@ CONTAINS
       DO WHILE( istp <= nitend .AND. nstop == 0 )
          !
          ncom_stp = istp
-#  if defined key_qco   ||   defined key_linssh
 #   if defined key_RK3
          CALL stp_RK3
 #   else
          CALL stp_MLF
 #   endif
-#  else
-         CALL stp
-#  endif
          istp = istp + 1
       END DO
       !
@@ -193,15 +185,11 @@ CONTAINS
                IF ( istp ==         nitend ) elapsed_time = zstptiming - elapsed_time
             ENDIF
             !
-#  if defined key_qco   ||   defined key_linssh
 #   if defined key_RK3
             CALL stp_RK3( istp )
 #   else
             CALL stp_MLF( istp )
 #   endif
-#  else
-            CALL stp    ( istp )
-#  endif
             istp = istp + 1
             !
             IF( lwp .AND. ln_timing )   WRITE(numtime,*) 'timing step ', istp-1, ' : ', MPI_Wtime() - zstptiming
diff --git a/src/OFF/nemogcm.F90 b/src/OFF/nemogcm.F90
index 6ea5e372ab2429be2d38a521af9b86c97a1a5505..9c011cfb7e758f387085cce0c7efce4b79add525 100644
--- a/src/OFF/nemogcm.F90
+++ b/src/OFF/nemogcm.F90
@@ -63,11 +63,7 @@ MODULE nemogcm
    USE prtctl         ! Print control                    (prt_ctl_init routine)
    USE timing         ! Timing
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
-#if defined key_qco || defined key_linssh
    USE stpmlf , ONLY : Nbb, Nnn, Naa, Nrhs   ! time level indices
-#else
-   USE step    , ONLY : Nbb, Nnn, Naa, Nrhs   ! time level indices
-#endif
    USE halo_mng
 
    IMPLICIT NONE
diff --git a/tests/DOME/EXPREF/1_namelist_cfg b/tests/DOME/EXPREF/1_namelist_cfg
index 7cafa15f9cff6de024b101ae5b15a6e9dcd43395..ebca85137d19906e3042779c238e5fabd51cb3be 100644
--- a/tests/DOME/EXPREF/1_namelist_cfg
+++ b/tests/DOME/EXPREF/1_namelist_cfg
@@ -13,9 +13,9 @@
 &namusr_def    !   User defined :   OVERFLOW configuration
 !-----------------------------------------------------------------------
    !                       !  type of vertical coordinate
-   ln_zco      = .false.      ! z-coordinate
-   ln_zps      = .true.       ! z-partial-step coordinate
-   ln_sco      = .false.      ! s-coordinate
+   l_zco       = .false.      ! z-coordinate
+   l_zps       = .true.       ! z-partial-step coordinate
+   l_sco       = .false.      ! s-coordinate
    rn_dx       =   5000.   !  horizontal resolution   [meters]
    rn_dz       =     60.   !  vertical   resolution   [meters]
    rn_f0       =  1.e-4    !  coriolis [s-1]
diff --git a/tests/DOME/EXPREF/namelist_cfg b/tests/DOME/EXPREF/namelist_cfg
index 2e9c8bea6a13c03a9feb3cf4c3e6665ce54a7c1d..4f081423fd569acc924cd8746736fab32b10061b 100644
--- a/tests/DOME/EXPREF/namelist_cfg
+++ b/tests/DOME/EXPREF/namelist_cfg
@@ -6,9 +6,9 @@
 &namusr_def    !   User defined :   OVERFLOW configuration
 !-----------------------------------------------------------------------
    !                       !  type of vertical coordinate
-   ln_zco      = .false.      ! z-coordinate
-   ln_zps      = .true.       ! z-partial-step coordinate
-   ln_sco      = .false.      ! s-coordinate
+   l_zco       = .false.      ! z-coordinate
+   l_zps       = .true.       ! z-partial-step coordinate
+   l_sco       = .false.      ! s-coordinate
    rn_dx       =   5000.   !  horizontal resolution   [meters]
    rn_dz       =     60.   !  vertical   resolution   [meters]
    rn_f0       =  1.e-4    !  coriolis [s-1]
diff --git a/tests/DOME/MY_SRC/usrdef_nam.F90 b/tests/DOME/MY_SRC/usrdef_nam.F90
index 30b61501a9650d5706f383efc59a2c452d7ee156..c73bf1158f997ff1ba234b9c2f0222e8846e8a1b 100644
--- a/tests/DOME/MY_SRC/usrdef_nam.F90
+++ b/tests/DOME/MY_SRC/usrdef_nam.F90
@@ -62,7 +62,7 @@ CONTAINS
       INTEGER ::   ighost_w, ighost_e, ighost_s, ighost_n
       REAL(wp)::   zlx, zly, zh ! Local scalars
       !!
-      NAMELIST/namusr_def/  ln_zco, ln_zps, ln_sco, rn_dx, rn_dz, rn_f0
+      NAMELIST/namusr_def/  l_zco, l_zps, l_sco, rn_dx, rn_dz, rn_f0
       !!----------------------------------------------------------------------
       !
       READ  ( numnam_cfg, namusr_def, IOSTAT = ios, ERR = 902 )
@@ -122,9 +122,9 @@ CONTAINS
          WRITE(numout,*) 'usr_def_nam  : read the user defined namelist (namusr_def) in namelist_cfg'
          WRITE(numout,*) '~~~~~~~~~~~ '
          WRITE(numout,*) '   Namelist namusr_def : DOME test case'
-         WRITE(numout,*) '      z-coordinate flag                 ln_zco = ', ln_zco
-         WRITE(numout,*) '      z-partial-step coordinate flag    ln_zps = ', ln_zps 
-         WRITE(numout,*) '      s-coordinate flag                 ln_sco = ', ln_sco  
+         WRITE(numout,*) '      z-coordinate flag                 l_zco  = ', l_zco
+         WRITE(numout,*) '      z-partial-step coordinate flag    l_zps  = ', l_zps 
+         WRITE(numout,*) '      s-coordinate flag                 l_sco  = ', l_sco  
          WRITE(numout,*) '      horizontal resolution             rn_dx  = ', rn_dx, ' m'
          WRITE(numout,*) '      vertical resolution               rn_dz  = ', rn_dz, ' m'
          WRITE(numout,*) '      resulting global domain size :    Ni0glo = ', kpi
diff --git a/tests/DOME/MY_SRC/usrdef_zgr.F90 b/tests/DOME/MY_SRC/usrdef_zgr.F90
index 0d4e2794ee695ab1db0e25772cd7970c81d12bff..19bdf020813230e4707661f3fd766add500ad38c 100644
--- a/tests/DOME/MY_SRC/usrdef_zgr.F90
+++ b/tests/DOME/MY_SRC/usrdef_zgr.F90
@@ -48,7 +48,7 @@ CONTAINS
       !! ** Purpose :   User defined the vertical coordinates
       !!
       !!----------------------------------------------------------------------
-      LOGICAL                             , INTENT(out) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags
+      LOGICAL                             , INTENT( in) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags
       LOGICAL                             , INTENT(out) ::   ld_isfcav                   ! under iceshelf cavity flag
       INTEGER , DIMENSION(:,:)            , INTENT(out) ::   k_top, k_bot                ! first & last ocean level
       REAL(wp), DIMENSION(:)              , INTENT(out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth     [m]
@@ -73,7 +73,6 @@ CONTAINS
       ! ---------------------------
       ! already set in usrdef_nam.F90 by reading the namusr_def namelist except for ISF
       ld_isfcav = .FALSE.
-      ld_sco = .TRUE.
       !
       !
       ! Build the vertical coordinate system
@@ -87,7 +86,7 @@ CONTAINS
       zhu(:,:) = 600_wp ; zhv(:,:) = 600_wp ; zhf(:,:) = 600_wp
       DO_2D( 0, 0, 0, 0 )
          zhu(ji,jj) =  0.5_wp * ( zht(ji,jj  ) + zht(ji+1,jj  ) )
-         zhv(jj,jj) =  0.5_wp * ( zht(ji,jj  ) + zht(ji  ,jj+1) )
+         zhv(ji,jj) =  0.5_wp * ( zht(ji,jj  ) + zht(ji  ,jj+1) )
          zhf(ji,jj) = 0.25_wp * ( zht(ji,jj  ) + zht(ji+1,jj  ) &
                        &        + zht(ji,jj+1) + zht(ji+1,jj+1) ) 
       END_2D
diff --git a/tests/DOME/cpp_DOME.fcm b/tests/DOME/cpp_DOME.fcm
index 63c8f248df6412197ae383980fd285179ae59c26..81fe7c6363b67cd6480e4d19ad202f6ed8000e82 100644
--- a/tests/DOME/cpp_DOME.fcm
+++ b/tests/DOME/cpp_DOME.fcm
@@ -1 +1 @@
- bld::tool::fppkeys key_xios key_agrif key_vco_3d
+ bld::tool::fppkeys key_xios key_agrif key_linssh key_vco_3d
diff --git a/tests/ISOMIP+/cpp_ISOMIP+.fcm b/tests/ISOMIP+/cpp_ISOMIP+.fcm
index e2ce910a38c73d422f72d322746ebcddc3c8ad42..772938ffcb020fb206165a5b5b15e931f93e6ae2 100644
--- a/tests/ISOMIP+/cpp_ISOMIP+.fcm
+++ b/tests/ISOMIP+/cpp_ISOMIP+.fcm
@@ -1 +1 @@
- bld::tool::fppkeys   key_xios key_qco key_vco_1d key_isf 
+ bld::tool::fppkeys   key_xios key_qco key_vco_1d3d key_isf 
diff --git a/tests/VORTEX/EXPREF/1_namelist_cfg b/tests/VORTEX/EXPREF/1_namelist_cfg
index 6d007b1abd2e9622ba28f25d533698a437ed68df..11aafe76d221af11dfae6ec8b27c6b340ec8f43e 100644
--- a/tests/VORTEX/EXPREF/1_namelist_cfg
+++ b/tests/VORTEX/EXPREF/1_namelist_cfg
@@ -1,5 +1,5 @@
 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-!! NEMO/OCE  Configuration namelist : overwrite defaults values defined in SHARED/namelist_ref
+!! NEMO/OCE  Configuration namelist : overwrite default values defined in SHARED/namelist_ref
 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 !!                         VORTEX configuration                       !!
 !!======================================================================
@@ -19,10 +19,12 @@
 !-----------------------------------------------------------------------
 &namusr_def    !   User defined :   VORTEX configuration: Flat bottom, beta-plane
 !-----------------------------------------------------------------------
-   rn_dx       =  30000.   !  x horizontal resolution   [meters]
-   rn_dy       =  30000.   !  y horizontal resolution   [meters]
+   rn_dx       =  10000.   !  x horizontal resolution   [meters]
+   rn_dy       =  10000.   !  y horizontal resolution   [meters]
    rn_dz       =    500.   !  z vertical resolution [meters]
    rn_ppgphi0  =   38.5    !  Reference latitude [degrees]
+   rn_ppumax   =    1.0    !  Max velocity scale [m/s]
+   nn_rot      =      0    !  = 0,1,2,3 domain quarter turn clockwise rotation 
 /
 !-----------------------------------------------------------------------
 &namrun        !   parameters of the run
@@ -39,11 +41,12 @@
 &namdom        !   time and space domain
 !-----------------------------------------------------------------------
    ln_linssh   =  .false.  !  =T  linear free surface  ==>>  model level are fixed in time
-   rn_Dt      =    480.   !  time step for the dynamics (and tracer if nn_acc=0)
+   rn_Dt       =    480.   !  time step for the dynamics (and tracer if nn_acc=0)
    rn_atfp     =   0.05    !  asselin time filter parameter
+   ln_meshmask = .false.   !  =T create a mesh file
 /
 !-----------------------------------------------------------------------
-&namcfg        !   parameters of the configuration                      (default: user defined GYRE)
+&namcfg        !   parameters of the configuration                      (default: use namusr_def in namelist_cfg)
 !-----------------------------------------------------------------------
 /
 !-----------------------------------------------------------------------
@@ -70,7 +73,7 @@
 !!======================================================================
 !
 !-----------------------------------------------------------------------
-&namsbc        !   Surface Boundary Condition (surface module)
+&namsbc        !   Surface Boundary Condition manager                   (default: NO selection)
 !-----------------------------------------------------------------------
   nn_fsbc     = 1         !  frequency of surface boundary condition computation
                           !     (also = the frequency of sea-ice & iceberg model call)
@@ -117,7 +120,7 @@
 !-----------------------------------------------------------------------
 &namdrg        !   top/bottom drag coefficient                          (default: NO selection)
 !-----------------------------------------------------------------------
-   ln_drg_OFF = .true.     !  free-slip       : Cd = 0                  (F => fill namdrg_bot
+   ln_drg_OFF  = .true.   !  free-slip       : Cd = 0                  (F => fill namdrg_bot
 /
 !!======================================================================
 !!                        Tracer (T-S) namelists                      !!
@@ -131,17 +134,17 @@
 !!======================================================================
 !
 !-----------------------------------------------------------------------
-&nameos        !   ocean Equation Of Seawater                           (default: OFF)
+&nameos        !   ocean Equation Of Seawater                           (default: NO selection)
 !-----------------------------------------------------------------------
    ln_seos     = .true.         !  = Use simplified equation of state (S-EOS)
-      !                            !  rd(T,S,Z)*rho0 = -a0*(1+.5*lambda*dT+mu*Z+nu*dS)*dT+b0*dS
-      rn_a0       =  0.28        !  thermal expension coefficient (for simplified equation of state)
-      rn_b0       =  0.          !  saline  expension coefficient (for simplified equation of state)
-      rn_lambda1  =  0.          !  cabbeling coeff in T^2  (=0 for linear eos)
-      rn_lambda2  =  0.          !  cabbeling coeff in S^2  (=0 for linear eos)
-      rn_mu1      =  0.          !  thermobaric coeff. in T (=0 for linear eos)
-      rn_mu2      =  0.          !  thermobaric coeff. in S (=0 for linear eos)
-      rn_nu       =  0.          !  cabbeling coeff in T*S  (=0 for linear eos)
+   !                            !  rd(T,S,Z)*rho0 = -a0*(1+.5*lambda*dT+mu*Z+nu*dS)*dT+b0*dS
+   rn_a0       =  0.28        !  thermal expension coefficient (for simplified equation of state)
+   rn_b0       =  0.          !  saline  expension coefficient (for simplified equation of state)
+   rn_lambda1  =  0.          !  cabbeling coeff in T^2  (=0 for linear eos)
+   rn_lambda2  =  0.          !  cabbeling coeff in S^2  (=0 for linear eos)
+   rn_mu1      =  0.          !  thermobaric coeff. in T (=0 for linear eos)
+   rn_mu2      =  0.          !  thermobaric coeff. in S (=0 for linear eos)
+   rn_nu       =  0.          !  cabbeling coeff in T*S  (=0 for linear eos)
 /
 !-----------------------------------------------------------------------
 &namtra_adv    !   advection scheme for tracer                          (default: NO selection)
@@ -149,7 +152,7 @@
    ln_traadv_cen = .false. !  2nd order centered scheme
       nn_cen_h   =  4            !  =2/4, horizontal 2nd order CEN / 4th order CEN
       nn_cen_v   =  4            !  =2/4, vertical   2nd order CEN / 4th order COMPACT
-   ln_traadv_fct = .true.  !  FCT scheme
+   ln_traadv_fct = .true. !  FCT scheme
       nn_fct_h   =  2            !  =2/4, horizontal 2nd / 4th order
       nn_fct_v   =  2            !  =2/4, vertical   2nd / COMPACT 4th order
    ln_traadv_mus = .false. !  MUSCL scheme
@@ -161,7 +164,6 @@
 !-----------------------------------------------------------------------
 &namtra_ldf    !   lateral diffusion scheme for tracers                 (default: NO selection)
 !-----------------------------------------------------------------------
-   !                       !  Operator type:
    ln_traldf_OFF   =  .true.  !  No explicit diffusion
    ln_traldf_lap   =  .false.  !    laplacian operator
    ln_traldf_blp   =  .false.  !  bilaplacian operator
@@ -208,7 +210,7 @@
 !-----------------------------------------------------------------------
 &namdyn_spg    !   surface pressure gradient                            (default: OFF)
 !-----------------------------------------------------------------------
-   ln_dynspg_exp  = .false. 
+   ln_dynspg_exp = .false. 
    ln_dynspg_ts  = .true.   ! split-explicit free surface
       ln_bt_fw      = .true.     ! Forward integration of barotropic Eqs.
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
@@ -249,7 +251,7 @@
 !!======================================================================
 !
 !-----------------------------------------------------------------------
-&namzdf        !   vertical physics                                     (default: NO selection)
+&namzdf        !   vertical physics manager                             (default: NO selection)
 !-----------------------------------------------------------------------
    !                       ! type of vertical closure
    ln_zdfcst   = .true.       !  constant mixing
diff --git a/tests/VORTEX/EXPREF/namelist_cfg b/tests/VORTEX/EXPREF/namelist_cfg
index 99b9b21a9e8ab0694c105fdda20a440e36516e0f..b44c62ab22e2bcc9ddcb1aeaea7ec3aee8a36379 100644
--- a/tests/VORTEX/EXPREF/namelist_cfg
+++ b/tests/VORTEX/EXPREF/namelist_cfg
@@ -23,6 +23,8 @@
    rn_dy       =  30000.   !  y horizontal resolution   [meters]
    rn_dz       =    500.   !  z vertical resolution [meters]
    rn_ppgphi0  =   38.5    !  Reference latitude [degrees]
+   rn_ppumax   =    1.0    !  Max velocity scale [m/s]
+   nn_rot      =      0    !  = 0,1,2,3 domain quarter turn clockwise rotation 
 /
 !-----------------------------------------------------------------------
 &namrun        !   parameters of the run
@@ -39,8 +41,9 @@
 &namdom        !   time and space domain
 !-----------------------------------------------------------------------
    ln_linssh   =  .false.  !  =T  linear free surface  ==>>  model level are fixed in time
-   rn_Dt      =   1440.   !  time step for the dynamics (and tracer if nn_acc=0)
+   rn_Dt       =   1440.   !  time step for the dynamics (and tracer if nn_acc=0)
    rn_atfp     =   0.05    !  asselin time filter parameter
+   ln_meshmask = .false.   !  =T create a mesh file
 /
 !-----------------------------------------------------------------------
 &namcfg        !   parameters of the configuration                      (default: use namusr_def in namelist_cfg)
diff --git a/tests/VORTEX/MY_SRC/usrdef_hgr.F90 b/tests/VORTEX/MY_SRC/usrdef_hgr.F90
index c6d0003e8d1ad92018c9dcd946dc0024fabf6767..60e09483c83c8337b566b20137817138ab53ab36 100644
--- a/tests/VORTEX/MY_SRC/usrdef_hgr.F90
+++ b/tests/VORTEX/MY_SRC/usrdef_hgr.F90
@@ -23,8 +23,6 @@ MODULE usrdef_hgr
    IMPLICIT NONE
    PRIVATE
 
-   REAL(wp) :: roffsetx, roffsety ! Offset in km to first f-point
-
    PUBLIC   usr_def_hgr   ! called by domhgr.F90
 
    !! * Substitutions
@@ -65,8 +63,9 @@ CONTAINS
       REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe1e2u, pe1e2v               ! u- & v-surfaces (if reduction in strait)   [m2]
       !
       INTEGER  ::   ji, jj     ! dummy loop indices
-      REAL(wp) ::   zbeta, zf0
-      REAL(wp) ::   zti, ztj   ! local scalars
+      REAL(wp) ::   zbeta, zf0 ! local scalars
+      REAL(wp) ::   zroffsetx, zroffsety
+      REAL(wp) ::   zti, ztj   
       !!-------------------------------------------------------------------------------
       !
       IF(lwp) WRITE(numout,*)
@@ -81,39 +80,88 @@ CONTAINS
       ! offset is given at first f-point, i.e. at (i,j) = (nn_hls+1, nn_hls+1)
       ! Here we assume the grid is centred around a T-point at the middle of
       ! of the domain (hence domain size is odd) 
-      roffsetx = (-REAL(Ni0glo-1, wp) + 1._wp) * 0.5 * 1.e-3 * rn_dx
-      roffsety = (-REAL(Nj0glo-1, wp) + 1._wp) * 0.5 * 1.e-3 * rn_dy
+      zroffsetx = (-REAL(Ni0glo-1, wp) + 1._wp) * 0.5_wp * 1.e-3 * rn_dx
+      zroffsety = (-REAL(Nj0glo-1, wp) + 1._wp) * 0.5_wp * 1.e-3 * rn_dy
 #if defined key_agrif
       IF( .NOT.Agrif_Root() ) THEN
          ! deduce offset from parent:
-         roffsetx = Agrif_Parent(roffsetx) &
-                  & + (-(nbghostcells_x_w - 1) + (Agrif_Parent(nbghostcells_x_w) + Agrif_Ix()-2)*Agrif_Rhox()) * 1.e-3 * rn_dx
-         roffsety = Agrif_Parent(roffsety) &
-                  & + (-(nbghostcells_y_s - 1) + (Agrif_Parent(nbghostcells_y_s) + Agrif_Iy()-2)*Agrif_Rhoy()) * 1.e-3 * rn_dy
+!         zroffsetx = Agrif_Parent(zroffsetx) &
+!                  & + REAL(-(nbghostcells_x_w - 1) + (Agrif_Parent(nbghostcells_x_w) + Agrif_Ix()-2)*Agrif_iRhox(), wp ) * 1.e-3 * rn_dx
+!         zroffsety = Agrif_Parent(zroffsety) &
+!                  & + REAL(-(nbghostcells_y_s - 1) + (Agrif_Parent(nbghostcells_y_s) + Agrif_Iy()-2)*Agrif_iRhoy(), wp ) * 1.e-3 * rn_dy
+! JC: More accurate positioning implied by 1:1 tests:
+!     (10-13 km shift between child and parent grids  otherwise)
+         zroffsetx = ((-REAL(Agrif_Parent(Ni0glo)-1, wp) + 1._wp) * 0.5_wp * Agrif_Rhox() &
+                   & + REAL(-(nbghostcells_x_w - 1) + (Agrif_Parent(nbghostcells_x_w)     & 
+                   & + Agrif_Ix()-2)*Agrif_iRhox(), wp )) * 1.e-3 * rn_dx
+         zroffsety = ((-REAL(Agrif_Parent(Nj0glo)-1, wp) + 1._wp) * 0.5_wp * Agrif_Rhoy() &
+                   & + REAL(-(nbghostcells_y_s - 1) + (Agrif_Parent(nbghostcells_y_s)     & 
+                   & + Agrif_Iy()-2)*Agrif_iRhoy(), wp )) * 1.e-3 * rn_dy
       ENDIF
 #endif         
       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          zti = REAL( mig(ji,0)-1, wp )  ! start at i=0 in the global grid without halos
          ztj = REAL( mjg(jj,0)-1, wp )  ! start at j=0 in the global grid without halos
          
-         plamt(ji,jj) = roffsetx + rn_dx * 1.e-3 * ( zti - 0.5_wp )
-         plamu(ji,jj) = roffsetx + rn_dx * 1.e-3 *   zti 
-         plamv(ji,jj) = plamt(ji,jj) 
-         plamf(ji,jj) = plamu(ji,jj) 
-         
-         pphit(ji,jj) = roffsety + rn_dy * 1.e-3 * ( ztj - 0.5_wp )
-         pphiv(ji,jj) = roffsety + rn_dy * 1.e-3 *   ztj 
-         pphiu(ji,jj) = pphit(ji,jj) 
-         pphif(ji,jj) = pphiv(ji,jj) 
+         ! Longitudes
+         IF ( nn_rot==0 ) THEN
+            plamt(ji,jj) = zroffsetx + rn_dx * 1.e-3 * ( zti - 0.5_wp )
+            plamu(ji,jj) = plamt(ji,jj) + rn_dx * 1.e-3 * 0.5_wp 
+            plamv(ji,jj) = plamt(ji,jj) 
+            plamf(ji,jj) = plamu(ji,jj) 
+         ELSEIF ( nn_rot==1 ) THEN
+            plamt(ji,jj) = -(zroffsety + rn_dy * 1.e-3 * ( ztj - 0.5_wp ))
+            plamu(ji,jj) = plamt(ji,jj)
+            plamv(ji,jj) = plamt(ji,jj) - rn_dy * 1.e-3 * 0.5_wp
+            plamf(ji,jj) = plamv(ji,jj) 
+         ELSEIF ( nn_rot==2 ) THEN
+            plamt(ji,jj) = -(zroffsetx + rn_dx * 1.e-3 * ( zti - 0.5_wp ))
+            plamu(ji,jj) = plamt(ji,jj) - rn_dx * 1.e-3 * 0.5_wp 
+            plamv(ji,jj) = plamt(ji,jj) 
+            plamf(ji,jj) = plamu(ji,jj) 
+         ELSEIF ( nn_rot==3 ) THEN
+            plamt(ji,jj) = (zroffsety + rn_dy * 1.e-3 * ( ztj - 0.5_wp ))
+            plamu(ji,jj) = plamt(ji,jj)
+            plamv(ji,jj) = plamt(ji,jj) + rn_dy * 1.e-3 * 0.5_wp
+            plamf(ji,jj) = plamv(ji,jj) 
+         ENDIF          
+         ! Latitudes:
+         IF ( nn_rot==0 ) THEN
+            pphit(ji,jj) = (zroffsety + rn_dy * 1.e-3 * ( ztj - 0.5_wp ))
+            pphiv(ji,jj) = pphit(ji,jj) + rn_dy * 1.e-3 * 0.5_wp 
+            pphiu(ji,jj) = pphit(ji,jj) 
+            pphif(ji,jj) = pphiv(ji,jj) 
+         ELSEIF ( nn_rot==1 ) THEN
+            pphit(ji,jj) = (zroffsetx + rn_dx * 1.e-3 * ( zti - 0.5_wp ))
+            pphiv(ji,jj) = pphit(ji,jj) 
+            pphiu(ji,jj) = pphit(ji,jj) + rn_dx * 1.e-3 * 0.5_wp 
+            pphif(ji,jj) = pphiu(ji,jj) 
+         ELSEIF ( nn_rot==2 ) THEN
+            pphit(ji,jj) = -(zroffsety + rn_dy * 1.e-3 * ( ztj - 0.5_wp ))
+            pphiv(ji,jj) = pphit(ji,jj) - rn_dy * 1.e-3 * 0.5_wp 
+            pphiu(ji,jj) = pphit(ji,jj) 
+            pphif(ji,jj) = pphiv(ji,jj) 
+         ELSEIF ( nn_rot==3 ) THEN
+            pphit(ji,jj) = -(zroffsetx + rn_dx * 1.e-3 * ( zti - 0.5_wp ))
+            pphiv(ji,jj) = pphit(ji,jj) 
+            pphiu(ji,jj) = pphit(ji,jj) - rn_dx * 1.e-3 * 0.5_wp 
+            pphif(ji,jj) = pphiu(ji,jj) 
+         ENDIF
       END_2D
       !     
       ! Horizontal scale factors (in meters)
       !                              ======
-      pe1t(:,:) = rn_dx  ;   pe2t(:,:) = rn_dy 
-      pe1u(:,:) = rn_dx  ;   pe2u(:,:) = rn_dy 
-      pe1v(:,:) = rn_dx  ;   pe2v(:,:) = rn_dy 
-      pe1f(:,:) = rn_dx  ;   pe2f(:,:) = rn_dy 
-
+      IF     ((nn_rot==0).OR.(nn_rot==2)) THEN
+         pe1t(:,:) = rn_dx  ;   pe2t(:,:) = rn_dy 
+         pe1u(:,:) = rn_dx  ;   pe2u(:,:) = rn_dy 
+         pe1v(:,:) = rn_dx  ;   pe2v(:,:) = rn_dy 
+         pe1f(:,:) = rn_dx  ;   pe2f(:,:) = rn_dy 
+      ELSEIF ((nn_rot==1).OR.(nn_rot==3)) THEN
+         pe1t(:,:) = rn_dy  ;   pe2t(:,:) = rn_dx 
+         pe1u(:,:) = rn_dy  ;   pe2u(:,:) = rn_dx 
+         pe1v(:,:) = rn_dy  ;   pe2v(:,:) = rn_dx 
+         pe1f(:,:) = rn_dy  ;   pe2f(:,:) = rn_dx 
+      ENDIF
       !                             ! NO reduction of grid size in some straits 
       ke1e2u_v = 0                  !    ==>> u_ & v_surfaces will be computed in dom_hgr routine
       pe1e2u(:,:) = 0._wp           !    CAUTION: set to zero to avoid error with some compilers that
diff --git a/tests/VORTEX/MY_SRC/usrdef_istate.F90 b/tests/VORTEX/MY_SRC/usrdef_istate.F90
index 7694c43c9c5ec788f0a943a2615b431b2e0a1b4c..b53fab3392df55880fab5e8a0bb59786b0c127e1 100644
--- a/tests/VORTEX/MY_SRC/usrdef_istate.F90
+++ b/tests/VORTEX/MY_SRC/usrdef_istate.F90
@@ -67,8 +67,8 @@ CONTAINS
       !
       !
       zf0   = 2._wp * omega * SIN( rad * rn_ppgphi0 )
-      zumax = 1._wp * SIGN(1._wp, zf0) ! Here Anticyclonic: set zumax=-1 for cyclonic
-      zlambda = SQRT(2._wp)*60.e3      ! Horizontal scale in meters
+      zumax = rn_ppumax * SIGN(1._wp, zf0) ! Here Anticyclonic: set zumax=-1 for cyclonic
+      zlambda = SQRT(2._wp)*60.e3          ! Horizontal scale in meters
       zn2 = 3.e-3**2
       zH = 0.5_wp * 5000._wp
       !
@@ -101,7 +101,15 @@ CONTAINS
             zdu = 0.5_wp * (pdept(ji  ,jj,jk) + pdept(ji+1,jj,jk))
             IF (zdu < zH) THEN
                zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH))
-               pu(ji,jj,jk) = (za * zf * zy * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk)
+               IF ( nn_rot==0 ) THEN
+                  pu(ji,jj,jk) = (za * zf * zy * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk)
+               ELSEIF ( nn_rot==1 ) THEN
+                  pu(ji,jj,jk) = -(za * zf * zx * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk)
+               ELSEIF ( nn_rot==2 ) THEN
+                  pu(ji,jj,jk) = -(za * zf * zy * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk)
+               ELSEIF ( nn_rot==3 ) THEN
+                  pu(ji,jj,jk) =  (za * zf * zx * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk)
+               ENDIF
             ELSE
                pu(ji,jj,jk) = 0._wp
             ENDIF
@@ -115,7 +123,15 @@ CONTAINS
             zdv = 0.5_wp * (pdept(ji  ,jj,jk) + pdept(ji,jj+1,jk))
             IF (zdv < zH) THEN
                zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH))
-               pv(ji,jj,jk) = -(za * zf * zx * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk)
+               IF ( nn_rot==0 ) THEN
+                  pv(ji,jj,jk) = -(za * zf * zx * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk)
+               ELSEIF ( nn_rot==1 ) THEN
+                  pv(ji,jj,jk) = -(za * zf * zy * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk)
+               ELSEIF ( nn_rot==2 ) THEN
+                  pv(ji,jj,jk) =  (za * zf * zx * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk)
+               ELSEIF ( nn_rot==3 ) THEN
+                  pv(ji,jj,jk) =  (za * zf * zy * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk)
+               ENDIF
             ELSE
                pv(ji,jj,jk) = 0._wp
             ENDIF
@@ -151,8 +167,8 @@ CONTAINS
       !
       !
       zf0   = 2._wp * omega * SIN( rad * rn_ppgphi0 )
-      zumax = 1._wp * SIGN(1._wp, zf0) ! Here Anticyclonic: set zumax=-1 for cyclonic
-      zlambda = SQRT(2._wp)*60.e3      ! Horizontal scale in meters 
+      zumax = rn_ppumax * SIGN(1._wp, zf0) ! Here Anticyclonic: set zumax=-1 for cyclonic
+      zlambda = SQRT(2._wp)*60.e3          ! Horizontal scale in meters 
       zH = 0.5_wp * 5000._wp
       !
       zP0 = rho0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp)
diff --git a/tests/VORTEX/MY_SRC/usrdef_nam.F90 b/tests/VORTEX/MY_SRC/usrdef_nam.F90
index 432fe596057892aac914e612b5e48e01642a2ed4..a8f5ff69d028e680106dcf6dcbc6c6c8a0a10d26 100644
--- a/tests/VORTEX/MY_SRC/usrdef_nam.F90
+++ b/tests/VORTEX/MY_SRC/usrdef_nam.F90
@@ -31,6 +31,8 @@ MODULE usrdef_nam
    REAL(wp), PUBLIC ::   rn_dy      ! resolution in meters defining the horizontal domain size
    REAL(wp), PUBLIC ::   rn_dz      ! vertical resolution 
    REAL(wp), PUBLIC ::   rn_ppgphi0 ! reference latitude for beta-plane 
+   REAL(wp), PUBLIC ::   rn_ppumax  ! vortex velocity 
+   INTEGER,  PUBLIC ::   nn_rot     ! vortex velocity 
 
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -62,12 +64,14 @@ CONTAINS
       INTEGER :: ighost_n, ighost_s, ighost_w, ighost_e
       REAL(wp)::   zlx, zly, zh ! Local scalars
       !!
-      NAMELIST/namusr_def/  rn_dx, rn_dy, rn_dz, rn_ppgphi0
+      NAMELIST/namusr_def/  rn_dx, rn_dy, rn_dz, rn_ppgphi0, rn_ppumax, nn_rot
       !!----------------------------------------------------------------------
       !
       READ  ( numnam_cfg, namusr_def, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namusr_def in configuration namelist' )
       !
+      IF( nn_rot  < 0   .OR.  nn_rot  > 3 )   CALL ctl_stop( 'bad namelist flag: nn_rot is 0, 1, 2 or 3' )
+      !
 #if defined key_agrif 
       ! Domain parameters are taken from parent:
       IF( .NOT. Agrif_Root() ) THEN
@@ -75,6 +79,8 @@ CONTAINS
          rn_dy = Agrif_Parent(rn_dy)/Agrif_Rhoy()
 !         rn_dz = Agrif_Parent(rn_dz)
          rn_ppgphi0 = Agrif_Parent(rn_ppgphi0)
+         rn_ppumax  = Agrif_Parent(rn_ppumax)
+         nn_rot     = Agrif_Parent(nn_rot)
       ENDIF
 #endif
       !
@@ -86,8 +92,13 @@ CONTAINS
 #if defined key_agrif 
       IF( Agrif_Root() ) THEN       ! Global Domain size:  VORTEX global domain is  1800 km x 1800 Km x 5000 m
 #endif
-         kpi = NINT( 1800.e3  / rn_dx ) + 3  
-         kpj = NINT( 1800.e3  / rn_dy ) + 3 
+         IF     ((nn_rot==0).OR.(nn_rot==2)) THEN
+            kpi = NINT( 1800.e3  / rn_dx ) + 3  
+            kpj = NINT( 1800.e3  / rn_dy ) + 3 
+         ELSEIF ((nn_rot==1).OR.(nn_rot==3)) THEN
+            kpi = NINT( 1800.e3  / rn_dy ) + 3  
+            kpj = NINT( 1800.e3  / rn_dx ) + 3 
+         ENDIF
 #if defined key_agrif 
       ELSE                          ! Global Domain size: add nbghostcells + 1 "land" point on each side
          ! At this stage, child ghosts have not been set
@@ -133,6 +144,8 @@ CONTAINS
          WRITE(numout,*) '         LY [km]: ', zly
          WRITE(numout,*) '          H [m] : ', zh
          WRITE(numout,*) '      Reference latitude            rn_ppgphi0 = ', rn_ppgphi0
+         WRITE(numout,*) '      Max. current amplitude        rn_ppumax  = ', rn_ppumax
+         WRITE(numout,*) '      Domain rotation               nn_rot     = ', nn_rot
          WRITE(numout,*) '   '
       ENDIF
       !
diff --git a/tools/DOMAINcfg/src/domhgr.F90 b/tools/DOMAINcfg/src/domhgr.F90
index fc30c5d1d92bea908dd08485750ff3d259148e1f..f768835728266fba9aef09ae878e37c180947206 100644
--- a/tools/DOMAINcfg/src/domhgr.F90
+++ b/tools/DOMAINcfg/src/domhgr.F90
@@ -195,6 +195,11 @@ CONTAINS
             gphi0 = gphi0 -  800._wp 
          ENDIF
          !
+         IF ( cp_cfg=='VORTEX' ) THEN
+            glam0 = glam0 + (-REAL(Ni0glo-1, wp) + 1._wp) * 0.5 * 1.e-3 * ppe1_m    
+            gphi0 = gphi0 + (-REAL(Nj0glo-1, wp) + 1._wp) * 0.5 * 1.e-3 * ppe2_m    
+         ENDIF
+         !
          DO jj = 1, jpj
             DO ji = 1, jpi
                glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 *  REAL( mig0(ji)-1   ) 
@@ -389,6 +394,7 @@ CONTAINS
          !
          zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
          zphi0   = ppgphi0 - REAL( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
+         IF ( cp_cfg=='VORTEX' ) zphi0 = ppgphi0 
          !
          zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
          !
diff --git a/tools/DOMAINcfg/src/domzgr.F90 b/tools/DOMAINcfg/src/domzgr.F90
index 065ac376e1576c6003b9da3890691912e31ea765..38c8556de8fd502f1afdb1db0fb4fb836b94deb6 100644
--- a/tools/DOMAINcfg/src/domzgr.F90
+++ b/tools/DOMAINcfg/src/domzgr.F90
@@ -574,7 +574,7 @@ CONTAINS
             CALL RANDOM_SEED()
             CALL RANDOM_NUMBER(zrand)
             DO_2D( 0, 0, 0, 0 )
-               bathy(ji,jj) = h_oce + 0.1_wp *h_oce * (zrand(mig(ji),mjg(jj))-1._wp) 
+               bathy(ji,jj) = h_oce + 0.3_wp *h_oce * (zrand(mig(ji),mjg(jj))-1._wp) 
             END_2D
             IF ( cp_cfg=='OVERFLOW' ) THEN
                DO jj=1,jpj
diff --git a/tools/DOMAINcfg/tests/VORTEX/1_namelist_cfg b/tools/DOMAINcfg/tests/VORTEX/1_namelist_cfg
new file mode 100644
index 0000000000000000000000000000000000000000..d0b194f2439b5a95571211d2485e421b699a3942
--- /dev/null
+++ b/tools/DOMAINcfg/tests/VORTEX/1_namelist_cfg
@@ -0,0 +1,94 @@
+!! NEMO/OCE :   Configuration namelist_cfg used to overwrite defaults value defined in namelist_ref
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!! NEMO/OCE  :  1 - Domain & run manager (namrun, namcfg, namdom, namzgr, namzgr_sco )
+!!              2 - diagnostics      (namnc4)
+!!              3 - miscellaneous    (nammpp, namctl)
+!!
+!! namelist skeleton : egrep -E '(^/ *$|^! *$|^ *$|&nam.*|!---.*|!! .*|!!==.*|!!>>>.*)' namelist_ref > namelist_skl
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!-----------------------------------------------------------------------
+&namrun        !   parameters of the run
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namdom        !   space and time domain (bathymetry, mesh, timestep)
+!-----------------------------------------------------------------------
+   ln_read_cfg = .false.
+   nn_bathy = -1 
+                           !  or compute (2) from external bathymetry
+   jphgr_msh   =       3   !  f-plan 
+   ppgphi0     =      38.5 !  reference latitude
+   ppe1_m      =   10000.0 !  zonal      grid-spacing (degrees)
+   ppe2_m      =   10000.0 !  meridional grid-spacing (degrees)
+   ppkth       =       0.0 ! if =0. assume uniform grid over pphmax meters
+   pphmax      =    5000.0 !  Maximum depth
+   rn_e3zps_min=    99999. !  partial step thickness is set larger than the minimum of
+   rn_e3zps_rat=    0.1    !  rn_e3zps_min and rn_e3zps_rat*e3t, with 0<rn_e3zps_rat<1
+/
+!-----------------------------------------------------------------------
+&namcfg        !   parameters of the configuration
+!-----------------------------------------------------------------------
+   !
+   ln_e3_dep   = .true.    ! =T : e3=dk[depth] in discret sens.
+   !                       !      ===>>> will become the only possibility in v4.0
+   !                       ! =F : e3 analytical derivative of depth function
+   !                       !      only there for backward compatibility test with v3.6
+      !                      ! if ln_e3_dep = T
+      ln_dept_mid = .true.   ! =T : set T points in the middle of cells
+   !                       !
+   cp_cfg      =   "VORTEX"   !  name of the configuration
+   jp_cfg      =      10   !  resolution of the configuration
+   jpidta = 65
+   jpjdta = 65
+   jpkdta      =      11   !  number of levels      ( >= jpk )
+   Ni0glo = 65
+   Nj0glo = 65
+   jpkglo      =      11
+   jperio = 0
+   ln_use_jattr = .false.  !  use (T) the file attribute: open_ocean_jstart, if present
+                           !  in netcdf input files, as the start j-row for reading
+   ln_domclo = .false.     ! computation of closed sea masks (see namclo)
+/
+!-----------------------------------------------------------------------
+&namzgr        !   vertical coordinate                                  (default: NO selection)
+!-----------------------------------------------------------------------
+!-----------------------------------------------------------------------
+   ln_zco      = .false.   !  z-coordinate - full    steps
+   ln_zps      = .true.    !  z-coordinate - partial steps
+   ln_sco      = .false.   !  s- or hybrid z-s-coordinate
+   ln_isfcav   = .false.   !  ice shelf cavity             (T: see namzgr_isf)
+/
+!-----------------------------------------------------------------------
+&namzgr_isf    !   isf cavity geometry definition
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namzgr_sco    !   s-coordinate or hybrid z-s-coordinate                (default F)
+!-----------------------------------------------------------------------
+   ln_s_sh94   = .true.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)|
+                           !  stretching coefficients for all functions
+   rn_sbot_min =  600.0    !  minimum depth of s-bottom surface (>0) (m)
+   rn_sbot_max = 3600.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
+                        !!!!!!!  Envelop bathymetry
+   rn_rmax     =    1.0    !  maximum cut-off r-value allowed (0<r_max<1)
+                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.)
+   rn_theta    =    0.0    !  surface control parameter (0<=theta<=20)
+   rn_bb       =    0.8    !  stretching with SH94 s-sigma
+/
+!-----------------------------------------------------------------------
+&namclo ! (closed sea : need ln_domclo = .true. in namcfg)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namlbc        !   lateral momentum boundary condition                  (default: NO selection)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namagrif      !  AGRIF zoom                                            ("key_agrif")
+!-----------------------------------------------------------------------
+npt_copy      = 4 
+/
+!-----------------------------------------------------------------------
+&namnc4        !   netcdf4 chunking and compression settings            ("key_netcdf4")
+!-----------------------------------------------------------------------
+/
diff --git a/tools/DOMAINcfg/tests/VORTEX/1_namelist_ref b/tools/DOMAINcfg/tests/VORTEX/1_namelist_ref
new file mode 120000
index 0000000000000000000000000000000000000000..10614842eabd73c6eb44bfa6f75cd7e69e82f7fd
--- /dev/null
+++ b/tools/DOMAINcfg/tests/VORTEX/1_namelist_ref
@@ -0,0 +1 @@
+../../namelist_ref
\ No newline at end of file
diff --git a/tools/DOMAINcfg/tests/VORTEX/AGRIF_FixedGrids.in b/tools/DOMAINcfg/tests/VORTEX/AGRIF_FixedGrids.in
new file mode 120000
index 0000000000000000000000000000000000000000..e8cf8bf1382bd3e00005731dd31c921dfd9a6b09
--- /dev/null
+++ b/tools/DOMAINcfg/tests/VORTEX/AGRIF_FixedGrids.in
@@ -0,0 +1 @@
+../../../../tests/VORTEX/EXPREF/AGRIF_FixedGrids.in
\ No newline at end of file
diff --git a/tools/DOMAINcfg/tests/VORTEX/make_domain_cfg.exe b/tools/DOMAINcfg/tests/VORTEX/make_domain_cfg.exe
new file mode 120000
index 0000000000000000000000000000000000000000..8d81caa2096415e5443d5837091c581ef8fa47d0
--- /dev/null
+++ b/tools/DOMAINcfg/tests/VORTEX/make_domain_cfg.exe
@@ -0,0 +1 @@
+../../make_domain_cfg.exe
\ No newline at end of file
diff --git a/tools/DOMAINcfg/tests/VORTEX/namelist_cfg b/tools/DOMAINcfg/tests/VORTEX/namelist_cfg
new file mode 100644
index 0000000000000000000000000000000000000000..7ebb9ddb40c33802ceea7d720b83c25b56101b8a
--- /dev/null
+++ b/tools/DOMAINcfg/tests/VORTEX/namelist_cfg
@@ -0,0 +1,94 @@
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!! NEMO/OCE :   Configuration namelist_cfg used to overwrite defaults value defined in namelist_ref
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!! NEMO/OCE  :  1 - Domain & run manager (namrun, namcfg, namdom, namzgr, namzgr_sco )
+!!              2 - diagnostics      (namnc4)
+!!              3 - miscellaneous    (nammpp, namctl)
+!!
+!! namelist skeleton : egrep -E '(^/ *$|^! *$|^ *$|&nam.*|!---.*|!! .*|!!==.*|!!>>>.*)' namelist_ref > namelist_skl
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!-----------------------------------------------------------------------
+&namrun        !   parameters of the run
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namdom        !   space and time domain (bathymetry, mesh, timestep)
+!-----------------------------------------------------------------------
+   ln_read_cfg = .false.
+   nn_bathy    =    -1     !  compute analyticaly (=0) or read (=1) the bathymetry file
+                           !  or compute (2) from external bathymetry
+   jphgr_msh   =       3   !  f-plan 
+   ppgphi0     =      38.5 !  reference latitude
+   ppe1_m      =   30000.0 !  zonal      grid-spacing (degrees)
+   ppe2_m      =   30000.0 !  meridional grid-spacing (degrees)
+   ppkth       =       0.0 ! if =0. assume uniform grid over pphmax meters
+   pphmax      =    5000.0 !  Maximum depth
+   rn_e3zps_min=    99999. !  partial step thickness is set larger than the minimum of
+   rn_e3zps_rat=    0.1    !  rn_e3zps_min and rn_e3zps_rat*e3t, with 0<rn_e3zps_rat<1
+/
+!-----------------------------------------------------------------------
+&namcfg        !   parameters of the configuration
+!-----------------------------------------------------------------------
+   !
+   ln_e3_dep   = .true.    ! =T : e3=dk[depth] in discret sens.
+   !                       !      ===>>> will become the only possibility in v4.0
+   !                       ! =F : e3 analytical derivative of depth function
+   !                       !      only there for backward compatibility test with v3.6
+      !                      ! if ln_e3_dep = T
+      ln_dept_mid = .true.   ! =T : set T points in the middle of cells
+   !                       !
+   cp_cfg      =   "VORTEX"   !  name of the configuration
+   jp_cfg      =      30   !  resolution of the configuration
+   jpidta      =      63   !  1st lateral dimension ( >= jpi )
+   jpjdta      =      63   !  2nd    "         "    ( >= jpj )
+   jpkdta      =      11   !  number of levels      ( >= jpk )
+   Ni0glo      =      63   !  1st dimension of global domain --> i =jpidta
+   Nj0glo      =      63   !  2nd    -                  -    --> j  =jpjdta
+   jpkglo      =      11
+   jperio      =       0   !  lateral cond. type (between 0 and 6)
+   ln_use_jattr = .false.  !  use (T) the file attribute: open_ocean_jstart, if present
+                           !  in netcdf input files, as the start j-row for reading
+   ln_domclo = .false.     ! computation of closed sea masks (see namclo)
+/
+!-----------------------------------------------------------------------
+&namzgr        !   vertical coordinate                                  (default: NO selection)
+!-----------------------------------------------------------------------
+!-----------------------------------------------------------------------
+   ln_zco      = .false.   !  z-coordinate - full    steps
+   ln_zps      = .true.    !  z-coordinate - partial steps
+   ln_sco      = .false.   !  s- or hybrid z-s-coordinate
+   ln_isfcav   = .false.   !  ice shelf cavity             (T: see namzgr_isf)
+/
+!-----------------------------------------------------------------------
+&namzgr_isf    !   isf cavity geometry definition
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namzgr_sco    !   s-coordinate or hybrid z-s-coordinate                (default F)
+!-----------------------------------------------------------------------
+   ln_s_sh94   = .true.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)|
+                           !  stretching coefficients for all functions
+   rn_sbot_min =  600.0    !  minimum depth of s-bottom surface (>0) (m)
+   rn_sbot_max = 3600.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
+                        !!!!!!!  Envelop bathymetry
+   rn_rmax     =    1.0    !  maximum cut-off r-value allowed (0<r_max<1)
+                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.)
+   rn_theta    =    0.0    !  surface control parameter (0<=theta<=20)
+   rn_bb       =    0.8    !  stretching with SH94 s-sigma
+/
+!-----------------------------------------------------------------------
+&namclo ! (closed sea : need ln_domclo = .true. in namcfg)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namlbc        !   lateral momentum boundary condition                  (default: NO selection)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namagrif      !  AGRIF zoom                                            ("key_agrif")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namnc4        !   netcdf4 chunking and compression settings            ("key_netcdf4")
+!-----------------------------------------------------------------------
+/
diff --git a/tools/DOMAINcfg/tests/VORTEX/namelist_ref b/tools/DOMAINcfg/tests/VORTEX/namelist_ref
new file mode 120000
index 0000000000000000000000000000000000000000..10614842eabd73c6eb44bfa6f75cd7e69e82f7fd
--- /dev/null
+++ b/tools/DOMAINcfg/tests/VORTEX/namelist_ref
@@ -0,0 +1 @@
+../../namelist_ref
\ No newline at end of file