Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ASM/asminc.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ASM/asminc.F90	(revision 12724)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ASM/asminc.F90	(revision 12731)
@@ -803,5 +803,8 @@
       ELSE
          ALLOCATE( ztim(jpi,jpj) )
-         ztim(:,:) = ssh_iau(:,:) / ( ht(:,:) + 1.0 - ssmask(:,:) )
+!!st         ztim(:,:) = ssh_iau(:,:) / ( ht_(:,:) + 1.0 - ssmask(:,:) )
+         DO_2D_11_11
+            ztim(ji,jj) = ssh_iau(ji,jj) / ( ht(ji,jj) + 1.0 - ssmask(ji,jj) )
+         END_2D
          DO jk = 1, jpkm1
             phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk)
@@ -900,5 +903,5 @@
          IF ( kt == nitdin_r ) THEN
             !
-            l_1st_euler = 0              ! Force Euler forward step
+            l_1st_euler = .TRUE.              ! Force Euler forward step
             !
             ! Sea-ice : SI3 case
@@ -974,6 +977,6 @@
 !           ! set to bottom of a level
 !                 DO jk = jpk-1, 2, -1
-!                   IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN
-!                     mld=gdepw(ji,jj,jk+1)
+!                   IF ((mld > gdepw(ji,jj,jk,Kmm)) .and. (mld < gdepw(ji,jj,jk+1,Kmm))) THEN
+!                     mld=gdepw(ji,jj,jk+1,Kmm)
 !                     jkmax=jk
 !                   ENDIF
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diaar5.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diaar5.F90	(revision 12724)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diaar5.F90	(revision 12731)
@@ -78,5 +78,5 @@
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: zarea_ssh , zbotpres       ! 2D workspace 
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: z2d, zpe                   ! 2D workspace 
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: z3d, zrhd , zrhop, ztpot   ! 3D workspace
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: z3d, zrhd , zrhop, ztpot, zgdept   ! 3D workspace
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ztsn                       ! 4D workspace
 
@@ -132,5 +132,9 @@
          ztsn(:,:,:,jp_tem) = ts(:,:,:,jp_tem,Kmm)                    ! thermosteric ssh
          ztsn(:,:,:,jp_sal) = sn0(:,:,:)
-         CALL eos( ztsn, zrhd, gdept(:,:,:,Kmm) )                       ! now in situ density using initial salinity
+!!st to be improved
+         DO jk = 1, jpkm1
+            zgdept(:,:,jk) = gdept(:,:,jk,Kmm)
+         END DO
+         CALL eos( ztsn, zrhd, zgdept)                       ! now in situ density using initial salinity
          !
          zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
@@ -159,5 +163,5 @@
       
          !                                         ! steric sea surface height
-         CALL eos( ts(:,:,:,:,Kmm), zrhd, zrhop, gdept(:,:,:,Kmm) )                 ! now in situ and potential density
+         CALL eos( ts(:,:,:,:,Kmm), zrhd, zrhop, zgdept )                 ! now in situ and potential density
          zrhop(:,:,jpk) = 0._wp
          CALL iom_put( 'rhop', zrhop )
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diawri.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diawri.F90	(revision 12724)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diawri.F90	(revision 12731)
@@ -435,5 +435,5 @@
       !
       REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d, ze3t       ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d, ze3t, zgdept       ! 3D workspace
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl   ! ABL 3D workspace
       !!----------------------------------------------------------------------
@@ -479,4 +479,5 @@
       DO jk = 1, jpk
          ze3t(:,:,jk) =  e3t(:,:,jk,Kmm)
+         zgdept(:,:,jk) =  gdept(:,:,jk,Kmm)
       END DO
 
@@ -805,5 +806,5 @@
          zw3d(:,:,:) = ( ( ze3t(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
          CALL histwrite( nid_T, "vovvle3t", it, ze3t (:,:,:)     , ndim_T , ndex_T  )   ! level thickness
-         CALL histwrite( nid_T, "vovvldep", it, gdept(:,:,:,Kmm) , ndim_T , ndex_T  )   ! t-point depth
+         CALL histwrite( nid_T, "vovvldep", it, zgdept , ndim_T , ndex_T  )   ! t-point depth !!st patch
          CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
       ENDIF
@@ -944,5 +945,5 @@
       !!
       INTEGER :: inum, jk
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t       ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, zgdept      ! 3D workspace
       !!----------------------------------------------------------------------
       ! 
@@ -954,6 +955,8 @@
       ENDIF 
       !
+!!st patch
       DO jk = 1, jpk
          ze3t(:,:,jk) =  e3t(:,:,jk,Kmm)
+         zgdept(:,:,jk) =  gdept(:,:,jk,Kmm)
       END DO
       !
@@ -962,7 +965,7 @@
       CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) )    ! now temperature
       CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) )    ! now salinity
-      CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm)              )    ! sea surface height
-      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm)                )    ! now i-velocity
-      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)                )    ! now j-velocity
+      CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:        ,Kmm) )    ! sea surface height
+      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:       ,Kmm) )    ! now i-velocity
+      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:       ,Kmm) )    ! now j-velocity
       IF( ln_zad_Aimp ) THEN
          CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww + wi        )    ! now k-velocity
@@ -971,5 +974,5 @@
       ENDIF
       CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep            )    ! now k-velocity
-      CALL iom_rstput( 0, 0, inum, 'ht'     , ht                 )    ! now water column height
+      CALL iom_rstput( 0, 0, inum, 'ht'     , ht(:,:)            )    ! now water column height
       !
       IF ( ln_isf ) THEN
@@ -1007,7 +1010,8 @@
       CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress
       CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
