Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/iceistate.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/iceistate.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/iceistate.F90	(revision 12680)
@@ -60,5 +60,5 @@
    INTEGER , PARAMETER ::   jp_hpd = 9           ! index of pnd depth        (m)
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read)
-   !
+
    !! * Substitutions
 #  include "do_loop_substitute.h90"
@@ -102,5 +102,5 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
       REAL(wp), DIMENSION(jpi,jpj)     ::   zapnd_ini, zhpnd_ini                       !data from namelist or nc file
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zti_3d , zts_3d                            !temporary arrays
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zti_3d , zts_3d                            !locak arrays
       !!
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d
@@ -452,6 +452,6 @@
       !!
       !!-----------------------------------------------------------------------------
-      INTEGER ::   ios   ! Local integer output status for namelist read
-      INTEGER ::   ifpr, ierror
+      INTEGER ::   ios, ifpr, ierror   ! Local integers
+
       !
       CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files
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 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ASM/asminc.F90	(revision 12680)
@@ -9,5 +9,5 @@
    !!                 ! 2007-04  (A. Weaver)  Merge with OPAVAR/NEMOVAR
    !!   NEMO     3.3  ! 2010-05  (D. Lea)  Update to work with NEMO v3.2
-   !!             -   ! 2010-05  (D. Lea)  add calc_month_len routine based on day_init 
+   !!             -   ! 2010-05  (D. Lea)  add calc_month_len routine based on day_init
    !!            3.4  ! 2012-10  (A. Weaver and K. Mogensen) Fix for direct initialization
    !!                 ! 2014-09  (D. Lea)  Local calc_date removed use routine from OBS
@@ -31,5 +31,5 @@
    USE zpshde          ! Partial step : Horizontal Derivative
    USE asmpar          ! Parameters for the assmilation interface
-   USE asmbkg          ! 
+   USE asmbkg          !
    USE c1d             ! 1D initialization
    USE sbc_oce         ! Surface boundary condition variables.
@@ -45,5 +45,5 @@
    IMPLICIT NONE
    PRIVATE
-   
+
    PUBLIC   asm_inc_init   !: Initialize the increment arrays and IAU weights
    PUBLIC   tra_asm_inc    !: Apply the tracer (T and S) increments
@@ -72,5 +72,5 @@
    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkg   , v_bkg      !: Background u- & v- velocity components
    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkginc, s_bkginc   !: Increment to the background T & S
-   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkginc, v_bkginc   !: Increment to the u- & v-components 
+   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkginc, v_bkginc   !: Increment to the u- & v-components
    REAL(wp), PUBLIC, DIMENSION(:)    , ALLOCATABLE ::   wgtiau               !: IAU weights for each time step
 #if defined key_asminc
@@ -80,9 +80,9 @@
    INTEGER , PUBLIC ::   nitbkg      !: Time step of the background state used in the Jb term
    INTEGER , PUBLIC ::   nitdin      !: Time step of the background state for direct initialization
-   INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval 
+   INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval
    INTEGER , PUBLIC ::   nitiaufin   !: Time step of the end of the IAU interval
-   ! 
+   !
    INTEGER , PUBLIC ::   niaufn      !: Type of IAU weighing function: = 0   Constant weighting
-   !                                 !: = 1   Linear hat-like, centred in middle of IAU interval 
+   !                                 !: = 1   Linear hat-like, centred in middle of IAU interval
    REAL(wp), PUBLIC ::   salfixmin   !: Ensure that the salinity is larger than this  value if (ln_salfix)
 
@@ -106,10 +106,10 @@
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE asm_inc_init  ***
-      !!          
+      !!
       !! ** Purpose : Initialize the assimilation increment and IAU weights.
       !!
       !! ** Method  : Initialize the assimilation increment and IAU weights.
       !!
-      !! ** Action  : 
+      !! ** Action  :
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::  Kbb, Kmm, Krhs  ! time level indices
@@ -263,5 +263,5 @@
          !
          !                                !---------------------------------------------------------
-         IF( niaufn == 0 ) THEN           ! Constant IAU forcing 
+         IF( niaufn == 0 ) THEN           ! Constant IAU forcing
             !                             !---------------------------------------------------------
             DO jt = 1, iiauper
@@ -269,10 +269,10 @@
             END DO
             !                             !---------------------------------------------------------
-         ELSEIF ( niaufn == 1 ) THEN      ! Linear hat-like, centred in middle of IAU interval 
+         ELSEIF ( niaufn == 1 ) THEN      ! Linear hat-like, centred in middle of IAU interval
             !                             !---------------------------------------------------------
             ! Compute the normalization factor
             znorm = 0._wp
             IF( MOD( iiauper, 2 ) == 0 ) THEN   ! Even number of time steps in IAU interval
-               imid = iiauper / 2 
+               imid = iiauper / 2
                DO jt = 1, imid
                   znorm = znorm + REAL( jt )
@@ -280,5 +280,5 @@
                znorm = 2.0 * znorm
             ELSE                                ! Odd number of time steps in IAU interval
-               imid = ( iiauper + 1 ) / 2        
+               imid = ( iiauper + 1 ) / 2
                DO jt = 1, imid - 1
                   znorm = znorm + REAL( jt )
@@ -307,11 +307,11 @@
              DO jt = 1, icycper
                 ztotwgt = ztotwgt + wgtiau(jt)
-                WRITE(numout,*) '         ', jt, '       ', wgtiau(jt) 
-             END DO   
+                WRITE(numout,*) '         ', jt, '       ', wgtiau(jt)
+             END DO
              WRITE(numout,*) '         ==================================='
              WRITE(numout,*) '         Time-integrated weight = ', ztotwgt
              WRITE(numout,*) '         ==================================='
           ENDIF
-         
+
       ENDIF
 
@@ -338,5 +338,5 @@
          CALL iom_open( c_asminc, inum )
          !
-         CALL iom_get( inum, 'time'       , zdate_inc   ) 
+         CALL iom_get( inum, 'time'       , zdate_inc   )
          CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb )
          CALL iom_get( inum, 'z_inc_datef', z_inc_datef )
@@ -345,5 +345,5 @@
          !
          IF(lwp) THEN
-            WRITE(numout,*) 
+            WRITE(numout,*)
             WRITE(numout,*) 'asm_inc_init : Assimilation increments valid between dates ', z_inc_dateb,' and ', z_inc_datef
             WRITE(numout,*) '~~~~~~~~~~~~'
@@ -359,5 +359,5 @@
             &                ' not agree with Direct Initialization time' )
 
-         IF ( ln_trainc ) THEN   
+         IF ( ln_trainc ) THEN
             CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 )
             CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 )
@@ -371,7 +371,7 @@
          ENDIF
 
-         IF ( ln_dyninc ) THEN   
-            CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 )              
-            CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 )              
+         IF ( ln_dyninc ) THEN
+            CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 )
+            CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 )
             ! Apply the masks
             u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:)
@@ -382,5 +382,5 @@
             WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0
          ENDIF
-        
+
          IF ( ln_sshinc ) THEN
             CALL iom_get( inum, jpdom_autoglo, 'bckineta', ssh_bkginc, 1 )
@@ -408,5 +408,5 @@
       IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN    ! Apply divergence damping filter
          !                                         !--------------------------------------
-         ALLOCATE( zhdiv(jpi,jpj) ) 
+         ALLOCATE( zhdiv(jpi,jpj) )
          !
          DO jt = 1, nn_divdmp
@@ -427,5 +427,5 @@
                      &               + 0.2_wp * ( zhdiv(ji+1,jj) - zhdiv(ji  ,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
                   v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk)                         &
-                     &               + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj  ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 
+                     &               + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj  ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
                END_2D
             END DO
@@ -433,5 +433,5 @@
          END DO
          !
-         DEALLOCATE( zhdiv ) 
+         DEALLOCATE( zhdiv )
          !
       ENDIF
@@ -454,8 +454,8 @@
          CALL iom_open( c_asmdin, inum )
          !
-         CALL iom_get( inum, 'rdastp', zdate_bkg ) 
+         CALL iom_get( inum, 'rdastp', zdate_bkg )
          !
          IF(lwp) THEN
-            WRITE(numout,*) 
+            WRITE(numout,*)
             WRITE(numout,*) '   ==>>>  Assimilation background state valid at : ', zdate_bkg
             WRITE(numout,*)
@@ -466,5 +466,5 @@
             &                ' not agree with Direct Initialization time' )
          !
-         IF ( ln_trainc ) THEN   
+         IF ( ln_trainc ) THEN
             CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg )
             CALL iom_get( inum, jpdom_autoglo, 'sn', s_bkg )
@@ -473,5 +473,5 @@
          ENDIF
          !
-         IF ( ln_dyninc ) THEN   
+         IF ( ln_dyninc ) THEN
             CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg )
             CALL iom_get( inum, jpdom_autoglo, 'vn', v_bkg )
@@ -501,15 +501,15 @@
       !
    END SUBROUTINE asm_inc_init
-   
-   
+
+
    SUBROUTINE tra_asm_inc( kt, Kbb, Kmm, pts, Krhs )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE tra_asm_inc  ***
-      !!          
+      !!
       !! ** Purpose : Apply the tracer (T and S) assimilation increments
       !!
       !! ** Method  : Direct initialization or Incremental Analysis Updating
       !!
-      !! ** Action  : 
+      !! ** Action  :
       !!----------------------------------------------------------------------
       INTEGER                                  , INTENT(in   ) :: kt             ! Current time step
@@ -523,6 +523,6 @@
       !!----------------------------------------------------------------------
       !
-      ! freezing point calculation taken from oc_fz_pt (but calculated for all depths) 
-      ! used to prevent the applied increments taking the temperature below the local freezing point 
+      ! freezing point calculation taken from oc_fz_pt (but calculated for all depths)
+      ! used to prevent the applied increments taking the temperature below the local freezing point
       DO jk = 1, jpkm1
         CALL eos_fzp( pts(:,:,jk,jp_sal,Kmm), fzptnz(:,:,jk), gdept(:,:,jk,Kmm) )
@@ -539,5 +539,5 @@
             !
             IF(lwp) THEN
-               WRITE(numout,*) 
+               WRITE(numout,*)
                WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
                WRITE(numout,*) '~~~~~~~~~~~~'
@@ -549,9 +549,9 @@
                   ! Do not apply negative increments if the temperature will fall below freezing
                   WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. &
-                     &   pts(:,:,jk,jp_tem,Kmm) + pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 
-                     pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt  
+                     &   pts(:,:,jk,jp_tem,Kmm) + pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) )
+                     pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt
                   END WHERE
                ELSE
-                  pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt  
+                  pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt
                ENDIF
                IF (ln_salfix) THEN
@@ -559,5 +559,5 @@
                   ! minimum value salfixmin
                   WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. &
-                     &   pts(:,:,jk,jp_sal,Kmm) + pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 
+                     &   pts(:,:,jk,jp_sal,Kmm) + pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin )
                      pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt
                   END WHERE
@@ -576,5 +576,5 @@
       ELSEIF ( ln_asmdin ) THEN        ! Direct Initialization
          !                             !--------------------------------------
-         !            
+         !
          IF ( kt == nitdin_r ) THEN
             !
@@ -584,18 +584,18 @@
             IF (ln_temnofreeze) THEN
                ! Do not apply negative increments if the temperature will fall below freezing
-               WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_tem,Kmm) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 
-                  pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
+               WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_tem,Kmm) + t_bkginc(:,:,:) > fzptnz(:,:,:) )
+                  pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:)
                END WHERE
             ELSE
-               pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
+               pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:)
             ENDIF
             IF (ln_salfix) THEN
                ! Do not apply negative increments if the salinity will fall below a specified
                ! minimum value salfixmin
-               WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_sal,Kmm) + s_bkginc(:,:,:) > salfixmin ) 
-                  pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
+               WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_sal,Kmm) + s_bkginc(:,:,:) > salfixmin )
+                  pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:)
                END WHERE
             ELSE
-               pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
+               pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:)
             ENDIF
 
@@ -619,5 +619,5 @@
             DEALLOCATE( s_bkg    )
          ENDIF
-         !  
+         !
       ENDIF
       ! Perhaps the following call should be in step
@@ -630,10 +630,10 @@
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE dyn_asm_inc  ***
-      !!          
+      !!
       !! ** Purpose : Apply the dynamics (u and v) assimilation increments.
       !!
       !! ** Method  : Direct initialization or Incremental Analysis Updating.
       !!
-      !! ** Action  : 
+      !! ** Action  :
       !!----------------------------------------------------------------------
       INTEGER                             , INTENT( in )  ::  kt             ! ocean time-step index
@@ -656,5 +656,5 @@
             !
             IF(lwp) THEN
-               WRITE(numout,*) 
+               WRITE(numout,*)
                WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
                WRITE(numout,*) '~~~~~~~~~~~~'
@@ -676,5 +676,5 @@
       ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization
          !                          !-----------------------------------------
-         !         
+         !
          IF ( kt == nitdin_r ) THEN
             !
@@ -683,5 +683,5 @@
             ! Initialize the now fields with the background + increment
             puu(:,:,:,Kmm) = u_bkg(:,:,:) + u_bkginc(:,:,:)
-            pvv(:,:,:,Kmm) = v_bkg(:,:,:) + v_bkginc(:,:,:)  
+            pvv(:,:,:,Kmm) = v_bkg(:,:,:) + v_bkginc(:,:,:)
             !
             puu(:,:,:,Kbb) = puu(:,:,:,Kmm)         ! Update before fields
@@ -702,10 +702,10 @@
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE ssh_asm_inc  ***
-      !!          
+      !!
       !! ** Purpose : Apply the sea surface height assimilation increment.
       !!
       !! ** Method  : Direct initialization or Incremental Analysis Updating.
       !!
-      !! ** Action  : 
+      !! ** Action  :
       !!----------------------------------------------------------------------
       INTEGER, INTENT(IN) :: kt         ! Current time step
@@ -727,5 +727,5 @@
             !
             IF(lwp) THEN
-               WRITE(numout,*) 
+               WRITE(numout,*)
                WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
                   &  kt,' with IAU weight = ', wgtiau(it)
@@ -779,5 +779,5 @@
       !!                  ***  ROUTINE ssh_asm_div  ***
       !!
-      !! ** Purpose :   ssh increment with z* is incorporated via a correction of the local divergence          
+      !! ** Purpose :   ssh increment with z* is incorporated via a correction of the local divergence
       !!                across all the water column
       !!
@@ -795,15 +795,15 @@
       REAL(wp), DIMENSION(:,:)  , POINTER       ::   ztim     ! local array
       !!----------------------------------------------------------------------
-      ! 
+      !
 #if defined key_asminc
       CALL ssh_asm_inc( kt, Kbb, Kmm ) !==   (calculate increments)
       !
-      IF( ln_linssh ) THEN 
+      IF( ln_linssh ) THEN
          phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t(:,:,1,Kmm) * tmask(:,:,1)
-      ELSE 
+      ELSE
          ALLOCATE( ztim(jpi,jpj) )
          ztim(:,:) = ssh_iau(:,:) / ( ht(:,:) + 1.0 - ssmask(:,:) )
-         DO jk = 1, jpkm1                                 
-            phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk) 
+         DO jk = 1, jpkm1
+            phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk)
          END DO
          !
@@ -818,10 +818,10 @@
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE seaice_asm_inc  ***
-      !!          
+      !!
       !! ** Purpose : Apply the sea ice assimilation increment.
       !!
       !! ** Method  : Direct initialization or Incremental Analysis Updating.
       !!
-      !! ** Action  : 
+      !! ** Action  :
       !!
       !!----------------------------------------------------------------------
@@ -844,9 +844,9 @@
             !
             it = kt - nit000 + 1
-            zincwgt = wgtiau(it)      ! IAU weight for the current time step 
+            zincwgt = wgtiau(it)      ! IAU weight for the current time step
             ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments)
             !
             IF(lwp) THEN
-               WRITE(numout,*) 
+               WRITE(numout,*)
                WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
                WRITE(numout,*) '~~~~~~~~~~~~'
@@ -866,6 +866,6 @@
             !
             ! Nudge sea ice depth to bring it up to a required minimum depth
-            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin ) 
-               zhicifinc(:,:) = (zhicifmin - hm_i(:,:)) * zincwgt    
+            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin )
+               zhicifinc(:,:) = (zhicifmin - hm_i(:,:)) * zincwgt
             ELSEWHERE
                zhicifinc(:,:) = 0.0_wp
@@ -907,8 +907,8 @@
             zofrld (:,:) = 1._wp - at_i(:,:)
             zohicif(:,:) = hm_i(:,:)
-            ! 
+            !
             ! Initialize the now fields the background + increment
             at_i(:,:) = 1. - MIN( MAX( 1.-at_i(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp)
-            at_i_b(:,:) = at_i(:,:) 
+            at_i_b(:,:) = at_i(:,:)
             fr_i(:,:) = at_i(:,:)        ! adjust ice fraction
             !
@@ -916,5 +916,5 @@
             !
             ! Nudge sea ice depth to bring it up to a required minimum depth
-            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin ) 
+            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin )
                zhicifinc(:,:) = zhicifmin - hm_i(:,:)
             ELSEWHERE
@@ -946,5 +946,5 @@
 !#if defined defined key_si3 || defined key_cice
 !
-!            IF (ln_seaicebal ) THEN       
+!            IF (ln_seaicebal ) THEN
 !             !! balancing salinity increments
 !             !! simple case from limflx.F90 (doesn't include a mass flux)
@@ -958,5 +958,5 @@
 !
 !             DO jj = 1, jpj
-!               DO ji = 1, jpi 
+!               DO ji = 1, jpi
 !           ! calculate change in ice and snow mass per unit area
 !           ! positive values imply adding salt to the ocean (results from ice formation)
@@ -969,10 +969,10 @@
 !
 !           ! prevent small mld
-!           ! less than 10m can cause salinity instability 
+!           ! less than 10m can cause salinity instability
 !                 IF (mld < 10) mld=10
 !
-!           ! set to bottom of a level 
+!           ! 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 
+!                   IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN
 !                     mld=gdepw(ji,jj,jk+1)
 !                     jkmax=jk
@@ -981,5 +981,5 @@
 !
 !            ! avoid applying salinity balancing in shallow water or on land
-!            ! 
+!            !
 !
 !            ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m)
@@ -992,9 +992,9 @@
 !
 !           ! put increments in for levels in the mixed layer
-!           ! but prevent salinity below a threshold value 
-!
-!                   DO jk = 1, jkmax              
-!
-!                     IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN 
+!           ! but prevent salinity below a threshold value
+!
+!                   DO jk = 1, jkmax
+!
+!                     IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN
 !                           sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn
 !                           sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn
@@ -1007,7 +1007,7 @@
 !      !
 !      !! Adjust fsalt. A +ve fsalt means adding salt to ocean
-!      !!           fsalt(ji,jj) =  fsalt(ji,jj) + zpmess     ! adjust fsalt  
-!      !!               
-!      !!           emps(ji,jj) = emps(ji,jj) + zpmess        ! or adjust emps (see icestp1d) 
+!      !!           fsalt(ji,jj) =  fsalt(ji,jj) + zpmess     ! adjust fsalt
+!      !!
+!      !!           emps(ji,jj) = emps(ji,jj) + zpmess        ! or adjust emps (see icestp1d)
 !      !!                                                     ! E-P (kg m-2 s-2)
 !      !            emp(ji,jj) = emp(ji,jj) + zpmess          ! E-P (kg m-2 s-2)
@@ -1022,5 +1022,5 @@
       !
    END SUBROUTINE seaice_asm_inc
-   
+
    !!======================================================================
 END MODULE asminc
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdydta.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdydta.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdydta.F90	(revision 12680)
@@ -251,6 +251,6 @@
                ij   = idx_bdy(jbdy)%nbj(ib,igrd)
                DO ik = 1, jpkm1
-                  dta_alias%u2d(ib) =   &
-                     & dta_alias%u2d(ib) + e3u(ii,ij,ik,Kmm) * umask(ii,ij,ik) * dta_alias%u3d(ib,ik)
+                  dta_alias%u2d(ib) = dta_alias%u2d(ib)   &
+                     &              + e3u(ii,ij,ik,Kmm) * umask(ii,ij,ik) * dta_alias%u3d(ib,ik)
                END DO
                dta_alias%u2d(ib) =  dta_alias%u2d(ib) * r1_hu(ii,ij,Kmm)
@@ -265,6 +265,6 @@
                ij   = idx_bdy(jbdy)%nbj(ib,igrd)
                DO ik = 1, jpkm1
-                  dta_alias%v2d(ib) =   &
-                     & dta_alias%v2d(ib) + e3v(ii,ij,ik,Kmm) * vmask(ii,ij,ik) * dta_alias%v3d(ib,ik)
+                  dta_alias%v2d(ib) = dta_alias%v2d(ib)   &
+                     &              + e3v(ii,ij,ik,Kmm) * vmask(ii,ij,ik) * dta_alias%v3d(ib,ik)
                END DO
                dta_alias%v2d(ib) =  dta_alias%v2d(ib) * r1_hv(ii,ij,Kmm)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdydyn.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdydyn.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdydyn.F90	(revision 12680)
@@ -37,5 +37,4 @@
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/C1D/step_c1d.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/C1D/step_c1d.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/C1D/step_c1d.F90	(revision 12680)
@@ -83,5 +83,5 @@
       IF(.NOT.ln_linssh )   CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn,      Naa )  ! after vertical scale factors 
 
-      IF(.NOT.ln_linssh )   CALL wzv           ( kstp, Nbb, Nnn, ww,  Naa )  ! now cross-level velocity 
+      IF(.NOT.ln_linssh )   CALL wzv           ( kstp, Nbb, Nnn, Naa, ww )  ! now cross-level velocity 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! diagnostics and outputs       
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/CRS/crsini.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/CRS/crsini.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/CRS/crsini.F90	(revision 12680)
@@ -35,5 +35,4 @@
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
-
 CONTAINS
    
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 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diaar5.F90	(revision 12680)
@@ -78,6 +78,6 @@
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: zarea_ssh , zbotpres       ! 2D workspace 
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: zpe, z2d                   ! 2D workspace 
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: zrhd , zrhop, ztpot, ze3t  ! 3D workspace
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: z2d, zpe                   ! 2D workspace 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: z3d, zrhd , zrhop, ztpot   ! 3D workspace
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ztsn                       ! 4D workspace
 
@@ -104,8 +104,8 @@
          END DO
          DO jk = 1, jpk
-            ze3t(:,:,jk) =  e3t(:,:,jk,Kmm)
+            z3d(:,:,jk) =  rau0 * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
          END DO 
          CALL iom_put( 'volcello'  , zrhd(:,:,:)  )  ! WARNING not consistent with CMIP DR where volcello is at ca. 2000
-         CALL iom_put( 'masscello' , rau0 * ze3t(:,:,:) * tmask(:,:,:) )  ! ocean mass
+         CALL iom_put( 'masscello' , z3d (:,:,:) )   ! ocean mass
       ENDIF 
       !
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diadct.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diadct.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diadct.F90	(revision 12680)
@@ -11,6 +11,6 @@
    !!            3.4  ! 09/2011 (C Bricaud)
    !!----------------------------------------------------------------------
-   !! does not work with agrif
-#if ! defined key_agrif
+#if ! defined key_agrif       
+   !!                        ==>>  CAUTION: does not work with agrif
    !!----------------------------------------------------------------------
    !!   dia_dct      :  Compute the transport through a sec.
@@ -66,6 +66,5 @@
    TYPE SECTION
       CHARACTER(len=60)                            :: name              ! name of the sec
-      LOGICAL                                      :: llstrpond         ! true if you want the computation of salt and
-                                                                       ! heat transports
+      LOGICAL                                      :: llstrpond         ! true if you want the computation of salt and heat transports
       LOGICAL                                      :: ll_ice_section    ! ice surface and ice volume computation
       LOGICAL                                      :: ll_date_line      ! = T if the section crosses the date-line
@@ -74,9 +73,9 @@
       INTEGER, DIMENSION(nb_point_max)             :: direction         ! vector direction of the point in the section
       CHARACTER(len=40),DIMENSION(nb_class_max)    :: classname         ! characteristics of the class