-      IF(  .NOT.ln_linssh  ) THEN             
-         CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept(:,:,:,Kmm)        )    !  T-cell depth 
-         CALL iom_rstput( 0, 0, inum, 'vovvle3t', ze3t(:,:,:)          )    !  T-cell thickness  
+      IF(  .NOT.ln_linssh  ) THEN
+!!st patch
+         CALL iom_rstput( 0, 0, inum, 'vovvldep', zgdept        )    !  T-cell depth 
+         CALL iom_rstput( 0, 0, inum, 'vovvle3t', ze3t          )    !  T-cell thickness  
       END IF
       IF( ln_wave .AND. ln_sdw ) THEN
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/dom_oce.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/dom_oce.F90	(revision 12724)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/dom_oce.F90	(revision 12731)
@@ -159,5 +159,7 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   hf_0, r1_hf_0   !: f-depth        [m] and [1/m]
    !                                                        ! time-dependent heights of ocean water column
+#if ! defined key_qco
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   ht          !: t-points           [m]
+#endif
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hu, r1_hu   !: u-depth            [m] and [1/m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hv, r1_hv   !: v-depth            [m] and [1/m]
@@ -293,7 +295,15 @@
       ii = ii+1
       ALLOCATE( ht_0(jpi,jpj) ,    hu_0(jpi,jpj)    ,    hv_0(jpi,jpj)     , hf_0(jpi,jpj) ,       &
-         &      ht  (jpi,jpj) ,    hu  (jpi,jpj,jpt),    hv  (jpi,jpj,jpt)                 ,       &
-         &                      r1_hu  (jpi,jpj,jpt), r1_hv  (jpi,jpj,jpt)                 ,       &
          &   r1_ht_0(jpi,jpj) , r1_hu_0(jpi,jpj) ,    r1_hv_0(jpi,jpj),   r1_hf_0(jpi,jpj) ,   STAT=ierr(ii)  )
+         !
+#if ! defined key_qco
+      ii = ii+1
+      ALLOCATE( ht  (jpi,jpj) ,    hu  (jpi,jpj,jpt),    hv  (jpi,jpj,jpt)                 ,       &
+         &                      r1_hu  (jpi,jpj,jpt), r1_hv  (jpi,jpj,jpt)                 ,   STAT=ierr(ii)  )
+#else
+      ii = ii+1
+      ALLOCATE(                    hu  (jpi,jpj,jpt),    hv  (jpi,jpj,jpt)                 ,       &
+         &                      r1_hu  (jpi,jpj,jpt), r1_hv  (jpi,jpj,jpt)                 ,   STAT=ierr(ii)  )
+#endif
          !
       ii = ii+1
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domain.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domain.F90	(revision 12724)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domain.F90	(revision 12731)
@@ -149,5 +149,5 @@
       CALL dom_msk( ik_top, ik_bot )    ! Masks
       !
-      ht_0(:,:) = 0._wp  ! Reference ocean thickness
+!!st      ht_0(:,:) = 0._wp  ! Reference ocean thickness
       hu_0(:,:) = 0._wp
       hv_0(:,:) = 0._wp
@@ -166,4 +166,12 @@
 
       !
+#if defined key_qco
+      !           !==  initialisation of time varying coordinate  ==!   Quasi-Euerian coordinate case
+      !
+      IF( .NOT.l_offline )   CALL dom_qe_init( Kbb, Kmm, Kaa )
+      !
+      IF( ln_linssh )        CALL ctl_stop('STOP','domain: key_qco and ln_linssh = T are incompatible')
+      !
+#else
       !           !==  time varying part of coordinate system  ==!
       !
@@ -176,7 +184,4 @@
             gde3w(:,:,:)    = gde3w_0(:,:,:)    ! = gdept as the sum of e3t
          !
-#if defined key_qco
-         !  Quasi-Euerian coordinate : no initialisation of e3. scale factors
-#else
          DO jt = 1, jpt                         ! vertical scale factors
             e3t(:,:,:,jt) =  e3t_0(:,:,:)
@@ -189,6 +194,4 @@
             e3f(:,:,:)    =  e3f_0(:,:,:)
          !
-#endif
-         !
          DO jt = 1, jpt                         ! water column thickness and its inverse
             hu(:,:,jt)    =    hu_0(:,:)
@@ -201,11 +204,9 @@
       ELSE                       != time varying : initialize before/now/after variables
          !
-#if defined key_qco
-         IF( .NOT.l_offline )   CALL dom_qe_init( Kbb, Kmm, Kaa )
-#else
          IF( .NOT.l_offline )   CALL dom_vvl_init( Kbb, Kmm, Kaa )
+         !
+      ENDIF
 #endif
-         !
-      ENDIF
+
       !
       IF( lk_c1d         )   CALL cor_c1d       ! 1D configuration: Coriolis set at T-point
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90	(revision 12724)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90	(revision 12731)
@@ -44,6 +44,6 @@
 !!st   PUBLIC  dom_qe_sf_nxt     ! called by steplf.F90
 !!st   PUBLIC  dom_qe_sf_update  ! called by steplf.F90
-   PUBLIC  dom_h_nxt         ! called by steplf.F90
-   PUBLIC  dom_h_update      ! called by steplf.F90
+!!st   PUBLIC  dom_h_nxt         ! called by steplf.F90
+!!st   PUBLIC  dom_h_update      ! called by steplf.F90
    PUBLIC  dom_qe_r3c        ! called by steplf.F90
 
@@ -175,50 +175,52 @@
 !!st end
       !
+!!st ATTENTION CAS ISF A GERER !!! 
       !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
-      IF( ln_isf ) THEN          !** IceShelF cavities
-      !                             ! to be created depending of the new names in isf
-      !                             ! it should be something like that :  (with h_isf = thickness of iceshelf)
-      !                             !  in fact currently, h_isf(:,:) is called : risfdep(:,:)
-!!gm - depth idea 0  : just realize that mask is not needed ===>>>> with ISF, rescale all grid point position below ISF : no mask !
-         gdept(:,:,1,Kmm) = gdept_0(:,:,1) * ( 1._wp + r3t(:,:,Kmm) )
-         gdepw(:,:,1,Kmm) = 0._wp                            ! Initialized to zero one for all
-         gde3w(:,:,1)     = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
-         DO jk = 2, jpk
-            gdept(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
-                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
-            gdepw(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
-                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
-            gde3w(:,:,jk)     = gdept(:,:,jk,Kmm) - ssh(:,:,Kmm)
-            gdept(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
-                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
-            gdepw(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
-                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
-         END DO
-         !
-      ELSE                       !** No cavities   (all depth rescaled, even inside topography: no mask)
-         !
-!!gm idea 0 :   just realize that mask is not needed ===>>>> without ISF, rescale all grid point position : no mask !
-         DO jk = 1, jpk
-            gdept(:,:,jk,Kmm) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
-            gdepw(:,:,jk,Kmm) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
-            gde3w(:,:,jk)     = gdept  (:,:,jk,Kmm)       - ssh(:,:,Kmm)
-            gdept(:,:,jk,Kbb) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
-            gdepw(:,:,jk,Kbb) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
-         END DO
-         !
-      ENDIF
+!!$      IF( ln_isf ) THEN          !** IceShelF cavities
+!!$      !                             ! to be created depending of the new names in isf
+!!$      !                             ! it should be something like that :  (with h_isf = thickness of iceshelf)
+!!$      !                             !  in fact currently, h_isf(:,:) is called : risfdep(:,:)
+!!$!!gm - depth idea 0  : just realize that mask is not needed ===>>>> with ISF, rescale all grid point position below ISF : no mask !
+!!$         gdept(:,:,1,Kmm) = gdept_0(:,:,1) * ( 1._wp + r3t(:,:,Kmm) )
+!!$         gdepw(:,:,1,Kmm) = 0._wp                            ! Initialized to zero one for all
+!!$         gde3w(:,:,1)     = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
+!!$         DO jk = 2, jpk
+!!$            gdept(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
+!!$                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
+!!$            gdepw(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
+!!$                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
+!!$            gde3w(:,:,jk)     = gdept(:,:,jk,Kmm) - ssh(:,:,Kmm)
+!!$            gdept(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
+!!$                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
+!!$            gdepw(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
+!!$                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
+!!$         END DO
+!!$         !
+!!$      ELSE                       !** No cavities   (all depth rescaled, even inside topography: no mask)
+!!$         !
+!!$!!gm idea 0 :   just realize that mask is not needed ===>>>> without ISF, rescale all grid point position : no mask !
+!!$         DO jk = 1, jpk
+!!$            gdept(:,:,jk,Kmm) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
+!!$            gdepw(:,:,jk,Kmm) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
+!!$            gde3w(:,:,jk)     = gdept  (:,:,jk,Kmm)       - ssh(:,:,Kmm)
+!!$            gdept(:,:,jk,Kbb) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
+!!$            gdepw(:,:,jk,Kbb) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
+!!$         END DO
+!!$         !
+!!$      ENDIF
       !
       !                    !==  thickness of the water column  !!   (ocean portion only)
-      ht(:,:)     = ht_0(:,:)           + ssh(:,:,Kmm)
-      hu(:,:,Kbb) = hu_0(:,:) * ( 1._wp + r3u(:,:,Kbb) )
-      hu(:,:,Kmm) = hu_0(:,:) * ( 1._wp + r3u(:,:,Kmm) )
-      hv(:,:,Kbb) = hv_0(:,:) * ( 1._wp + r3v(:,:,Kbb) )
-      hv(:,:,Kmm) = hv_0(:,:) * ( 1._wp + r3v(:,:,Kmm) )
-      !
-      !                    !==  inverse of water column thickness   ==!   (u- and v- points)
-      r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
-      r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) )
-      r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) )
-      r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) )
+!!st ht(:,:)     = ht_0(:,:)           + ssh(:,:,Kmm)
+!!$      hu(:,:,Kbb) = hu_0(:,:) * ( 1._wp + r3u(:,:,Kbb) )
+!!$      hu(:,:,Kmm) = hu_0(:,:) * ( 1._wp + r3u(:,:,Kmm) )
+!!$      hv(:,:,Kbb) = hv_0(:,:) * ( 1._wp + r3v(:,:,Kbb) )
+!!$      hv(:,:,Kmm) = hv_0(:,:) * ( 1._wp + r3v(:,:,Kmm) )
+!!$      !
+!!$      !                    !==  inverse of water column thickness   ==!   (u- and v- points)
+!!$      r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
+!!$      r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) )
+!!$      r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) )
+!!$      r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) )
+!!st end
       !
    END SUBROUTINE dom_qe_zgr