-      REAL(wp), DIMENSION(nb_class_max)            :: zsigi           ,&! in-situ   density classes    (99 if you don't want)
-                                                      zsigp           ,&! potential density classes    (99 if you don't want)
-                                                      zsal            ,&! salinity classes   (99 if you don't want)
-                                                      ztem            ,&! temperature classes(99 if you don't want)
-                                                      zlay              ! level classes      (99 if you don't want)
+      REAL(wp), DIMENSION(nb_class_max)            :: zsigi             ! in-situ   density classes    (99 if you don't want)
+      REAL(wp), DIMENSION(nb_class_max)            :: zsigp             ! potential density classes    (99 if you don't want)
+      REAL(wp), DIMENSION(nb_class_max)            :: zsal              ! salinity classes   (99 if you don't want)
+      REAL(wp), DIMENSION(nb_class_max)            :: ztem              ! temperature classes(99 if you don't want)
+      REAL(wp), DIMENSION(nb_class_max)            :: zlay              ! level classes      (99 if you don't want)
       REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output
       REAL(wp)                                         :: slopeSection  ! slope of the section
@@ -89,4 +88,5 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d 
    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d  
+
 
    !! * Substitutions
@@ -1122,6 +1122,6 @@
   !!    |               |                  |       interpolation between ptab(I,J,K) and ptab(I,J,K+1)
   !!    |               |                  |       zbis = 
-  !!    |               |                  |      [ e3w_n(I+1,J,K)*ptab(I,J,K) + ( e3w_n(I,J,K) - e3w_n(I+1,J,K) ) * ptab(I,J,K-1) ]
-  !!    |               |                  |      /[ e3w_n(I+1,J,K) + e3w_n(I,J,K) - e3w_n(I+1,J,K) ] 
+  !!    |               |                  |      [ e3w_n(I+1,J,K,NOW)*ptab(I,J,K) + ( e3w_n(I,J,K,NOW) - e3w_n(I+1,J,K,NOW) ) * ptab(I,J,K-1) ]
+  !!    |               |                  |     /[ e3w_n(I+1,J,K,NOW)             +   e3w_n(I,J,K,NOW) - e3w_n(I+1,J,K,NOW) ] 
   !!    |               |                  | 
   !!    |               |                  |    2. Horizontal interpolation: compute value at U/V point
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diahsb.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diahsb.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diahsb.F90	(revision 12680)
@@ -57,5 +57,4 @@
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -159,16 +158,16 @@
       !
       DO jk = 1, jpkm1           ! volume variation (calculated with scale factors)
-         zwrk(:,:,jk) =   surf(:,:) * e3t(:,:,jk,Kmm)*tmask(:,:,jk)   &
+         zwrk(:,:,jk) =   surf    (:,:) * e3t    (:,:,jk,Kmm)*tmask    (:,:,jk)   &
             &           - surf_ini(:,:) * e3t_ini(:,:,jk    )*tmask_ini(:,:,jk)
       END DO
       zdiff_v2 = glob_sum_full( 'diahsb', zwrk(:,:,:) )     ! glob_sum_full needed as tmask and tmask_ini could be different
       DO jk = 1, jpkm1           ! heat content variation
-         zwrk(:,:,jk) = ( surf(:,:)*e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_tem,Kmm)   &
-            &           - surf_ini(:,:)*hc_loc_ini(:,:,jk) )
+         zwrk(:,:,jk) = ( surf    (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_tem,Kmm)   &
+            &           - surf_ini(:,:) *         hc_loc_ini(:,:,jk) )
       END DO
       zdiff_hc = glob_sum_full( 'diahsb', zwrk(:,:,:) )
       DO jk = 1, jpkm1           ! salt content variation
-         zwrk(:,:,jk) = ( surf(:,:)*e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_sal,Kmm)   &
-            &           - surf_ini(:,:)*sc_loc_ini(:,:,jk) )
+         zwrk(:,:,jk) = ( surf    (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_sal,Kmm)   &
+            &           - surf_ini(:,:) *         sc_loc_ini(:,:,jk) )
       END DO
       zdiff_sc = glob_sum_full( 'diahsb', zwrk(:,:,:) )
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diamlr.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diamlr.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diamlr.F90	(revision 12680)
@@ -4,22 +4,24 @@
    !! Management of the IOM context for multiple-linear-regression analysis
    !!======================================================================
-   !! History :       !  2019  (S. Mueller)
+   !! History :  4.0  !  2019  (S. Mueller)   Original code
    !!----------------------------------------------------------------------
 
    USE par_oce        , ONLY :   wp, jpi, jpj
    USE phycst         , ONLY :   rpi
+   USE dom_oce        , ONLY :   adatrj
+   USE tide_mod
+   !
    USE in_out_manager , ONLY :   lwp, numout, ln_timing
    USE iom            , ONLY :   iom_put, iom_use, iom_update_file_name
-   USE dom_oce        , ONLY :   adatrj
    USE timing         , ONLY :   timing_start, timing_stop
 #if defined key_iomput
    USE xios
 #endif
-   USE tide_mod
 
    IMPLICIT NONE
    PRIVATE
 
-   LOGICAL, PUBLIC ::   lk_diamlr = .FALSE.
+   LOGICAL, PUBLIC ::   lk_diamlr = .FALSE.   !:         ===>>>   NOT a DOCTOR norm name :  use l_diamlr
+   !                                                              lk_  is used only for logical controlled by a CPP key
 
    PUBLIC ::   dia_mlr_init, dia_mlr_iom_init, dia_mlr
@@ -42,7 +44,7 @@
       !!
       !!----------------------------------------------------------------------
-
+      !
       lk_diamlr = .TRUE.
-
+      !
       IF(lwp) THEN
          WRITE(numout, *)
@@ -50,6 +52,7 @@
          WRITE(numout, *) '~~~~~~~~~~~~   multiple-linear-regression analysis'
       END IF
-
+      !
    END SUBROUTINE dia_mlr_init
+
 
    SUBROUTINE dia_mlr_iom_init
@@ -396,4 +399,5 @@
    END SUBROUTINE dia_mlr_iom_init
 
+
    SUBROUTINE dia_mlr
       !!----------------------------------------------------------------------
@@ -403,6 +407,6 @@
       !!
       !!----------------------------------------------------------------------
-
       REAL(wp), DIMENSION(jpi,jpj) ::   zadatrj2d
+      !!----------------------------------------------------------------------
 
       IF( ln_timing )   CALL timing_start('dia_mlr')
@@ -411,13 +415,13 @@
       ! (value of adatrj converted to time in units of seconds)
       !
-      ! A 2-dimensional field of constant value is sent, and subsequently used
-      ! directly or transformed to a scalar or a constant 3-dimensional field as
-      ! required.
+      ! A 2-dimensional field of constant value is sent, and subsequently used directly 
+      ! or transformed to a scalar or a constant 3-dimensional field as required.
       zadatrj2d(:,:) = adatrj*86400.0_wp
       IF ( iom_use('diamlr_time') ) CALL iom_put('diamlr_time', zadatrj2d)
-
+      !
       IF( ln_timing )   CALL timing_stop('dia_mlr')
-
+      !
    END SUBROUTINE dia_mlr
 
+   !!======================================================================
 END MODULE diamlr
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diaptr.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diaptr.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diaptr.F90	(revision 12680)
@@ -60,4 +60,5 @@
 
    LOGICAL ::   ll_init = .TRUE.        !: tracers  trend flag (set from namelist in trdini)
+   
    !! * Substitutions
 #  include "do_loop_substitute.h90"
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 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diawri.F90	(revision 12680)
@@ -121,5 +121,5 @@
       REAL(wp)::   zztmp2, zztmpy   !   -      -
       REAL(wp), DIMENSION(jpi,jpj)     ::   z2d   ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d, ze3t, ze3u, ze3v, ze3w   ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace
       !!----------------------------------------------------------------------
       ! 
@@ -137,20 +137,32 @@
       CALL iom_put("e3v_0", e3v_0(:,:,:) )
       !
-      DO jk = 1, jpk
-         ze3t(:,:,jk) =  e3t(:,:,jk,Kmm)
-         ze3u(:,:,jk) =  e3u(:,:,jk,Kmm)
-         ze3v(:,:,jk) =  e3v(:,:,jk,Kmm)
-         ze3w(:,:,jk) =  e3w(:,:,jk,Kmm)
-      END DO 
-      !
-      CALL iom_put( "e3t" , ze3t(:,:,:) )
-      CALL iom_put( "e3u" , ze3u(:,:,:) )
-      CALL iom_put( "e3v" , ze3v(:,:,:) )
-      CALL iom_put( "e3w" , ze3w(:,:,:) )
-      IF( iom_use("e3tdef") )   &
-         CALL iom_put( "e3tdef"  , ( ( ze3t(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) !!st r3t
-
-      IF( ll_wd ) THEN
-         CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying)
+      IF ( iom_use("e3t") .OR. iom_use("e3tdef") ) THEN  ! time-varying e3t
+         DO jk = 1, jpk
+            z3d(:,:,jk) =  e3t(:,:,jk,Kmm)
+         END DO
+         CALL iom_put( "e3t"     ,     z3d(:,:,:) )
+         CALL iom_put( "e3tdef"  , ( ( z3d(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) !!st r3t
+      ENDIF 
+      IF ( iom_use("e3u") ) THEN                         ! time-varying e3u
+         DO jk = 1, jpk
+            z3d(:,:,jk) =  e3u(:,:,jk,Kmm)
+         END DO 
+         CALL iom_put( "e3u" , z3d(:,:,:) )
+      ENDIF
+      IF ( iom_use("e3v") ) THEN                         ! time-varying e3v
+         DO jk = 1, jpk
+            z3d(:,:,jk) =  e3v(:,:,jk,Kmm)
+         END DO 
+         CALL iom_put( "e3v" , z3d(:,:,:) )
+      ENDIF
+      IF ( iom_use("e3w") ) THEN                         ! time-varying e3w
+         DO jk = 1, jpk
+            z3d(:,:,jk) =  e3w(:,:,jk,Kmm)
+         END DO 
+         CALL iom_put( "e3w" , z3d(:,:,:) )
+      ENDIF
+
+      IF( ll_wd ) THEN                                   ! sea surface height (brought back to the reference used for wetting and drying)
+         CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) )
       ELSE
          CALL iom_put( "ssh" , ssh(:,:,Kmm) )              ! sea surface height
@@ -216,6 +228,6 @@
 
       IF( ln_zad_Aimp ) ww = ww + wi               ! Recombine explicit and implicit parts of vertical velocity for diagnostic output
-      !
       CALL iom_put( "woce", ww )                   ! vertical velocity
+
       IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
          ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
@@ -777,6 +789,6 @@
          CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T  )   ! heat content
          CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T  )   ! salt content
-         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface heat content
-         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface salinity content
+         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
+         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
       ELSE
          CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T  )   ! temperature
@@ -930,13 +942,15 @@
       !!----------------------------------------------------------------------
       ! 
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
-      IF(lwp) WRITE(numout,*) '                and named :', cdfile_name, '...nc'
-
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
+         WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
+         WRITE(numout,*) '                and named :', cdfile_name, '...nc'
+      ENDIF 
+      !
       DO jk = 1, jpk
          ze3t(:,:,jk) =  e3t(:,:,jk,Kmm)
       END DO
-
+      !
 #if defined key_si3
      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
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 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/dom_oce.F90	(revision 12680)
@@ -2,5 +2,4 @@
    !!======================================================================
    !!                       ***  MODULE dom_oce  ***
-   !!
    !! ** Purpose :   Define in memory all the ocean space domain variables
    !!======================================================================
@@ -75,5 +74,5 @@
    LOGICAL, PUBLIC ::   l_Iperio, l_Jperio   !   should we explicitely take care I/J periodicity
 
-   !                                 !  domain MPP decomposition parameters
+   !                             !: domain MPP decomposition parameters
    INTEGER             , PUBLIC ::   nimpp, njmpp     !: i- & j-indexes for mpp-subdomain left bottom
    INTEGER             , PUBLIC ::   nreci, nrecj     !: overlap region in i and j
@@ -138,29 +137,30 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3vw_0   !: vw-vert. scale factor [m]
    !                                                        !  time-dependent scale factors
+#if ! defined key_qco
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e3t, e3u, e3v, e3w, e3uw, e3vw  !: vert. scale factor [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   e3f                             !: F-point vert. scale factor [m]
+#endif
    !                                                        !  time-dependent ratio ssh / h_0
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   r3t, r3u, r3v                   !: [-]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   r3f                             !: [-]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   r3t_f, r3u_f, r3v_f             !: [-]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   r3t, r3u, r3v                   !: time-dependent    ratio at t-, u- and v-point [-]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   r3f                             !: mid-time-level    ratio at f-point            [-]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   r3t_f, r3u_f, r3v_f             !: now time-filtered ratio at t-, u- and v-point [-]
 
    !                                                        !  reference depths of cells
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_0  !: t- depth              [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdepw_0  !: w- depth              [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gde3w_0  !: w- depth (sum of e3w) [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   gdept_0  !: t- depth              [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   gdepw_0  !: w- depth              [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   gde3w_0  !: w- depth (sum of e3w) [m]
    !                                                        !  time-dependent depths of cells
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  gdept, gdepw
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  gde3w
-
-   !                                                      !  reference heights of water column
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0  !: t-depth              [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0  !: u-depth              [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hv_0  !: v-depth              [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hf_0  !: f-depth              [m]
-   !                                                      !  reciprocal reference heights of water column
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   r1_ht_0, r1_hu_0, r1_hv_0, r1_hf_0   !: t-depth   [1/m]
-                                                          ! time-dependent heights of water column
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ht                     !: height of water column at T-points [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hu, hv, r1_hu, r1_hv   !: height of water column [m] and reciprocal [1/m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   gdept, gdepw
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   gde3w
+
+   !                                                        !  reference heights of ocean water column and its inverse
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   ht_0, r1_ht_0   !: t-depth        [m] and [1/m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   hu_0, r1_hu_0   !: u-depth        [m] and [1/m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   hv_0, r1_hv_0   !: v-depth        [m] and [1/m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   hf_0, r1_hf_0   !: f-depth        [m] and [1/m]
+   !                                                        ! time-dependent heights of ocean water column
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   ht          !: t-points           [m]
+   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]
 
    INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1)
@@ -178,6 +178,6 @@
    !! ---------------------------------------------------------------------
 !!gm Proposition of new name for top/bottom vertical indices
-!   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mtk_t, mtk_u, mtk_v   !: top first wet T-, U-, V-, F-level (ISF)
-!   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbk_t, mbk_u, mbk_v   !: bottom last wet T-, U- and V-level
+!   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mtk_t, mtk_u, mtk_v   !: top    first wet T-, U-, and V-level (ISF)
+!   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbk_t, mbk_u, mbk_v   !: bottom last  wet T-, U-, and V-level
 !!gm
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbkt, mbku, mbkv   !: bottom last wet T-, U- and V-level
@@ -187,7 +187,7 @@
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mikt, miku, mikv, mikf  !: top first wet T-, U-, V-, F-level           (ISF)
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssumask, ssvmask, ssfmask    !: surface mask at T-,U-, V- and F-pts
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask        !: land/ocean mask at WT-, WU- and WV-pts
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)           ::   ssmask, ssumask, ssvmask, ssfmask   !: surface mask at T-,U-, V- and F-pts
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask            !: land/ocean mask at T-, U-, V- and F-pts
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask                 !: land/ocean mask at WT-, WU- and WV-pts
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tpol, fpol          !: north fold mask (jperio= 3 or 4)
@@ -211,10 +211,10 @@
    REAL(wp), PUBLIC ::   adatrj        !: number of elapsed days since the begining of the whole simulation
    !                                   !: (cumulative duration of previous runs that may have used different time-step size)
-   INTEGER , PUBLIC, DIMENSION(  0: 2) ::   nyear_len     !: length in days of the previous/current/next year
-   INTEGER , PUBLIC, DIMENSION(-11:25) ::   nmonth_len    !: length in days of the months of the current year
-   INTEGER , PUBLIC, DIMENSION(-11:25) ::   nmonth_beg    !: second since Jan 1st 0h of the current year and the half of the months
-   INTEGER , PUBLIC                  ::   nsec1jan000     !: second since Jan 1st 0h of nit000 year and Jan 1st 0h the current year
-   INTEGER , PUBLIC                  ::   nsec000_1jan000   !: second since Jan 1st 0h of nit000 year and nit000
-   INTEGER , PUBLIC                  ::   nsecend_1jan000   !: second since Jan 1st 0h of nit000 year and nitend
+   INTEGER , PUBLIC, DIMENSION(  0: 2) ::   nyear_len         !: length in days of the previous/current/next year
+   INTEGER , PUBLIC, DIMENSION(-11:25) ::   nmonth_len        !: length in days of the months of the current year
+   INTEGER , PUBLIC, DIMENSION(-11:25) ::   nmonth_beg        !: second since Jan 1st 0h of the current year and the half of the months
+   INTEGER , PUBLIC                    ::   nsec1jan000       !: second since Jan 1st 0h of nit000 year and Jan 1st 0h the current year
+   INTEGER , PUBLIC                    ::   nsec000_1jan000   !: second since Jan 1st 0h of nit000 year and nit000
+   INTEGER , PUBLIC                    ::   nsecend_1jan000   !: second since Jan 1st 0h of nit000 year and nitend
 
    !!----------------------------------------------------------------------
@@ -249,13 +249,17 @@
    INTEGER FUNCTION dom_oce_alloc()
       !!----------------------------------------------------------------------
-      INTEGER, DIMENSION(12) :: ierr
+      INTEGER                ::   ii
+      INTEGER, DIMENSION(30) :: ierr
       !!----------------------------------------------------------------------
-      ierr(:) = 0
+      ii = 0   ;   ierr(:) = 0
       !
-      ALLOCATE( mig(jpi), mjg(jpj), STAT=ierr(1) )
-         !
-      ALLOCATE( mi0(jpiglo)   , mi1 (jpiglo),  mj0(jpjglo)   , mj1 (jpjglo) ,     &
-         &      tpol(jpiglo)  , fpol(jpiglo)                                , STAT=ierr(2) )
-         !
+      ii = ii+1
+      ALLOCATE( mig(jpi), mjg(jpj), STAT=ierr(ii) )
+         !
+      ii = ii+1
+      ALLOCATE( mi0 (jpiglo) , mi1 (jpiglo),  mj0(jpjglo) , mj1 (jpjglo) ,     &
+         &      tpol(jpiglo) , fpol(jpiglo)                              , STAT=ierr(ii) )
+         !
+      ii = ii+1
       ALLOCATE( glamt(jpi,jpj) ,    glamu(jpi,jpj) ,  glamv(jpi,jpj) ,  glamf(jpi,jpj) ,     &
          &      gphit(jpi,jpj) ,    gphiu(jpi,jpj) ,  gphiv(jpi,jpj) ,  gphif(jpi,jpj) ,     &
@@ -268,35 +272,50 @@
          &      e1e2v(jpi,jpj) , r1_e1e2v(jpi,jpj) , e1_e2v(jpi,jpj)                   ,     &
          &      e1e2f(jpi,jpj) , r1_e1e2f(jpi,jpj)                                     ,     &
-         &      ff_f (jpi,jpj) ,    ff_t (jpi,jpj)                                     , STAT=ierr(3) )
-         !
+         &      ff_f (jpi,jpj) ,    ff_t (jpi,jpj)                                     , STAT=ierr(ii) )
+         !
+      ii = ii+1
       ALLOCATE( gdept_0(jpi,jpj,jpk)     , gdepw_0(jpi,jpj,jpk)     , gde3w_0(jpi,jpj,jpk) ,      &
-         &      gdept  (jpi,jpj,jpk,jpt) , gdepw  (jpi,jpj,jpk,jpt) , gde3w  (jpi,jpj,jpk) , STAT=ierr(4) )
-         !
-      ALLOCATE( e3t_0(jpi,jpj,jpk)     , e3u_0(jpi,jpj,jpk)     , e3v_0(jpi,jpj,jpk)     , e3f_0(jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk)     ,   &
-         &      e3t  (jpi,jpj,jpk,jpt) , e3u  (jpi,jpj,jpk,jpt) , e3v  (jpi,jpj,jpk,jpt) , e3f  (jpi,jpj,jpk) , e3w  (jpi,jpj,jpk,jpt) ,   &
-         &      e3uw_0(jpi,jpj,jpk)     , e3vw_0(jpi,jpj,jpk)     ,         &
-         &      e3uw  (jpi,jpj,jpk,jpt) , e3vw  (jpi,jpj,jpk,jpt) ,         &
-         &      r3t  (jpi,jpj,jpt)     , r3u  (jpi,jpj,jpt)     , r3v  (jpi,jpj,jpt)     , r3f  (jpi,jpj) ,  &
-         &      r3t_f(jpi,jpj)         , r3u_f(jpi,jpj)         , r3v_f(jpi,jpj)                          ,  STAT=ierr(5) )
-         !
-      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(6)  )
-         !
-      ALLOCATE( risfdep(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(7)  ) 
-         !
-      ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(8) )
-         !
-      ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                        &
-         &      ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) ,    &
-         &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) , STAT=ierr(9) )
-         !
-      ALLOCATE( mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj), mikf(jpi,jpj), STAT=ierr(10) )
-         !
+         &      gdept  (jpi,jpj,jpk,jpt) , gdepw  (jpi,jpj,jpk,jpt) , gde3w  (jpi,jpj,jpk) , STAT=ierr(ii) )
+         !
+      ii = ii+1
+      ALLOCATE(  e3t_0(jpi,jpj,jpk) , e3u_0 (jpi,jpj,jpk) , e3v_0 (jpi,jpj,jpk) , e3f_0(jpi,jpj,jpk) ,      &
+         &       e3w_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk)                      ,  STAT=ierr(ii) )
+         !
+#if ! defined key_qco
+      ii = ii+1
+      ALLOCATE( e3t(jpi,jpj,jpk,jpt) , e3u (jpi,jpj,jpk,jpt) , e3v (jpi,jpj,jpk,jpt) , e3f(jpi,jpj,jpk) ,      &
+         &      e3w(jpi,jpj,jpk,jpt) , e3uw(jpi,jpj,jpk,jpt) , e3vw(jpi,jpj,jpk,jpt)                    ,  STAT=ierr(ii) )
+#endif  
+         !
+      ii = ii+1
+      ALLOCATE( r3t  (jpi,jpj,jpt)   , r3u  (jpi,jpj,jpt)    , r3v  (jpi,jpj,jpt)    , r3f  (jpi,jpj) ,  &
+         &      r3t_f(jpi,jpj)       , r3u_f(jpi,jpj)        , r3v_f(jpi,jpj)                         ,  STAT=ierr(ii) )       
+         !
+      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)  )
+         !
+      ii = ii+1
+      ALLOCATE( risfdep(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(ii)  ) 
+         !
+      ii = ii+1
+      ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(ii) )
+         !
+      ii = ii+1
+      ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                                           &
+         &      ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) ,     &
+         &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) ,                    STAT=ierr(ii) )
+         !
+      ii = ii+1
+      ALLOCATE( mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj), mikf(jpi,jpj), STAT=ierr(ii) )
+         !
+      ii = ii+1
       ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) ,     &
-         &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk) , STAT=ierr(11) )
-         !
-      ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(12) )
+         &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk) , STAT=ierr(ii) )
+         !
+      ii = ii+1
+      ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(ii) )
       !
       dom_oce_alloc = MAXVAL(ierr)
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 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domain.F90	(revision 12680)
@@ -35,5 +35,9 @@
    USE dommsk         ! domain: set the mask system
    USE domwri         ! domain: write the meshmask file