@@ -294,46 +296,46 @@
 !    END SUBROUTINE dom_qe_sf_nxt
 !!st end
-
-   SUBROUTINE dom_h_nxt( kt, Kbb, Kmm, Kaa, kcall )
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE dom_qe_sf_nxt  ***
-      !!
-      !! ** Purpose :  - compute the after water heigh used in tra_zdf, dynnxt,
-      !!                 tranxt and dynspg routines
-      !!
-      !! ** Method  :  - z_star case:  Proportionnaly to the water column thickness.
-      !!
-      !! ** Action  :  - h(u/v) update wrt ssh/h(u/v)_0
-      !!
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT( in )           ::   kt             ! time step
-      INTEGER, INTENT( in )           ::   Kbb, Kmm, Kaa  ! time step
-      INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence
-      !
-      !!----------------------------------------------------------------------
-      !
-      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
-      !
-      IF( ln_timing )   CALL timing_start('dom_h_nxt')
-      !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dom_h_nxt : compute after scale factors'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
-      ENDIF
-      !
-      ! *********************************** !
-      ! After depths at u- v points         !
-      ! *********************************** !
-      hu(:,:,Kaa) = hu_0(:,:) * ( 1._wp + r3u(:,:,Kaa) )
-      hv(:,:,Kaa) = hv_0(:,:) * ( 1._wp + r3v(:,:,Kaa) )
-      !                                        ! Inverse of the local depth
-      r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
-      r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
-      !
-      IF( ln_timing )   CALL timing_stop('dom_h_nxt')
-      !
-   END SUBROUTINE dom_h_nxt
-
+!!st 
+!!$   SUBROUTINE dom_h_nxt( kt, Kbb, Kmm, Kaa, kcall )
+!!$      !!----------------------------------------------------------------------
+!!$      !!                ***  ROUTINE dom_qe_sf_nxt  ***
+!!$      !!
+!!$      !! ** Purpose :  - compute the after water heigh used in tra_zdf, dynnxt,
+!!$      !!                 tranxt and dynspg routines
+!!$      !!
+!!$      !! ** Method  :  - z_star case:  Proportionnaly to the water column thickness.
+!!$      !!
+!!$      !! ** Action  :  - h(u/v) update wrt ssh/h(u/v)_0
+!!$      !!
+!!$      !!----------------------------------------------------------------------
+!!$      INTEGER, INTENT( in )           ::   kt             ! time step
+!!$      INTEGER, INTENT( in )           ::   Kbb, Kmm, Kaa  ! time step
+!!$      INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence
+!!$      !
+!!$      !!----------------------------------------------------------------------
+!!$      !
+!!$      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
+!!$      !
+!!$      IF( ln_timing )   CALL timing_start('dom_h_nxt')
+!!$      !
+!!$      IF( kt == nit000 ) THEN
+!!$         IF(lwp) WRITE(numout,*)
+!!$         IF(lwp) WRITE(numout,*) 'dom_h_nxt : compute after scale factors'
+!!$         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
+!!$      ENDIF
+!!$      !
+!!$      ! *********************************** !
+!!$      ! After depths at u- v points         !
+!!$      ! *********************************** !
+!!$      hu(:,:,Kaa) = hu_0(:,:) * ( 1._wp + r3u(:,:,Kaa) )
+!!$      hv(:,:,Kaa) = hv_0(:,:) * ( 1._wp + r3v(:,:,Kaa) )
+!!$      !                                        ! Inverse of the local depth
+!!$      r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
+!!$      r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
+!!$      !
+!!$      IF( ln_timing )   CALL timing_stop('dom_h_nxt')
+!!$      !
+!!$   END SUBROUTINE dom_h_nxt
+!!st end
 ! !!st
 !    SUBROUTINE dom_qe_sf_update( kt, Kbb, Kmm, Kaa )