+#if ! defined key_qco
    USE domvvl         ! variable volume
+#else
+   USE domqe          ! variable volume
+#endif
    USE c1d            ! 1D configuration
    USE dyncor_c1d     ! 1D configuration: Coriolis term    (cor_c1d routine)
@@ -77,5 +81,5 @@
       CHARACTER (len=*), INTENT(in) :: cdstr                  ! model: NEMO or SAS. Determines core restart variables
       !
-      INTEGER ::   ji, jj, jk, ik   ! dummy loop indices
+      INTEGER ::   ji, jj, jk, jt   ! dummy loop indices
       INTEGER ::   iconf = 0    ! local integers
       CHARACTER (len=64) ::   cform = "(A12, 3(A13, I7))"
@@ -111,5 +115,5 @@
          CASE( 7 )   ;   WRITE(numout,*) '         (i.e. cyclic east-west and north-south)'
          CASE DEFAULT
-            CALL ctl_stop( 'jperio is out of range' )
+            CALL ctl_stop( 'dom_init:   jperio is out of range' )
          END SELECT
          WRITE(numout,*)     '      Ocean model configuration used:'
@@ -141,5 +145,5 @@
       IF( ln_closea ) CALL dom_clo      ! Read in masks to define closed seas and lakes
 
-      CALL dom_zgr( ik_top, ik_bot )    ! Vertical mesh and bathymetry
+      CALL dom_zgr( ik_top, ik_bot )    ! Vertical mesh and bathymetry (return top and bottom ocean t-level indices)
 
       CALL dom_msk( ik_top, ik_bot )    ! Masks
@@ -160,37 +164,46 @@
       r1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp -  ssvmask(:,:) )
       r1_hf_0(:,:) = ssfmask(:,:) / ( hf_0(:,:) + 1._wp -  ssfmask(:,:) )
+
       !
       !           !==  time varying part of coordinate system  ==!
       !
       IF( ln_linssh ) THEN       != Fix in time : set to the reference one for all
-      !
-         !       before        !          now          !       after         !
-            gdept(:,:,:,Kbb) = gdept_0  ;   gdept(:,:,:,Kmm) = gdept_0   ;   gdept(:,:,:,Kaa) = gdept_0   ! depth of grid-points
-            gdepw(:,:,:,Kbb) = gdepw_0  ;   gdepw(:,:,:,Kmm) = gdepw_0   ;   gdepw(:,:,:,Kaa) = gdepw_0   !
-                                            gde3w            = gde3w_0   !        ---          !
-         !
-              e3t(:,:,:,Kbb) =   e3t_0  ;     e3t(:,:,:,Kmm) =   e3t_0   ;   e3t(:,:,:,Kaa) =  e3t_0    ! scale factors
-              e3u(:,:,:,Kbb) =   e3u_0  ;     e3u(:,:,:,Kmm) =   e3u_0   ;   e3u(:,:,:,Kaa) =  e3u_0    !
-              e3v(:,:,:,Kbb) =   e3v_0  ;     e3v(:,:,:,Kmm) =   e3v_0   ;   e3v(:,:,:,Kaa) =  e3v_0    !
-                                     e3f =   e3f_0   !        ---          !
-              e3w(:,:,:,Kbb) =   e3w_0  ;     e3w(:,:,:,Kmm) =   e3w_0   ;    e3w(:,:,:,Kaa) =   e3w_0   !
-             e3uw(:,:,:,Kbb) =  e3uw_0  ;    e3uw(:,:,:,Kmm) =  e3uw_0   ;   e3uw(:,:,:,Kaa) =  e3uw_0   !
-             e3vw(:,:,:,Kbb) =  e3vw_0  ;    e3vw(:,:,:,Kmm) =  e3vw_0   ;   e3vw(:,:,:,Kaa) =  e3vw_0   !
-         !
-! !!st new variable h1_hu_0 h1_hv_0
-!          z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )     ! _i mask due to ISF
-!          z1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) )
-         !
-         !        before       !          now          !       after         !
-                                      ht =    ht_0   !                     ! water column thickness
-               hu(:,:,Kbb) =    hu_0  ;      hu(:,:,Kmm) =    hu_0   ;    hu(:,:,Kaa) =    hu_0   !
-               hv(:,:,Kbb) =    hv_0  ;      hv(:,:,Kmm) =    hv_0   ;    hv(:,:,Kaa) =    hv_0   !
-            r1_hu(:,:,Kbb) = r1_hu_0  ;   r1_hu(:,:,Kmm) = r1_hu_0   ; r1_hu(:,:,Kaa) = r1_hu_0   ! inverse of water column thickness
-            r1_hv(:,:,Kbb) = r1_hv_0  ;   r1_hv(:,:,Kmm) = r1_hv_0   ; r1_hv(:,:,Kaa) = r1_hv_0   !
-         !
+         !
+         DO jt = 1, jpt                         ! depth of t- and w-grid-points
+            gdept(:,:,:,jt) = gdept_0(:,:,:)
+            gdepw(:,:,:,jt) = gdepw_0(:,:,:)
+         END DO
+            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(:,:,:)
+            e3u(:,:,:,jt) =  e3u_0(:,:,:)
+            e3v(:,:,:,jt) =  e3v_0(:,:,:)
+            e3w(:,:,:,jt) =  e3w_0(:,:,:)
+            e3uw(:,:,:,jt) = e3uw_0(:,:,:)
+            e3vw(:,:,:,jt) = e3vw_0(:,:,:)
+         END DO
+            e3f(:,:,:)    =  e3f_0(:,:,:)
+         !
+#endif
+         !
+         DO jt = 1, jpt                         ! water column thickness and its inverse
+            hu(:,:,jt)    =    hu_0(:,:)
+            hv(:,:,jt)    =    hv_0(:,:)
+            r1_hu(:,:,jt) = r1_hu_0(:,:)
+            r1_hv(:,:,jt) = r1_hv_0(:,:)
+         END DO
+            ht(:,:) =    ht_0(:,:)
          !
       ELSE                       != time varying : initialize before/now/after variables
          !
-         IF( .NOT.l_offline )  CALL dom_vvl_init( Kbb, Kmm, Kaa )
+#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
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/dommsk.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/dommsk.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/dommsk.F90	(revision 12680)
@@ -143,5 +143,5 @@
       ! it masks boundaries (bathy=0) where needed depending on the configuration (closed, periodic...)
       CALL lbc_lnk( 'dommsk', tmask  , 'T', 1._wp )      ! Lateral boundary conditions
-
+      
      ! Mask corrections for bdy (read in mppini2)
       READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903)
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 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90	(revision 12680)
@@ -9,31 +9,31 @@
    !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
    !!            4.1  !  2019-08  (A. Coward, D. Storkey) rename dom_vvl_sf_swp -> dom_vvl_sf_update for new timestepping
-   !!            4.x  ! 2020-02  (G. Madec, S. Techene) pure z* (quasi-eulerian) coordinate
+   !!            4.x  !  2020-02  (G. Madec, S. Techene) pure z* (quasi-eulerian) coordinate
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   dom_qe_init     : define initial vertical scale factors, depths and column thickness
-   !!   dom_qe_sf_nxt   : Compute next vertical scale factors
-   !!   dom_qe_sf_update   : Swap vertical scale factors and update the vertical grid
+   !!   dom_qe_init   : define initial vertical scale factors, depths and column thickness
+   !!   dom_qe_sf_nxt : Compute next vertical scale factors
+   !!   dom_qe_sf_update: Swap vertical scale factors and update the vertical grid
    !!   dom_qe_interpol : Interpolate vertical scale factors from one grid point to another
-   !!   dom_qe_r3c      : Compute ssh/h_0 ratioat t-, u-, v-, and optionally f-points
-   !!   dom_qe_rst      : read/write restart file
-   !!   dom_qe_ctl      : Check the vvl options
+   !!   dom_qe_r3c    : Compute ssh/h_0 ratioat t-, u-, v-, and optionally f-points
+   !!       qe_rst_read : read/write restart file
+   !!   dom_qe_ctl    : Check the vvl options
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers
-   USE phycst          ! physical constant
-   USE dom_oce         ! ocean space and time domain
+   USE oce            ! ocean dynamics and tracers
+   USE phycst         ! physical constant
+   USE dom_oce        ! ocean space and time domain
    USE dynadv  , ONLY : ln_dynadv_vec
-   USE isf_oce         ! iceshelf cavities
-   USE sbc_oce         ! ocean surface boundary condition
-   USE wet_dry         ! wetting and drying
-   USE usrdef_istate   ! user defined initial state (wad only)
-   USE restart         ! ocean restart
+   USE isf_oce        ! iceshelf cavities
+   USE sbc_oce        ! ocean surface boundary condition
+   USE wet_dry        ! wetting and drying
+   USE usrdef_istate  ! user defined initial state (wad only)
+   USE restart        ! ocean restart
    !
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O manager library
-   USE lib_mpp         ! distributed memory computing library
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE timing          ! Timing
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
+   USE lib_mpp        ! distributed memory computing library
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE timing         ! Timing
 
    IMPLICIT NONE
@@ -42,7 +42,8 @@
    PUBLIC  dom_qe_init       ! called by domain.F90
    PUBLIC  dom_qe_zgr        ! called by isfcpl.F90
-   PUBLIC  dom_qe_sf_nxt     ! called by steplf.F90
-   PUBLIC  dom_qe_sf_update  ! called by steplf.F90
+!!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
    PUBLIC  dom_qe_r3c        ! called by steplf.F90
 
@@ -102,13 +103,13 @@
       !
       !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
-      CALL dom_qe_rst( nit000, Kbb, Kmm, 'READ' )
-      e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
+      CALL qe_rst_read( nit000, Kbb, Kmm )
+!!st      e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
       !
       CALL dom_qe_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column
       !
-      IF(lwxios) THEN   ! define variables in restart file when writing with XIOS
-         CALL iom_set_rstw_var_active('e3t_b')
-         CALL iom_set_rstw_var_active('e3t_n')
-      ENDIF
+      ! IF(lwxios) THEN   ! define variables in restart file when writing with XIOS
+      !    CALL iom_set_rstw_var_active('e3t_b')
+      !    CALL iom_set_rstw_var_active('e3t_n')
+      ! ENDIF
       !
    END SUBROUTINE dom_qe_init
@@ -147,28 +148,30 @@
       CALL dom_qe_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) )
       !
-      DO jk = 1, jpkm1                    ! Horizontal interpolation of e3t
-         e3t(:,:,jk,Kbb) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) * tmask(:,:,jk) )   ! Kbb time level
-         e3u(:,:,jk,Kbb) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) * umask(:,:,jk) )
-         e3v(:,:,jk,Kbb) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) * vmask(:,:,jk) )
-         e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) * tmask(:,:,jk) )   ! Kmm time level
-         e3u(:,:,jk,Kmm) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) * umask(:,:,jk) )
-         e3v(:,:,jk,Kmm) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) * vmask(:,:,jk) )
-         e3f(:,:,jk)     = e3f_0(:,:,jk) * ( 1._wp + r3f(:,:)     * fmask(:,:,jk) )
-      END DO
-      !
-      DO jk = 1, jpk                      ! Vertical interpolation of e3t,u,v
-         !                                   ! The ratio does not have to be masked at w-level
-         e3w (:,:,jk,Kbb) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )   ! Kbb time level
-         e3uw(:,:,jk,Kbb) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) )
-         e3vw(:,:,jk,Kbb) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) )
-         e3w (:,:,jk,Kmm) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )   ! Kmm time level
-         e3uw(:,:,jk,Kmm) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) )
-         e3vw(:,:,jk,Kmm) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) )
-      END DO
-      !
-      ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
-      e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm)
-      e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm)
-      e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm)
+! !!st
+!       DO jk = 1, jpkm1                    ! Horizontal interpolation of e3t
+!          e3t(:,:,jk,Kbb) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) * tmask(:,:,jk) )   ! Kbb time level
+!          e3u(:,:,jk,Kbb) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) * umask(:,:,jk) )
+!          e3v(:,:,jk,Kbb) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) * vmask(:,:,jk) )
+!          e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) * tmask(:,:,jk) )   ! Kmm time level
+!          e3u(:,:,jk,Kmm) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) * umask(:,:,jk) )
+!          e3v(:,:,jk,Kmm) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) * vmask(:,:,jk) )
+!          e3f(:,:,jk)     = e3f_0(:,:,jk) * ( 1._wp + r3f(:,:)     * fmask(:,:,jk) )
+!       END DO
+!       !
+!       DO jk = 1, jpk                      ! Vertical interpolation of e3t,u,v
+!          !                                   ! The ratio does not have to be masked at w-level
+!          e3w (:,:,jk,Kbb) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )   ! Kbb time level
+!          e3uw(:,:,jk,Kbb) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) )
+!          e3vw(:,:,jk,Kbb) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) )
+!          e3w (:,:,jk,Kmm) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )   ! Kmm time level
+!          e3uw(:,:,jk,Kmm) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) )
+!          e3vw(:,:,jk,Kmm) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) )
+!       END DO
+!       !
+!       ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
+!       e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm)
+!       e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm)
+!       e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm)
+!!st end
       !
       !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
@@ -221,20 +224,86 @@
    END SUBROUTINE dom_qe_zgr
 