@@ -445,89 +447,89 @@
 !!st end
 
-   SUBROUTINE dom_h_update( kt, Kbb, Kmm, Kaa )
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE dom_qe_sf_update  ***
-      !!
-      !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors
-      !!               compute all depths and related variables for next time step
-      !!               write outputs and restart file
-      !!
-      !! ** Method  :  - reconstruct scale factor at other grid points (interpolate)
-      !!               - recompute depths and water height fields
-      !!
-      !! ** Action  :  - Recompute:
-      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w
-      !!                    h(u/v) and h(u/v)r
-      !!
-      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
-      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt              ! time step
-      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices
-      !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zcoef        ! local scalar
-      !!----------------------------------------------------------------------
-      !
-      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
-      !
-      IF( ln_timing )   CALL timing_start('dom_qe_sf_update')
-      !
-      IF( kt == nit000 )   THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dom_qe_sf_update : - interpolate scale factors and compute depths for next time step'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
-      ENDIF
-      !
-      ! Compute all missing vertical scale factor and depths
-      ! ====================================================
-      ! Horizontal scale factor interpolations
-      ! --------------------------------------
-      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
-
-      IF( ln_isf ) THEN          !** IceShelF cavities
-      !                             ! to be created depending of the new names in isf
-      !                             ! it should be something like that :  (with h_isf = thickness of iceshelf)
-      !                             !  in fact currently, h_isf(:,:) is called : risfdep(:,:)
-   !!gm - depth idea 0  : just realize that mask is not needed ===>>>> with ISF, rescale all grid point position below ISF : no mask !
-         gdept(:,:,1,Kmm) = gdept_0(:,:,1) * ( 1._wp + r3t(:,:,Kmm) )
-         gdepw(:,:,1,Kmm) = 0._wp                            ! Initialized to zero one for all
-         gde3w(:,:,1)     = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
-         DO jk = 2, jpk
-            gdept(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
-                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
-            gdepw(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
-                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
-            gde3w(:,:,jk)     = gdept(:,:,jk,Kmm) - ssh(:,:,Kmm)
-            gdept(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
-                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
-            gdepw(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
-                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
-         END DO
-         !
-      ELSE                       !** No cavities   (all depth rescaled, even inside topography: no mask)
-         !
-   !!gm idea 0 :   just realize that mask is not needed ===>>>> without ISF, rescale all grid point position : no mask !
-         DO jk = 1, jpk
-            gdept(:,:,jk,Kmm) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
-            gdepw(:,:,jk,Kmm) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
-            gde3w(:,:,jk)     = gdept  (:,:,jk,Kmm)       - ssh(:,:,Kmm)
-            gdept(:,:,jk,Kbb) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
-            gdepw(:,:,jk,Kbb) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
-         END DO
-         !
-      ENDIF
-
-      ! Local depth and Inverse of the local depth of the water
-      ! -------------------------------------------------------
-      !
-      ht(:,:) = ht_0(:,:) + ssh(:,:,Kmm)
-
-      ! write restart file
-      ! ==================
-      IF( ln_timing )   CALL timing_stop('dom_qe_sf_update')
-      !
-   END SUBROUTINE dom_h_update
-
+!!$   SUBROUTINE dom_h_update( kt, Kbb, Kmm, Kaa )
+!!$      !!----------------------------------------------------------------------
+!!$      !!                ***  ROUTINE dom_qe_sf_update  ***
+!!$      !!
+!!$      !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors
+!!$      !!               compute all depths and related variables for next time step
+!!$      !!               write outputs and restart file
+!!$      !!
+!!$      !! ** Method  :  - reconstruct scale factor at other grid points (interpolate)
+!!$      !!               - recompute depths and water height fields
+!!$      !!
+!!$      !! ** Action  :  - Recompute:
+!!$      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w
+!!$      !!                    h(u/v) and h(u/v)r
+!!$      !!
+!!$      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
+!!$      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
+!!$      !!----------------------------------------------------------------------
+!!$      INTEGER, INTENT( in ) ::   kt              ! time step
+!!$      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices
+!!$      !
+!!$      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+!!$      REAL(wp) ::   zcoef        ! local scalar
+!!$      !!----------------------------------------------------------------------
+!!$      !
+!!$      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
+!!$      !
+!!$      IF( ln_timing )   CALL timing_start('dom_qe_sf_update')
+!!$      !
+!!$      IF( kt == nit000 )   THEN
+!!$         IF(lwp) WRITE(numout,*)
+!!$         IF(lwp) WRITE(numout,*) 'dom_qe_sf_update : - interpolate scale factors and compute depths for next time step'
+!!$         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
+!!$      ENDIF
+!!$      !
+!!$      ! Compute all missing vertical scale factor and depths
+!!$      ! ====================================================
+!!$      ! Horizontal scale factor interpolations
+!!$      ! --------------------------------------
+!!$      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
+!!$
+!!$      IF( ln_isf ) THEN          !** IceShelF cavities
+!!$      !                             ! to be created depending of the new names in isf
+!!$      !                             ! it should be something like that :  (with h_isf = thickness of iceshelf)
+!!$      !                             !  in fact currently, h_isf(:,:) is called : risfdep(:,:)
+!!$   !!gm - depth idea 0  : just realize that mask is not needed ===>>>> with ISF, rescale all grid point position below ISF : no mask !
+!!$         gdept(:,:,1,Kmm) = gdept_0(:,:,1) * ( 1._wp + r3t(:,:,Kmm) )
+!!$         gdepw(:,:,1,Kmm) = 0._wp                            ! Initialized to zero one for all
+!!$         gde3w(:,:,1)     = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
+!!$         DO jk = 2, jpk
+!!$            gdept(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
+!!$                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
+!!$            gdepw(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
+!!$                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
+!!$            gde3w(:,:,jk)     = gdept(:,:,jk,Kmm) - ssh(:,:,Kmm)
+!!$            gdept(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
+!!$                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
+!!$            gdepw(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
+!!$                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
+!!$         END DO
+!!$         !
+!!$      ELSE                       !** No cavities   (all depth rescaled, even inside topography: no mask)
+!!$         !
+!!$   !!gm idea 0 :   just realize that mask is not needed ===>>>> without ISF, rescale all grid point position : no mask !
+!!$         DO jk = 1, jpk
+!!$            gdept(:,:,jk,Kmm) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
+!!$            gdepw(:,:,jk,Kmm) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
+!!$            gde3w(:,:,jk)     = gdept  (:,:,jk,Kmm)       - ssh(:,:,Kmm)
+!!$            gdept(:,:,jk,Kbb) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
+!!$            gdepw(:,:,jk,Kbb) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
+!!$         END DO
+!!$         !
+!!$      ENDIF
+!!$
+!!$      ! Local depth and Inverse of the local depth of the water
+!!$      ! -------------------------------------------------------
+!!$      !
+!!$!!st      ht(:,:) = ht_0(:,:) + ssh(:,:,Kmm)
+!!$
+!!$      ! write restart file
+!!$      ! ==================
+!!$      IF( ln_timing )   CALL timing_stop('dom_qe_sf_update')
+!!$      !
+!!$   END SUBROUTINE dom_h_update
+!!st end
 
    SUBROUTINE dom_qe_r3c( pssh, pr3t, pr3u, pr3v, pr3f )
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domzgr_substitute.h90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domzgr_substitute.h90	(revision 12724)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domzgr_substitute.h90	(revision 12731)
@@ -19,4 +19,12 @@
 #   define  e3uw(i,j,k,t)  (e3uw_0(i,j,k)*(1._wp+r3u(i,j,t)))
 #   define  e3vw(i,j,k,t)  (e3vw_0(i,j,k)*(1._wp+r3v(i,j,t)))
+#   define  ht(i,j)        (ht_0(i,j)+ssh(i,j,Kmm))
+#   define  hu(i,j,t)      (hu_0(i,j)*(1._wp+r3u(i,j,t)))
+#   define  hv(i,j,t)      (hv_0(i,j)*(1._wp+r3v(i,j,t)))
+#   define  r1_hu(i,j,t)   (r1_hu_0(i,j)/(1._wp+r3u(i,j,t)))
+#   define  r1_hv(i,j,t)   (r1_hv_0(i,j)/(1._wp+r3v(i,j,t)))
+#   define  gdept(i,j,k,t) (gdept_0(i,j,k)*(1._wp+r3t(i,j,t)))
+#   define  gdepw(i,j,k,t) (gdepw_0(i,j,k)*(1._wp+r3t(i,j,t)))
+#   define  gde3w(i,j,k)   (gdept_0(i,j,k)*(1._wp+r3t(i,j,Kmm))-ssh(i,j,Kmm))
 #endif
 !!----------------------------------------------------------------------
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/istate.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/istate.F90	(revision 12724)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/istate.F90	(revision 12731)
@@ -60,4 +60,5 @@
       !
       INTEGER ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zgdept     !!st temporary modif to be able to use gdept subtitute
 !!gm see comment further down
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace
@@ -116,5 +117,9 @@
             !
          ELSE                                 ! user defined initial T and S
-            CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )         
+!!st zgdept 
+            DO jk = 1, jpkm1
+               zgdept(:,:,jk) = gdept(:,:,jk,Kbb)
+            END DO
+            CALL usr_def_istate( zgdept, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )         
          ENDIF
          ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
@@ -128,5 +133,5 @@
 !!gm POTENTIAL BUG :
 !!gm  ISSUE :  if ssh(:,:,Kbb) /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
-!!             as well as gdept and gdepw....   !!!!! 
+!!             as well as gdept_ and gdepw_....   !!!!! 
 !!      ===>>>>   probably a call to domvvl initialisation here....
 
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatf.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatf.F90	(revision 12724)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatf.F90	(revision 12731)
@@ -215,8 +215,14 @@
             zwfld(:,:) = emp_b(:,:) - emp(:,:)
             IF ( ln_rnf ) zwfld(:,:) =  zwfld(:,:) - ( rnf_b(:,:) - rnf(:,:) )
-            DO jk = 1, jpkm1
-               ze3t_f(:,:,jk) = ze3t_f(:,:,jk) - zcoef * zwfld(:,:) * tmask(:,:,jk) &
-                              &                        * pe3t(:,:,jk,Kmm) / ( ht(:,:) + 1._wp - ssmask(:,:) )
-            END DO
+!!st
+!!$            DO jk = 1, jpkm1
+!!$               ze3t_f(:,:,jk) = ze3t_f(:,:,jk) - zcoef * zwfld(:,:) * tmask(:,:,jk) &
+!!$                              &                        * pe3t(:,:,jk,Kmm) / ( ht_(:,:) + 1._wp - ssmask(:,:) )
+!!$            END DO
+!!st end            
+            DO_3D_11_11( 1, jpkm1 )
+               ze3t_f(ji,jj,jk) = ze3t_f(ji,jj,jk) - zcoef * zwfld(ji,jj) * tmask(ji,jj,jk) &
+                              &                * pe3t(ji,jj,jk,Kmm) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )
+            END_3D
             !
             ! ice shelf melting (deal separately as it can be in depth)
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatfQCO.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatfQCO.F90	(revision 12724)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatfQCO.F90	(revision 12731)
@@ -221,14 +221,15 @@
       ! integration
       !
-      IF(.NOT.ln_linssh ) THEN
-         hu(:,:,Kmm) = e3u(:,:,1,Kmm ) * umask(:,:,1)
-         hv(:,:,Kmm) = e3v(:,:,1,Kmm ) * vmask(:,:,1)
-         DO jk = 2, jpkm1
-            hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm ) * umask(:,:,jk)
-            hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm ) * vmask(:,:,jk)
-         END DO
-         r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) )
-         r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) )
-      ENDIF
+!!st
+!!$      IF(.NOT.ln_linssh ) THEN
+!!$         hu(:,:,Kmm) = e3u(:,:,1,Kmm ) * umask(:,:,1)
+!!$         hv(:,:,Kmm) = e3v(:,:,1,Kmm ) * vmask(:,:,1)
+!!$         DO jk = 2, jpkm1
+!!$            hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm ) * umask(:,:,jk)
+!!$            hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm ) * vmask(:,:,jk)
+!!$         END DO
+!!$         r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) )
+!!$         r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) )
+!!$      ENDIF
       !
       uu_b(:,:,Kaa) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynhpg.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynhpg.F90	(revision 12724)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynhpg.F90	(revision 12731)