-
-   SUBROUTINE dom_qe_sf_nxt( kt, Kbb, Kmm, Kaa, kcall )
+! !!st
+!    SUBROUTINE dom_qe_sf_nxt( kt, Kbb, Kmm, Kaa, kcall )
+!       !!----------------------------------------------------------------------
+!       !!                ***  ROUTINE dom_qe_sf_nxt  ***
+!       !!
+!       !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
+!       !!                 tranxt and dynspg routines
+!       !!
+!       !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
+!       !!
+!       !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
+!       !!               - tilde_e3t_a: after increment of vertical scale factor
+!       !!                              in z_tilde case
+!       !!               - e3(t/u/v)_a
+!       !!
+!       !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
+!       !!----------------------------------------------------------------------
+!       INTEGER, INTENT( in )           ::   kt             ! time step
+!       INTEGER, INTENT( in )           ::   Kbb, Kmm, Kaa  ! time step
+!       INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence
+!       !
+!       INTEGER                ::   ji, jj, jk            ! dummy loop indices
+!       INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
+!       REAL(wp)               ::   z2dt, z_tmin, z_tmax  ! local scalars
+!       LOGICAL                ::   ll_do_bclinic         ! local logical
+!       REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
+!       !!----------------------------------------------------------------------
+!       !
+!       IF( ln_linssh )   RETURN      ! No calculation in linear free surface
+!       !
+!       IF( ln_timing )   CALL timing_start('dom_qe_sf_nxt')
+!       !
+!       IF( kt == nit000 ) THEN
+!          IF(lwp) WRITE(numout,*)
+!          IF(lwp) WRITE(numout,*) 'dom_qe_sf_nxt : compute after scale factors'
+!          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
+!       ENDIF
+!
+!
+!       ! ******************************* !
+!       ! After acale factors at t-points !
+!       ! ******************************* !
+!       !                                                ! --------------------------------------------- !
+!       !                                                ! z_star coordinate and barotropic z-tilde part !
+!       !                                                ! --------------------------------------------- !
+!       !
+!       !
+!       ! *********************************** !
+!       ! After scale factors at u- v- points !
+!       ! *********************************** !
+!       !
+!       DO jk = 1, jpkm1
+!          e3t(:,:,jk,Kaa) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kaa) * tmask(:,:,jk) )
+!          e3u(:,:,jk,Kaa) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kaa) * umask(:,:,jk) )
+!          e3v(:,:,jk,Kaa) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kaa) * vmask(:,:,jk) )
+!       END DO
+!       !
+!       ! *********************************** !
+!       ! 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_qe_sf_nxt')
+!       !
+!    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 scale factors used in tra_zdf, dynnxt,
+      !! ** Purpose :  - compute the after water heigh used in tra_zdf, dynnxt,
       !!                 tranxt and dynspg routines
       !!
-      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
-      !!
-      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
-      !!               - tilde_e3t_a: after increment of vertical scale factor
-      !!                              in z_tilde case
-      !!               - e3(t/u/v)_a
-      !!
-      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
+      !! ** 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
@@ -242,39 +311,15 @@
       INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence
       !
-      INTEGER                ::   ji, jj, jk            ! dummy loop indices
-      INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
-      REAL(wp)               ::   z2dt, z_tmin, z_tmax  ! local scalars
-      LOGICAL                ::   ll_do_bclinic         ! local logical
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
       !!----------------------------------------------------------------------
       !
       IF( ln_linssh )   RETURN      ! No calculation in linear free surface
       !
-      IF( ln_timing )   CALL timing_start('dom_qe_sf_nxt')
+      IF( ln_timing )   CALL timing_start('dom_h_nxt')
       !
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dom_qe_sf_nxt : compute after scale factors'
+         IF(lwp) WRITE(numout,*) 'dom_h_nxt : compute after scale factors'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
       ENDIF
-
-
-      ! ******************************* !
-      ! After acale factors at t-points !
-      ! ******************************* !
-      !                                                ! --------------------------------------------- !
-      !                                                ! z_star coordinate and barotropic z-tilde part !
-      !                                                ! --------------------------------------------- !
-      !
-      !
-      ! *********************************** !
-      ! After scale factors at u- v- points !
-      ! *********************************** !
-      !
-      DO jk = 1, jpkm1
-         e3t(:,:,jk,Kaa) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kaa) * tmask(:,:,jk) )
-         e3u(:,:,jk,Kaa) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kaa) * umask(:,:,jk) )
-         e3v(:,:,jk,Kaa) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kaa) * vmask(:,:,jk) )
-      END DO
       !
       ! *********************************** !
@@ -287,52 +332,118 @@
       r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
       !
-      IF( ln_timing )   CALL timing_stop('dom_qe_sf_nxt')
-      !
-   END SUBROUTINE dom_qe_sf_nxt
-
-
-   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
 
-
-   SUBROUTINE dom_qe_sf_update( kt, Kbb, Kmm, Kaa )
+! !!st
+!    SUBROUTINE dom_qe_sf_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:
+!       !!                    e3(u/v)_b
+!       !!                    e3w(:,:,:,Kmm)
+!       !!                    e3(u/v)w_b
+!       !!                    e3(u/v)w_n
+!       !!                    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
+!       ! --------------------------------------
+!       ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt
+!       ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
+!
+!
+!       ! Scale factor computation
+!       DO jk = 1, jpk             ! Horizontal interpolation
+!          e3f(:,:,jk)      =  e3f_0(:,:,jk) * ( 1._wp + r3f(:,:) * fmask(:,:,jk) )   ! Kmm time level
+!          !                       ! Vertical interpolation
+! 	      !                                   ! The ratio does not have to be masked at w-level
+!          e3w (:,:,jk,Kmm) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )   ! Kmm time level
+!          e3uw(:,:,jk,Kmm) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) )
+!          e3vw(:,:,jk,Kmm) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) )
+!          e3w (:,:,jk,Kbb) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )   ! Kbb time level
+!          e3uw(:,:,jk,Kbb) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) )
+!          e3vw(:,:,jk,Kbb) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) )
+!       END DO
+!
+!
+!       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( lrst_oce  )   CALL dom_qe_rst( kt, Kbb, Kmm, 'WRITE' )
+!       !
+!       IF( ln_timing )   CALL timing_stop('dom_qe_sf_update')
+!       !
+!    END SUBROUTINE dom_qe_sf_update
+!!st end
+
+   SUBROUTINE dom_h_update( kt, Kbb, Kmm, Kaa )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE dom_qe_sf_update  ***
@@ -346,8 +457,4 @@
       !!
       !! ** Action  :  - Recompute:
-      !!                    e3(u/v)_b
-      !!                    e3w(:,:,:,Kmm)
-      !!                    e3(u/v)w_b
-      !!                    e3(u/v)w_n
       !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w
       !!                    h(u/v) and h(u/v)r
@@ -377,21 +484,5 @@
       ! Horizontal scale factor interpolations
       ! --------------------------------------
-      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt
       ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
-
-
-      ! Scale factor computation
-      DO jk = 1, jpk             ! Horizontal interpolation
-         e3f(:,:,jk)      =  e3f_0(:,:,jk) * ( 1._wp + r3f(:,:) * fmask(:,:,jk) )   ! Kmm time level
-         !                       ! Vertical interpolation
-	      !                                   ! The ratio does not have to be masked at w-level
-         e3w (:,:,jk,Kmm) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )   ! Kmm time level
-         e3uw(:,:,jk,Kmm) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) )
-         e3vw(:,:,jk,Kmm) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) )
-         e3w (:,:,jk,Kbb) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )   ! Kbb time level
-         e3uw(:,:,jk,Kbb) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) )
-         e3vw(:,:,jk,Kbb) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) )
-      END DO
-
 
       IF( ln_isf ) THEN          !** IceShelF cavities
@@ -399,5 +490,5 @@
       !                             ! 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 !
+   !!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
@@ -417,5 +508,5 @@
       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 !
+   !!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) )
@@ -435,9 +526,7 @@
       ! write restart file
       ! ==================
-      IF( lrst_oce  )   CALL dom_qe_rst( kt, Kbb, Kmm, 'WRITE' )
-      !
       IF( ln_timing )   CALL timing_stop('dom_qe_sf_update')
       !
-   END SUBROUTINE dom_qe_sf_update
+   END SUBROUTINE dom_h_update
 
 
@@ -507,19 +596,16 @@
 
 
-   SUBROUTINE dom_qe_rst( kt, Kbb, Kmm, cdrw )
+   SUBROUTINE qe_rst_read( kt, Kbb, Kmm )
       !!---------------------------------------------------------------------
-      !!                   ***  ROUTINE dom_qe_rst  ***
-      !!
-      !! ** Purpose :   Read or write VVL file in restart file
+      !!                   ***  ROUTINE qe_rst_read  ***
+      !!
+      !! ** Purpose :   Read ssh in restart file
       !!
       !! ** Method  :   use of IOM library
-      !!                if the restart does not contain vertical scale factors,
-      !!                they are set to the _0 values
-      !!                if the restart does not contain vertical scale factors increments (z_tilde),
-      !!                they are set to 0.
+      !!                if the restart does not contain ssh,
+      !!                it is set to the _0 values.
       !!----------------------------------------------------------------------
       INTEGER         , INTENT(in) ::   kt        ! ocean time-step
       INTEGER         , INTENT(in) ::   Kbb, Kmm  ! ocean time level indices
-      CHARACTER(len=*), INTENT(in) ::   cdrw      ! "READ"/"WRITE" flag
       !
       INTEGER ::   ji, jj, jk
@@ -527,12 +613,9 @@
       !!----------------------------------------------------------------------
       !
-      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
-         !                                   ! ===============
          IF( ln_rstart ) THEN                   !* Read the restart file
             CALL rst_read_open                  !  open the restart file if necessary
-            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    )
             !
-            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
-            id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
+            id1 = iom_varid( numror, 'sshb', ldstop = .FALSE. )
+            id2 = iom_varid( numror, 'sshn', ldstop = .FALSE. )
             !
             !                             ! --------- !
@@ -541,39 +624,34 @@
             !
             IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
+               CALL iom_get( numror, jpdom_autoglo, 'sshb'   , ssh(:,:,Kbb), ldxios = lrxios    )
+               CALL iom_get( numror, jpdom_autoglo, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    )
                ! needed to restart if land processor not computed
-               IF(lwp) write(numout,*) 'dom_qe_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files'
-               WHERE ( tmask(:,:,:) == 0.0_wp )
-                  e3t(:,:,:,Kmm) = e3t_0(:,:,:)
-                  e3t(:,:,:,Kbb) = e3t_0(:,:,:)
+               IF(lwp) write(numout,*) 'qe_rst_read : ssh(:,:,Kbb) and ssh(:,:,Kmm) found in restart files'
+               WHERE ( ssmask(:,:) == 0.0_wp )   !!gm/st ==> sm should not be necessary on ssh when it was required on e3
+                  ssh(:,:,Kmm) = 0._wp
+                  ssh(:,:,Kbb) = 0._wp
                END WHERE
                IF( neuler == 0 ) THEN
-                  e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
+                  ssh(:,:,Kbb) = ssh(:,:,Kmm)
                ENDIF
             ELSE IF( id1 > 0 ) THEN
-               IF(lwp) write(numout,*) 'dom_qe_rst WARNING : e3t(:,:,:,Kmm) not found in restart files'
-               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
+               IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kmm) not found in restart files'
+               IF(lwp) write(numout,*) 'sshn set equal to sshb.'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
-               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
+               CALL iom_get( numror, jpdom_autoglo, 'sshb', ssh(:,:,Kbb), ldxios = lrxios )
+               ssh(:,:,Kmm) = ssh(:,:,Kbb)
                neuler = 0
             ELSE IF( id2 > 0 ) THEN
-               IF(lwp) write(numout,*) 'dom_qe_rst WARNING : e3t(:,:,:,Kbb) not found in restart files'
-               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
+               IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kbb) not found in restart files'
+               IF(lwp) write(numout,*) 'sshb set equal to sshn.'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
-               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
+               CALL iom_get( numror, jpdom_autoglo, 'sshn', ssh(:,:,Kmm), ldxios = lrxios )
+               ssh(:,:,Kbb) = ssh(:,:,Kmm)
                neuler = 0
             ELSE
-               IF(lwp) write(numout,*) 'dom_qe_rst WARNING : e3t(:,:,:,Kmm) not found in restart file'
-               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
+               IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kmm) not found in restart file'
+               IF(lwp) write(numout,*) 'ssh_b and ssh_n set to zero'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
-               DO jk = 1, jpk
-                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
-                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
-                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
-               END DO
-               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
+               ssh(:,:,:) = 0._wp
                neuler = 0
             ENDIF
@@ -583,16 +661,14 @@
             IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
                !
-               IF( cn_cfg == 'wad' ) THEN
-                  ! Wetting and drying test case
+               IF( cn_cfg == 'wad' ) THEN            ! Wetting and drying test case
                   CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
-                  ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
-                  ssh (:,:,Kmm)     = ssh(:,:,Kbb)
-                  uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
-                  vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
-               ELSE
-                  ! if not test case
+                  ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
+                  ssh(:,:    ,Kmm) = ssh(:,:    ,Kbb)
+                  uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb)
+                  vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb)
+               ELSE                                  ! if not test case
                   ssh(:,:,Kmm) = -ssh_ref
                   ssh(:,:,Kbb) = -ssh_ref
-
+                  !
                   DO_2D_11_11
                      IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
@@ -601,16 +677,10 @@
                      ENDIF
                   END_2D
-               ENDIF !If test case else
-
-               ! Adjust vertical metrics for all wad
-               DO jk = 1, jpk
-                  e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) * tmask(:,:,jk) )
-               END DO
-               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
+               ENDIF
 
                DO ji = 1, jpi
                   DO jj = 1, jpj
                      IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
-                       CALL ctl_stop( 'dom_qe_rst: ht_0 must be positive at potentially wet points' )
+                       CALL ctl_stop( 'qe_rst_read: ht_0 must be positive at potentially wet points' )
                      ENDIF
                   END DO
@@ -623,31 +693,11 @@
 !               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
 !               !
-!               DO jk=1,jpk
-!                  e3t(:,:,jk,Kbb) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) &
-!                     &            / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk)
-!               END DO
-!               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
-                ssh(:,:,Kmm)=0._wp
-                e3t(:,:,:,Kmm)=e3t_0(:,:,:)
-                e3t(:,:,:,Kbb)=e3t_0(:,:,:)
+                ssh(:,:,:) = 0._wp
                !
             ENDIF           ! end of ll_wd edits
             !
          ENDIF
-         !
-      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
-         !                                   ! ===================
-         IF(lwp) WRITE(numout,*) '---- dom_qe_rst ----'
-         IF( lwxios ) CALL iom_swap(      cwxios_context          )
-         !                                           ! --------- !
-         !                                           ! all cases !
-         !                                           ! --------- !
-         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios )
-         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios )
-         !
-         IF( lwxios ) CALL iom_swap(      cxios_context          )
-      ENDIF
-      !
-   END SUBROUTINE dom_qe_rst
+      !
+   END SUBROUTINE qe_rst_read
 
 
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domvvl.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domvvl.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domvvl.F90	(revision 12680)
@@ -12,12 +12,4 @@
    !!----------------------------------------------------------------------
 
-   !!----------------------------------------------------------------------
-   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
-   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors
-   !!   dom_vvl_sf_update   : Swap vertical scale factors and update the vertical grid
-   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
-   !!   dom_vvl_rst      : read/write restart file
-   !!   dom_vvl_ctl      : Check the vvl options
-   !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
    USE phycst          ! physical constant
@@ -36,11 +28,5 @@
    IMPLICIT NONE
    PRIVATE
-
-   PUBLIC  dom_vvl_init       ! called by domain.F90
-   PUBLIC  dom_vvl_zgr        ! called by isfcpl.F90
-   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90
-   PUBLIC  dom_vvl_sf_update  ! called by step.F90
-   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90
-
+   
    !                                                      !!* Namelist nam_vvl
    LOGICAL , PUBLIC :: ln_vvl_zstar           = .FALSE.    ! zstar  vertical coordinate
@@ -64,4 +50,28 @@
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! retoring period for low freq. divergence
 
+#if defined key_qco
+   !!----------------------------------------------------------------------
+   !!   'key_qco'      EMPTY MODULE      Quasi-Eulerian vertical coordonate
+   !!----------------------------------------------------------------------
+#else
+   !!----------------------------------------------------------------------
+   !!   Default key      Old management of time varying vertical coordinate
+   !!----------------------------------------------------------------------
+   
+   !!----------------------------------------------------------------------
+   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
+   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors
+   !!   dom_vvl_sf_update   : Swap vertical scale factors and update the vertical grid
+   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
+   !!   dom_vvl_rst      : read/write restart file
+   !!   dom_vvl_ctl      : Check the vvl options
+   !!----------------------------------------------------------------------
+
+   PUBLIC  dom_vvl_init       ! called by domain.F90
+   PUBLIC  dom_vvl_zgr        ! called by isfcpl.F90
+   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90
+   PUBLIC  dom_vvl_sf_update  ! called by step.F90
+   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90
+   
    !! * Substitutions
 #  include "do_loop_substitute.h90"
@@ -136,5 +146,6 @@
       !
    END SUBROUTINE dom_vvl_init
-   !
+
+
    SUBROUTINE dom_vvl_zgr(Kbb, Kmm, Kaa)
       !!----------------------------------------------------------------------
@@ -1039,4 +1050,6 @@
    END SUBROUTINE dom_vvl_ctl
 
+#endif
+
    !!======================================================================
 END MODULE domvvl
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 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domzgr_substitute.h90	(revision 12680)
@@ -12,11 +12,11 @@
 !!----------------------------------------------------------------------
 #if defined key_qco
-#   define  e3t(i,j,k,t)   (e3t_0(i,j,k)*(1.+r3t(i,j,t)*tmask(i,j,k)))
-#   define  e3u(i,j,k,t)   (e3u_0(i,j,k)*(1.+r3u(i,j,t)*umask(i,j,k)))
-#   define  e3v(i,j,k,t)   (e3v_0(i,j,k)*(1.+r3v(i,j,t)*vmask(i,j,k)))
-#   define  e3f(i,j,k)     (e3f_0(i,j,k)*(1.+r3f(i,j)*fmask(i,j,k)))
-#   define  e3w(i,j,k,t)   (e3w_0(i,j,k)*(1.+r3t(i,j,t)))
-#   define  e3uw(i,j,k,t)  (e3uw_0(i,j,k)*(1.+r3u(i,j,t)))
-#   define  e3vw(i,j,k,t)  (e3vw_0(i,j,k)*(1.+r3v(i,j,t)))
+#   define  e3t(i,j,k,t)   (e3t_0(i,j,k)*(1._wp+r3t(i,j,t)*tmask(i,j,k)))
+#   define  e3u(i,j,k,t)   (e3u_0(i,j,k)*(1._wp+r3u(i,j,t)*umask(i,j,k)))
+#   define  e3v(i,j,k,t)   (e3v_0(i,j,k)*(1._wp+r3v(i,j,t)*vmask(i,j,k)))
+#   define  e3f(i,j,k)     (e3f_0(i,j,k)*(1._wp+r3f(i,j)*fmask(i,j,k)))
+#   define  e3w(i,j,k,t)   (e3w_0(i,j,k)*(1._wp+r3t(i,j,t)))
+#   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)))
 #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 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/istate.F90	(revision 12680)
@@ -43,4 +43,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
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 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatf.F90	(revision 12680)
@@ -58,4 +58,21 @@
 
    PUBLIC    dyn_atf   ! routine called by step.F90
+
+#if defined key_qco
+   !!----------------------------------------------------------------------
+   !!   'key_qco'      EMPTY ROUTINE     Quasi-Eulerian vertical coordonate
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dyn_atf ( kt, Kbb, Kmm, Kaa, puu, pvv, pe3t, pe3u, pe3v )
+      INTEGER                             , INTENT(in   ) :: kt               ! ocean time-step index
+      INTEGER                             , INTENT(in   ) :: Kbb, Kmm, Kaa    ! before and after time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv         ! velocities to be time filtered
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: pe3t, pe3u, pe3v ! scale factors to be time filtered
+
+      WRITE(*,*) 'dyn_atf: You should not have seen this print! error?', kt
+   END SUBROUTINE dyn_atf
+
+#else
 
    !! * Substitutions
@@ -314,4 +331,6 @@
    END SUBROUTINE dyn_atf
 
+#endif
+
    !!=========================================================================
 END MODULE dynatf
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/sshwzv.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/sshwzv.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/sshwzv.F90	(revision 12680)
@@ -114,5 +114,6 @@
       !
 #if defined key_agrif
-      Kbb_a = Kbb; Kmm_a = Kmm; Krhs_a = Kaa; CALL agrif_ssh( kt )
+      Kbb_a = Kbb   ;   Kmm_a = Kmm   ;   Krhs_a = Kaa
+      CALL agrif_ssh( kt )
 #endif
       !
@@ -134,5 +135,5 @@
 
 
-   SUBROUTINE wzv( kt, Kbb, Kmm, pww, Kaa )
+   SUBROUTINE wzv( kt, Kbb, Kmm, Kaa, pww )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE wzv  ***
@@ -151,5 +152,5 @@
       INTEGER                         , INTENT(in)    ::   kt             ! time step
       INTEGER                         , INTENT(in)    ::   Kbb, Kmm, Kaa  ! time level indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pww            ! now vertical velocity
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pww            ! vertical velocity at Kmm
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -165,13 +166,13 @@
          IF(lwp) WRITE(numout,*) '~~~~~ '
          !
-         pww(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
-      ENDIF
-      !                                           !------------------------------!
-      !                                           !     Now Vertical Velocity    !
-      !                                           !------------------------------!
-      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
+         pww(:,:,jpk) = 0._wp           ! bottom boundary condition: w=0 (set once for all)
+      ENDIF
+      !
+      z1_2dt = 1. / ( 2. * rdt )        ! set time step size (Euler/Leapfrog)
       IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
       !
-      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
+      !                                               !===============================!
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      !==  z_tilde and layer cases  ==!
+         !                                            !===============================!
          ALLOCATE( zhdiv(jpi,jpj,jpk) )
          !
@@ -188,14 +189,21 @@
          DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
             ! computation of w
-            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)    &
-               &                         + zhdiv(:,:,jk)                       &
-               &                         + z1_2dt * ( e3t(:,:,jk,Kaa)          &
-               &                                    - e3t(:,:,jk,Kbb) )  ) * tmask(:,:,jk)
+            pww(:,:,jk) = pww(:,:,jk+1) - (   e3t(:,:,jk,Kmm) * hdiv(:,:,jk)   &
+               &                            +                  zhdiv(:,:,jk)   &
+               &                            + z1_2dt * ( e3t(:,:,jk,Kaa)       &
+               &                                       - e3t(:,:,jk,Kbb) )   ) * tmask(:,:,jk)
          END DO
          !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
          DEALLOCATE( zhdiv )
-      ELSE   ! z_star and linear free surface cases
-         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
-            ! computation of w
+         !                                            !=================================!
+      ELSEIF( ln_linssh )   THEN                      !==  linear free surface cases  ==!
+         !                                            !=================================!
+         DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence
+            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)  ) * tmask(:,:,jk)
+         END DO
+         !                                            !==========================================!
+      ELSE                                            !==  Quasi-Eulerian vertical coordinate  ==!   ('key_qco')
+         !                                            !==========================================!
+         DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence
             pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)    &
                &                         + z1_2dt * ( e3t(:,:,jk,Kaa)          &
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/IOM/restart.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/IOM/restart.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/IOM/restart.F90	(revision 12680)
@@ -285,15 +285,8 @@
       !
       IF( neuler == 0 ) THEN                                  ! Euler restart (neuler=0)
-         ts   (:,:,:,:,Kbb) = ts   (:,:,:,:,Kmm)              ! all before fields set to now values
-         uu   (:,:,:  ,Kbb) = uu   (:,:,:  ,Kmm)
-         vv   (:,:,:  ,Kbb) = vv   (:,:,:  ,Kmm)
-         ssh  (:,:    ,Kbb) = ssh  (:,:    ,Kmm)
-         !
-         IF( .NOT.ln_linssh ) THEN
-            DO jk = 1, jpk
-               e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm)
-            END DO
-         ENDIF
-         !
+         ts (:,:,:,:,Kbb) = ts (:,:,:,:,Kmm)                  ! all before fields set to now values
+         uu (:,:,:  ,Kbb) = uu (:,:,:  ,Kmm)
+         vv (:,:,:  ,Kbb) = vv (:,:,:  ,Kmm)
+         ssh(:,:    ,Kbb) = ssh(:,:    ,Kmm)
       ENDIF
       !
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 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfcpl.F90	(revision 12680)
@@ -226,7 +226,5 @@
 #if ! defined key_qco
       DO jk = 1, jpk