@@ -781,6 +781,6 @@
       !-------------------------------------------------------------
 
-!!bug gm   :  e3w-gde3w = 0.5*e3w  ....  and gde3w(2)-gde3w(1)=e3w(:,:,2,:) ....   to be verified
-!          true if gde3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be
+!!bug gm   :  e3w-gde3w(:,:,:) = 0.5*e3w  ....  and gde3w(:,:,2)-gde3w(:,:,1)=e3w(:,:,2,Kmm) ....   to be verified
+!          true if gde3w(:,:,:) is really defined as the sum of the e3w scale factors as, it seems to me, it should be
 
       DO_2D_00_00
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynspg_ts.F90	(revision 12724)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynspg_ts.F90	(revision 12731)
@@ -257,5 +257,5 @@
       zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:)        ! NB: FULL domain : put a value in last row and column
       !
-      CALL dyn_cor_2d( hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV,  &   ! <<== in
+      CALL dyn_cor_2d( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV,  &   ! <<== in
          &                                                                     zu_trd, zv_trd   )   ! ==>> out
       !
@@ -388,5 +388,5 @@
          zhup2_e(:,:) = hu(:,:,Kmm)
          zhvp2_e(:,:) = hv(:,:,Kmm)
-         zhtp2_e(:,:) = ht(:,:)
+!!st not used ?         zhtp2_e(:,:) = ht_(:,:)
       ENDIF
       !
@@ -574,5 +574,5 @@
          ! at each time step. We however keep them constant here for optimization.
          ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated)