-         e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
-             &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
-             &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
+         e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + (ht_0(:,:) + ssh(:,:,Kmm)) * r1_ht_0(:,:) )
       END DO
       e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcice_cice.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcice_cice.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcice_cice.F90	(revision 12680)
@@ -239,7 +239,7 @@
                !
                DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
-                  e3t(:,:,jk,Kmm) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kmm)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
-                  e3t(:,:,jk,Kbb) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kbb)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
-               ENDDO
+                  e3t(:,:,jk,Kmm) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kmm)*r1_ht_0(:,:)*tmask(:,:,jk) )
+                  e3t(:,:,jk,Kbb) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kbb)*r1_ht_0(:,:)*tmask(:,:,jk) )
+               END DO
                e3t(:,:,:,Krhs) = e3t(:,:,:,Kbb)
                ! Reconstruction of all vertical scale factors at now and before time-steps
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatf.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatf.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatf.F90	(revision 12680)
@@ -132,10 +132,6 @@
          ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from pts(Kmm) terms
          DO jk = 1, jpkm1
-            ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa)   &
-               &                                   / e3t(:,:,jk,Kmm)   &
-               &            - pts(:,:,jk,jp_tem,Kmm)) * zfact
-            ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa)   &
-               &                                   / e3t(:,:,jk,Kmm)   &
-               &            - pts(:,:,jk,jp_sal,Kmm)) * zfact
+            ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - pts(:,:,jk,jp_tem,Kmm)) * zfact
+            ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - pts(:,:,jk,jp_sal,Kmm)) * zfact
          END DO
          CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_tot, ztrdt )
@@ -336,6 +332,5 @@
                   IF ( jk == misfkb_cav(ji,jj) ) THEN
                      ztc_f  = ztc_f  - zfact1 * ( risf_cav_tsc(ji,jj,jn) - risf_cav_tsc_b(ji,jj,jn) )  &
-                            &                 * e3t(ji,jj,jk,Kmm) / rhisf_tbl_cav(ji,jj)   &
-                            &                 * rfrac_tbl_cav(ji,jj)
+                            &                 * e3t(ji,jj,jk,Kmm) / rhisf_tbl_cav(ji,jj) * rfrac_tbl_cav(ji,jj)
                   END IF
                END IF
@@ -351,6 +346,5 @@
                   IF ( jk == misfkb_par(ji,jj) ) THEN
                      ztc_f  = ztc_f  - zfact1 * ( risf_par_tsc(ji,jj,jn) - risf_par_tsc_b(ji,jj,jn) )  &
-                            &                 * e3t(ji,jj,jk,Kmm) / rhisf_tbl_par(ji,jj)   &
-                            &                 * rfrac_tbl_par(ji,jj)
+                            &                 * e3t(ji,jj,jk,Kmm) / rhisf_tbl_par(ji,jj) * rfrac_tbl_par(ji,jj)
                   END IF
                END IF
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/step.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/step.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/step.F90	(revision 12680)
@@ -33,5 +33,9 @@
    !!            4.1  !  2019-08  (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme
    !!----------------------------------------------------------------------
-
+#if defined key_qco
+   !!----------------------------------------------------------------------
+   !!   'key_qco'      EMPTY MODULE      Quasi-Eulerian vertical coordonate
+   !!----------------------------------------------------------------------
+#else
    !!----------------------------------------------------------------------
    !!   stp             : OPA system time-stepping
@@ -171,5 +175,5 @@
                             CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor)
       IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn,      Naa )    ! after vertical scale factors 
-                            CALL wzv           ( kstp, Nbb, Nnn, ww,  Naa )    ! now cross-level velocity 
+                            CALL wzv           ( kstp, Nbb, Nnn, Naa, ww  )    ! now 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
@@ -200,5 +204,5 @@
                             CALL dyn_zdf    ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion
       IF( ln_dynspg_ts ) THEN                                                       ! vertical scale factors and vertical velocity need to be updated
-                            CALL wzv        ( kstp, Nbb, Nnn, ww, Naa )             ! now cross-level velocity 
+                            CALL wzv        ( kstp, Nbb, Nnn, Naa, ww )             ! now cross-level velocity 
          IF( ln_zad_Aimp )  CALL wAimp      ( kstp,      Nnn )                      ! Adaptive-implicit vertical advection partitioning
       ENDIF
@@ -339,4 +343,5 @@
    END SUBROUTINE stp
    !
+#endif
    !!======================================================================
 END MODULE step
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 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/stepLF.F90	(revision 12680)
@@ -55,5 +55,4 @@
    !!----------------------------------------------------------------------
    INTEGER, PUBLIC :: Nbb, Nnn, Naa, Nrhs          !! used by nemo_init
-
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
@@ -185,9 +184,9 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
-                            CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor)
+                            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
+      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, ww,  Naa )    ! now cross-level velocity
+                            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
@@ -213,12 +212,12 @@
                                                       ! With split-explicit free surface, since now transports have been updated and ssh(:,:,Nrhs) as well
       IF( ln_dynspg_ts ) THEN                         ! vertical scale factors and vertical velocity need to be updated
-                            CALL div_hor       ( kstp, Nbb, Nnn )                ! Horizontal divergence  (2nd call in time-split case)
-         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)
+                            CALL div_hor    ( kstp, Nbb, Nnn )                ! Horizontal divergence  (2nd call in time-split case)
+         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)
       ENDIF
                             CALL dyn_zdf    ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion
       IF( ln_dynspg_ts ) THEN                                                       ! vertical scale factors and vertical velocity need to be updated
-                            CALL wzv        ( kstp, Nbb, Nnn, ww, Naa )             ! now cross-level velocity
+                            CALL wzv        ( kstp, Nbb, Nnn, Naa, ww )             ! Nnn cross-level velocity
          IF( ln_zad_Aimp )  CALL wAimp      ( kstp,      Nnn )                      ! Adaptive-implicit vertical advection partitioning
       ENDIF
@@ -294,11 +293,10 @@
 !!    place.
 !!
-!!jc2: dynnxt must be the latest call. e3t(:,:,:,Nbb) are indeed updated in that routine
-                         CALL zdyn_ts       ( Nnn, Naa, e3u, e3v, uu, vv )                   ! barotrope ajustment
+                         CALL zdyn_ts       ( Nnn, Naa, uu, vv )                   ! barotrope ajustment
                          CALL finalize_sbc  ( kstp, Nbb, Naa, uu, vv, ts )                   ! boundary condifions
                          CALL ssh_atf       ( kstp, Nbb, Nnn, Naa, ssh )                     ! time filtering of "now" sea surface height
                          CALL dom_qe_r3c    ( ssh(:,:,Nnn), r3t_f, r3u_f, r3v_f )            ! "now" ssh/h_0 ratio from filtrered ssh
                          CALL tra_atf_qco   ( kstp, Nbb, Nnn, Naa, ts )                      ! time filtering of "now" tracer arrays
-                         CALL dyn_atf_qco   ( kstp, Nbb, Nnn, Naa, uu, vv )  ! time filtering of "now" velocities and scale factors
+                         CALL dyn_atf_qco   ( kstp, Nbb, Nnn, Naa, uu, vv  )  ! time filtering of "now" velocities and scale factors
                          r3t(:,:,Nnn) = r3t_f(:,:)
                          r3u(:,:,Nnn) = r3u_f(:,:)
@@ -312,5 +310,6 @@
       Naa = Nrhs
       !
-      IF(.NOT.ln_linssh) CALL dom_qe_sf_update( kstp, Nbb, Nnn, Naa )  ! recompute vertical scale factors
+      !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
       !
       IF( ln_diahsb  )   CALL dia_hsb       ( kstp, Nbb, Nnn )  ! - ML - global conservation diagnostics
@@ -367,5 +366,5 @@
 
 
-   SUBROUTINE zdyn_ts (Kmm, Kaa, pe3u, pe3v, puu, pvv)
+   SUBROUTINE zdyn_ts (Kmm, Kaa, puu, pvv)
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdyn_ts  ***
@@ -380,5 +379,4 @@
       INTEGER                             , INTENT(in   ) :: Kmm, Kaa    ! before and after time level indices
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv         ! velocities
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in   ) :: pe3u, pe3v   ! scale factors
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
@@ -391,9 +389,9 @@
          ! Ensure below that barotropic velocities match time splitting estimate
          ! Compute actual transport and replace it with ts estimate at "after" time step
-         zue(:,:) = pe3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
-         zve(:,:) = pe3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
+         zue(:,:) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
+         zve(:,:) = e3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
          DO jk = 2, jpkm1
-            zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
-            zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
+            zue(:,:) = zue(:,:) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
+            zve(:,:) = zve(:,:) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
          END DO
          DO jk = 1, jpkm1
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/C14/trcsms_c14.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/C14/trcsms_c14.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/C14/trcsms_c14.F90	(revision 12680)
@@ -28,4 +28,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/CFC/trcsms_cfc.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/CFC/trcsms_cfc.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/CFC/trcsms_cfc.F90	(revision 12680)
@@ -49,4 +49,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P2Z/p2zbio.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P2Z/p2zbio.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P2Z/p2zbio.F90	(revision 12680)
@@ -58,4 +58,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P2Z/p2zexp.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P2Z/p2zexp.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P2Z/p2zexp.F90	(revision 12680)
@@ -39,4 +39,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P2Z/p2zopt.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P2Z/p2zopt.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P2Z/p2zopt.F90	(revision 12680)
@@ -40,4 +40,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P2Z/p2zsed.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P2Z/p2zsed.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P2Z/p2zsed.F90	(revision 12680)
@@ -33,4 +33,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zbc.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zbc.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zbc.F90	(revision 12680)
@@ -48,4 +48,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zflx.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zflx.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zflx.F90	(revision 12680)
@@ -54,4 +54,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zopt.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zopt.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zopt.F90	(revision 12680)
@@ -44,4 +44,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zpoc.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zpoc.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zpoc.F90	(revision 12680)
@@ -39,4 +39,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zrem.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zrem.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zrem.F90	(revision 12680)
@@ -44,4 +44,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zsed.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zsed.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zsed.F90	(revision 12680)
@@ -39,4 +39,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zsms.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 12680)
@@ -41,4 +41,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/SED/oce_sed.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/SED/oce_sed.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/SED/oce_sed.F90	(revision 12680)
@@ -13,5 +13,5 @@
    USE dom_oce , ONLY :   glamt     =>   glamt          !: longitude of t-point (degre)
    USE dom_oce , ONLY :   gphit     =>   gphit          !: latitude  of t-point (degre)
-   USE dom_oce , ONLY :   e3t       =>   e3t            !: latitude  of t-point (degre)
+!!st USE dom_oce , ONLY :   e3t       =>   e3t            !: latitude  of t-point (degre)
    USE dom_oce , ONLY :   e3t_1d    =>   e3t_1d         !: reference depth of t-points (m)
    USE dom_oce , ONLY :   gdepw_0   =>   gdepw_0        !: reference depth of t-points (m)
@@ -53,4 +53,2 @@
 
 END MODULE oce_sed
-
-
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/SED/seddta.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/SED/seddta.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/SED/seddta.F90	(revision 12680)
@@ -24,4 +24,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !! $Id$
 CONTAINS
@@ -164,5 +165,5 @@
       CALL pack_arr ( jpoce,  rainrm_dta(1:jpoce,jscal), trc_data(1:jpi,1:jpj,14), iarroce(1:jpoce) )
       rainrm_dta(1:jpoce,jscal) = rainrm_dta(1:jpoce,jscal) * 1e-4
-      ! vector temperature [°C] and salinity 
+      ! vector temperature [�C] and salinity 
       CALL pack_arr ( jpoce,  temp(1:jpoce), trc_data(1:jpi,1:jpj,15), iarroce(1:jpoce) )
       CALL pack_arr ( jpoce,  salt(1:jpoce), trc_data(1:jpi,1:jpj,16), iarroce(1:jpoce) )
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/TRP/trcsink.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/TRP/trcsink.F90	(revision 12679)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/TRP/trcsink.F90	(revision 12680)
@@ -26,4 +26,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