-         CALL dyn_cor_2d( zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV,    zu_trd, zv_trd   )
+         CALL dyn_cor_2d( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV,    zu_trd, zv_trd   )
          !
          ! Add tidal astronomical forcing if defined
@@ -1186,5 +1186,5 @@
 
 
-   SUBROUTINE dyn_cor_2d( phu, phv, punb, pvnb, zhU, zhV,    zu_trd, zv_trd   )
+   SUBROUTINE dyn_cor_2d( pht, phu, phv, punb, pvnb, zhU, zhV,    zu_trd, zv_trd   )
       !!---------------------------------------------------------------------
       !!                   ***  ROUTINE dyn_cor_2d  ***
@@ -1194,5 +1194,5 @@
       INTEGER  ::   ji ,jj                             ! dummy loop indices
       REAL(wp) ::   zx1, zx2, zy1, zy2, z1_hu, z1_hv   !   -      -
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: phu, phv, punb, pvnb, zhU, zhV
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pht, phu, phv, punb, pvnb, zhU, zhV
       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: zu_trd, zv_trd
       !!----------------------------------------------------------------------
@@ -1203,10 +1203,10 @@
             z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) )
             zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu                    &
-               &               * (  e1e2t(ji+1,jj)*ht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) )   &
-               &                  + e1e2t(ji  ,jj)*ht(ji  ,jj)*ff_t(ji  ,jj) * ( pvnb(ji  ,jj) + pvnb(ji  ,jj-1) )   )
+               &               * (  e1e2t(ji+1,jj)*pht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) )   &
+               &                  + e1e2t(ji  ,jj)*pht(ji  ,jj)*ff_t(ji  ,jj) * ( pvnb(ji  ,jj) + pvnb(ji  ,jj-1) )   )
                !
             zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    &
-               &               * (  e1e2t(ji,jj+1)*ht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) )   & 
-               &                  + e1e2t(ji,jj  )*ht(ji,jj  )*ff_t(ji,jj  ) * ( punb(ji,jj  ) + punb(ji-1,jj  ) )   ) 
+               &               * (  e1e2t(ji,jj+1)*pht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) )   & 
+               &                  + e1e2t(ji,jj  )*pht(ji,jj  )*ff_t(ji,jj  ) * ( punb(ji,jj  ) + punb(ji-1,jj  ) )   ) 
          END_2D
          !         
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfcpl.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfcpl.F90	(revision 12724)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfcpl.F90	(revision 12731)
@@ -142,11 +142,13 @@
       !!----------------------------------------------------------------------
       INTEGER :: jk                               ! loop index
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, ze3u, ze3v    ! e3t , e3u, e3v
-      !!----------------------------------------------------------------------
-      !
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, ze3u, ze3v, zgdepw  ! e3t , e3u, e3v
+      !!----------------------------------------------------------------------
+      !
+!!st patch to be able to use substitution
       DO jk = 1, jpk
          ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
          ze3u(:,:,jk) = e3u(:,:,jk,Kmm)
          ze3v(:,:,jk) = e3v(:,:,jk,Kmm)
+         zgdepw(:,:,jk) = gdepw(:,:,jk,Kmm)
       END DO 
       !
@@ -157,5 +159,5 @@
       CALL iom_rstput( kt, nitrst, numrow, 'e3u_n'  , ze3u , ldxios = lwxios )
       CALL iom_rstput( kt, nitrst, numrow, 'e3v_n'  , ze3v , ldxios = lwxios )
-      CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', gdepw(:,:,:,Kmm) , ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', zgdepw , ldxios = lwxios )
       IF( lwxios ) CALL iom_swap( cxios_context )
       !
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfdynatf.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfdynatf.F90	(revision 12724)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfdynatf.F90	(revision 12731)
@@ -15,5 +15,6 @@
    USE phycst , ONLY: r1_rho0         ! physical constant
    USE dom_oce                        ! time and space domain
-
+   USE oce, ONLY : ssh                ! sea-surface height
+!!st
    USE in_out_manager
 
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfstp.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfstp.F90	(revision 12724)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfstp.F90	(revision 12731)
@@ -11,17 +11,17 @@
 
    !!----------------------------------------------------------------------
-   !!   isfstp       : compute iceshelf melt and heat flux
-   !!----------------------------------------------------------------------
-   !
-   USE isf_oce                                      ! isf variables
-   USE isfload, ONLY: isf_load                      ! ice shelf load
-   USE isftbl , ONLY: isf_tbl_lvl                   ! ice shelf boundary layer
-   USE isfpar , ONLY: isf_par, isf_par_init         ! ice shelf parametrisation
-   USE isfcav , ONLY: isf_cav, isf_cav_init         ! ice shelf cavity
-   USE isfcpl , ONLY: isfcpl_rst_write, isfcpl_init ! isf variables
-
-   USE dom_oce                                          ! ocean space and time domain
-   USE domvvl,  ONLY: ln_vvl_zstar                      ! zstar logical
-   USE zdfdrg,  ONLY: r_Cdmin_top, r_ke0_top            ! vertical physics: top/bottom drag coef.
+   !!   isfstp        : compute iceshelf melt and heat flux
+   !!----------------------------------------------------------------------
+   USE isf_oce        ! isf variables
+   USE isfload  , ONLY: isf_load                      ! ice shelf load
+   USE isftbl   , ONLY: isf_tbl_lvl                   ! ice shelf boundary layer
+   USE isfpar   , ONLY: isf_par, isf_par_init         ! ice shelf parametrisation
+   USE isfcav   , ONLY: isf_cav, isf_cav_init         ! ice shelf cavity
+   USE isfcpl   , ONLY: isfcpl_rst_write, isfcpl_init ! isf variables
+
+   USE dom_oce        ! ocean space and time domain
+   USE oce      , ONLY: ssh                           ! sea surface height
+   USE domvvl   , ONLY: ln_vvl_zstar                  ! zstar logical
+   USE zdfdrg   , ONLY: r_Cdmin_top, r_ke0_top        ! vertical physics: top/bottom drag coef.
    !
    USE lib_mpp, ONLY: ctl_stop, ctl_nam
@@ -31,5 +31,4 @@
 
    IMPLICIT NONE
-
    PRIVATE
 
@@ -87,5 +86,5 @@
             ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
          END DO 
-         CALL isf_tbl_lvl(ht, ze3t, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav)
+         CALL isf_tbl_lvl(ht(:,:), ze3t, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav)
          !
          ! 1.3: compute ice shelf melt
@@ -112,5 +111,5 @@
             ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
          END DO
-         CALL isf_tbl_lvl(ht, ze3t, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par)
+         CALL isf_tbl_lvl(ht(:,:), ze3t, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par)
          !
          ! 2.3: compute ice shelf melt
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/LDF/ldfslp.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/LDF/ldfslp.F90	(revision 12724)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/LDF/ldfslp.F90	(revision 12731)
@@ -298,5 +298,5 @@
 !               !                                         ! jk must be >= ML level for zfk=1. otherwise  zfk=0.
 !               zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp )
-!               zck = gdepw(ji,jj,jk)    / MAX( hmlp(ji,jj), 10. )
+!               zck = gdepw(ji,jj,jk,Kmm)    / MAX( hmlp(ji,jj), 10. )
 !               zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk)
 !               zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk)
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdfosm.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdfosm.F90	(revision 12724)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdfosm.F90	(revision 12731)
@@ -475,5 +475,8 @@
 
       zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need wb here, so subtract it
-      zhbl_t(:,:) = MIN(zhbl_t(:,:), ht(:,:))
+!!st      zhbl_t(:,:) = MIN(zhbl_t(:,:), ht_(:,:))
+      DO_2D_11_11
+         zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), ht(ji,jj))
+      END_2D
       zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_Dt + ww(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom
 
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/stepLF.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/stepLF.F90	(revision 12724)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/stepLF.F90	(revision 12731)
@@ -90,4 +90,6 @@
 !!gm kcall can be removed, I guess
       INTEGER ::   kcall        ! optional integer argument (dom_vvl_sf_nxt)
+!!st patch
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgdept  ! st patch
       !! ---------------------------------------------------------------------
 #if defined key_agrif
@@ -194,12 +196,17 @@
       !  Ocean dynamics : hdiv, ssh, e3, u, v, w
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
+!!st patch
+!!st patch to be able to use substitution
+      DO jk = 1, jpk
+         zgdept(:,:,jk) = gdept(:,:,jk,Nnn)
+      END DO
+!!st end
                             CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh,  Naa )    ! after ssh (includes call to div_hor)
       IF( .NOT.ln_linssh )  CALL dom_qe_r3c    ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa) )
-      IF( .NOT.ln_linssh )  CALL dom_h_nxt     ( kstp, Nbb, Nnn,       Naa )    ! after vertical scale factors
+!!st      IF( .NOT.ln_linssh )  CALL dom_h_nxt     ( kstp, Nbb, Nnn,       Naa )    ! after vertical scale factors
       !IF( .NOT.ln_linssh )  CALL dom_qe_sf_nxt ( kstp, Nbb, Nnn,      Naa )    ! after vertical scale factors
                             CALL wzv           ( kstp, Nbb, Nnn, Naa, ww  )    ! Nnn cross-level velocity
       IF( ln_zad_Aimp )     CALL wAimp         ( kstp,      Nnn           )  ! Adaptive-implicit vertical advection partitioning
-                            CALL eos    ( ts(:,:,:,:,Nnn), rhd, rhop, gdept(:,:,:,Nnn) )  ! now in situ density for hpg computation
+                            CALL eos    ( ts(:,:,:,:,Nnn), rhd, rhop, zgdept )  ! now in situ density for hpg computation
 
 
@@ -226,5 +233,5 @@
          IF(.NOT.ln_linssh) CALL dom_qe_r3c ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )
          !IF(.NOT.ln_linssh) CALL dom_qe_sf_nxt ( kstp, Nbb, Nnn, Naa, kcall=2 )  ! after vertical scale factors (update depth average component)
-         IF(.NOT.ln_linssh) CALL dom_h_nxt  ( kstp, Nbb, Nnn, Naa, kcall=2 )  ! after vertical scale factors (update depth average component)
+!!st         IF(.NOT.ln_linssh) CALL dom_h_nxt  ( kstp, Nbb, Nnn, Naa, kcall=2 )  ! after vertical scale factors (update depth average component)
       ENDIF
                             CALL dyn_zdf    ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion
@@ -322,5 +329,5 @@
       !
       !IF(.NOT.ln_linssh) CALL dom_qe_sf_update( kstp, Nbb, Nnn, Naa )  ! recompute vertical scale factors
-      IF(.NOT.ln_linssh) CALL dom_h_update  ( kstp, Nbb, Nnn, Naa )  ! recompute vertical scale factors
+!!st      IF(.NOT.ln_linssh) CALL dom_h_update  ( kstp, Nbb, Nnn, Naa )  ! recompute vertical scale factors
       !
       IF( ln_diahsb  )   CALL dia_hsb       ( kstp, Nbb, Nnn )  ! - ML - global conservation diagnostics
