Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DOM/domain.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DOM/domain.F90	(revision 14428)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DOM/domain.F90	(revision 14547)
@@ -91,5 +91,4 @@
       INTEGER ::   ji, jj, jk, jt   ! dummy loop indices
       INTEGER ::   iconf = 0    ! local integers
-      REAL(wp)::   zrdt
       CHARACTER (len=64) ::   cform = "(A12, 3(A13, I7))"
       INTEGER , DIMENSION(jpi,jpj) ::   ik_top , ik_bot       ! top and bottom ocean level
@@ -372,6 +371,25 @@
       !
       ! set current model timestep rDt = 2*rn_Dt if MLF or rDt = rn_Dt if RK3
+#if defined key_RK3
+      rDt   =         rn_Dt
+      r1_Dt = 1._wp / rDt
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '           ===>>>   Runge Kutta 3rd order (RK3) :   rDt = ', rDt
+         WRITE(numout,*)
+      ENDIF
+      !
+#else
       rDt   = 2._wp * rn_Dt
       r1_Dt = 1._wp / rDt
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '           ===>>>   Modified Leap-Frog (MLF) :   rDt = ', rDt
+         WRITE(numout,*)
+      ENDIF
+      !
+#endif
       !
       IF( l_SAS .AND. .NOT.ln_linssh ) THEN
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/divhor.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/divhor.F90	(revision 14428)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/divhor.F90	(revision 14547)
@@ -13,4 +13,5 @@
    !!             -   ! 2014-12  (G. Madec) suppression of cross land advection option
    !!             -   ! 2015-10  (G. Madec) add velocity and rnf flag in argument of div_hor
+   !!            4.5  ! 2015-10  (S. Techene, G. Madec)  hdiv replaced by e3divUh
    !!----------------------------------------------------------------------
 
@@ -36,5 +37,10 @@
    PRIVATE
 
-   PUBLIC   div_hor    ! routine called by step.F90 and istate.F90
+   !                  !! * Interface
+   INTERFACE div_hor
+      MODULE PROCEDURE div_hor_RK3, div_hor_old
+   END INTERFACE
+
+   PUBLIC   div_hor    ! routine called by ssh_nxt.F90 and istate.F90
 
    !! * Substitutions
@@ -48,5 +54,65 @@
 CONTAINS
 
-   SUBROUTINE div_hor( kt, Kbb, Kmm )
+   SUBROUTINE div_hor_RK3( kt, Kbb, Kmm, puu, pvv, pe3divUh )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE div_hor_RK3  ***
+      !!                    
+      !! ** Purpose :   compute the horizontal divergence at now time-step
+      !!
+      !! ** Method  :   the now divergence is computed as :
+      !!         hdiv = 1/(e1e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
+      !!      and correct with runoff inflow (div_rnf) and cross land flow (div_cla) 
+      !!
+      !! ** Action  : - thickness weighted horizontal divergence of in input velocity (puu,pvv)
+      !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in   ) ::   kt, Kbb, Kmm   ! ocean time-step & time-level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   puu, pvv       ! horizontal velocity
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pe3divUh       ! e3t*div[Uh]
+      !
+      INTEGER  ::   ji, jj, jk    ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('div_hor_RK3')
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'div_hor_RK3 : thickness weighted horizontal divergence '
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+         hdiv    (:,:,:) = 0._wp    ! initialize hdiv & pe3divUh for the halos and jpk level at the first time step
+      ENDIF
+      ! 
+      pe3divUh(:,:,:) = 0._wp    !!gm to be applied to the halos only
+      !
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )                                                  
+         hdiv(ji,jj,jk) = (   e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * puu(ji  ,jj,jk)      &
+            &               - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * puu(ji-1,jj,jk)      &
+            &               + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * pvv(ji,jj  ,jk)      &
+            &               - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * pvv(ji,jj-1,jk)  )   &
+            &            * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+      END_3D
+      !
+      IF( ln_rnf )   CALL sbc_rnf_div( hdiv, Kmm )             !==  + runoffs divergence  ==!
+      !
+#if defined key_asminc 
+      IF( ln_sshinc .AND. ln_asmiau )   &                      !==  + SSH assimilation increment  ==!
+         &           CALL ssh_asm_div( kt, Kbb, Kmm, hdiv )
+#endif
+      !
+      IF( ln_isf )   CALL isf_hdiv( kt, Kmm, hdiv )            !==  + ice-shelf mass exchange ==!
+      !
+!!gm Patch before suppression of hdiv from all modules that use it
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )                            !==  e3t * Horizontal divergence  ==!
+         pe3divUh(ji,jj,jk) = hdiv(ji,jj,jk) * e3t(ji,jj,jk,Kmm)
+      END_3D
+!!gm end
+      !
+      CALL lbc_lnk( 'divhor', hdiv, 'T', 1._wp )   !   (no sign change)
+      !
+      IF( ln_timing )   CALL timing_stop('div_hor_RK3')
+      !
+   END SUBROUTINE div_hor_RK3
+
+
+   SUBROUTINE div_hor_old( kt, Kbb, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE div_hor  ***
@@ -64,6 +130,4 @@
       !
       INTEGER  ::   ji, jj, jk    ! dummy loop indices
-      REAL(wp) ::   zraur, zdep   ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj) :: ztmp
       !!----------------------------------------------------------------------
       !
@@ -98,5 +162,5 @@
       IF( ln_timing )   CALL timing_stop('div_hor')
       !
-   END SUBROUTINE div_hor
+   END SUBROUTINE div_hor_old
    
    !!======================================================================
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynspg.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynspg.F90	(revision 14428)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynspg.F90	(revision 14547)
@@ -57,5 +57,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_spg( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa, k_only_ADV )
+   SUBROUTINE dyn_spg( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_spg  ***
@@ -79,8 +79,7 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation
       REAL(wp), DIMENSION(jpi,jpj,jpt)    , INTENT(inout) ::  pssh, puu_b, pvv_b  ! SSH and barotropic velocities at main time levels
-      INTEGER , OPTIONAL                  , INTENT( in )  ::  k_only_ADV          ! only Advection in the RHS
       !
       INTEGER  ::   ji, jj, jk                   ! dummy loop indices
-      REAL(wp) ::   z2dt, zg_2, zintp, zgrho0r, zld   ! local scalars
+      REAL(wp) ::   zg_2, zintp, zgrho0r, zld    ! local scalars
       REAL(wp)             , DIMENSION(jpi,jpj) ::   zpgu, zpgv   ! 2D workspace
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zpice
@@ -151,6 +150,6 @@
          IF( ln_wave .and. ln_bern_srfc ) THEN          !== Add J terms: depth-independent Bernoulli head
             DO_2D( 0, 0, 0, 0 )
-               zpgu(ji,jj) = zpgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) / e1u(ji,jj)   !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3]
-               zpgv(ji,jj) = zpgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) / e2v(ji,jj)
+               zpgu(ji,jj) = zpgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) * r1_e1u(ji,jj)   !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3]
+               zpgv(ji,jj) = zpgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) * r1_e2v(ji,jj)
             END_2D
          ENDIF
@@ -167,5 +166,5 @@
       SELECT CASE ( nspg )                   !== surface pressure gradient computed and add to the general trend ==!
       CASE ( np_EXP )   ;   CALL dyn_spg_exp( kt,      Kmm,       puu, pvv, Krhs )                    ! explicit
-      CASE ( np_TS  )   ;   CALL dyn_spg_ts ( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa, k_only_ADV ) ! time-splitting
+      CASE ( np_TS  )   ;   CALL dyn_spg_ts ( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa ) ! time-splitting
       END SELECT
       !
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynspg_ts.F90	(revision 14428)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynspg_ts.F90	(revision 14547)
@@ -68,4 +68,5 @@
    PUBLIC dyn_spg_ts        ! called by dyn_spg 
    PUBLIC dyn_spg_ts_init   !    -    - dyn_spg_init
+   PUBLIC dyn_drg_init      ! called by rk3ssh
 
    !! Time filtered arrays at baroclinic time step:
@@ -80,4 +81,8 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ftsw, ftse         ! (only used with een vorticity scheme)
 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshe_rhs           ! RHS of ssh Eq.
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ue_rhs, Ve_rhs    ! RHS of barotropic velocity Eq.
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   CdU_u, CdU_v       ! top/bottom stress at u- & v-points
+
    REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! local ratios
    REAL(wp) ::   r1_8  = 0.125_wp         !
@@ -117,5 +122,5 @@
 
 
-   SUBROUTINE dyn_spg_ts( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa, k_only_ADV )
+   SUBROUTINE dyn_spg_ts( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa )
       !!----------------------------------------------------------------------
       !!
@@ -145,6 +150,5 @@
       INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs, Kaa ! ocean time level indices
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation
-      REAL(wp), DIMENSION(jpi,jpj,jpt)    , INTENT(inout) ::  pssh, puu_b, pvv_b  ! SSH and barotropic velocities at main time levels
-      INTEGER , OPTIONAL                  , INTENT( in )  ::  k_only_ADV          ! only Advection in the RHS
+      REAL(wp), DIMENSION(jpi,jpj    ,jpt), INTENT(inout) ::  pssh, puu_b, pvv_b  ! SSH and barotropic velocities at main time levels
       !
       INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices
@@ -152,9 +156,9 @@
       LOGICAL  ::   ll_init               ! =T : special startup of 2d equations
       INTEGER  ::   noffset               ! local integers  : time offset for bdy update
-      REAL(wp) ::   r1_Dt_b, z1_hu, z1_hv          ! local scalars
-      REAL(wp) ::   za0, za1, za2, za3              !   -      -
+      REAL(wp) ::   z1_hu, z1_hv              ! local scalars
+      REAL(wp) ::   za0, za1, za2, za3        !   -      -
       REAL(wp) ::   zztmp, zldg               !   -      -
-      REAL(wp) ::   zhu_bck, zhv_bck, zhdiv         !   -      -
-      REAL(wp) ::   zun_save, zvn_save              !   -      -
+      REAL(wp) ::   zhu_bck, zhv_bck, zhdiv   !   -      -
+      REAL(wp) ::   zun_save, zvn_save        !   -      -
       REAL(wp), DIMENSION(jpi,jpj) :: zu_trd, zu_frc, zu_spg, zssh_frc
       REAL(wp), DIMENSION(jpi,jpj) :: zv_trd, zv_frc, zv_spg
@@ -163,7 +167,4 @@
       REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points
       REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV         ! fluxes
-!!st#if defined key_qco 
-!!st      REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v
-!!st#endif
       !
       REAL(wp) ::   zwdramp                     ! local scalar - only used if ln_wd_dl = .True. 
@@ -175,4 +176,5 @@
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2    ! averages over the sub-steps of zuwdmask and zvwdmask
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2d          ! 2D workspace
       REAL(wp) ::   zt0substep !   Time of day at the beginning of the time substep
       !!----------------------------------------------------------------------
@@ -185,5 +187,4 @@
 !     zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp
       !                                         ! inverse of baroclinic time step 
-      r1_Dt_b = 1._wp / rDt 
       !
       ll_init     = ln_bt_av                    ! if no time averaging, then no specific restart 
@@ -228,15 +229,42 @@
       !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
       ! -----------------------------------------------------------------------------
+#if defined key_RK3
+      !                    !========================================!
+      !                    !==  Phase 1 for RK3 time integration  ==!
+      !                    !========================================!
+      !
+      !                          ! Currently, RK3 requires the forward mode and NO time filtering of barotropic variables
+      IF( kt == nit000 ) THEN
+         IF( .NOT.ln_bt_fw  .OR. ln_bt_av )   CALL ctl_stop( 'dyn_spg_ts: RK3 requires ln_bt_fw=T AND ln_bt_av=F')
+      ENDIF
+      !
+      !                          ! set values computed in RK3_ssh
+      zssh_frc(:,:) = sshe_rhs(:,:)
+        zu_frc(:,:) =   Ue_rhs(:,:)
+        zv_frc(:,:) =   Ve_rhs(:,:)
+      zCdU_u  (:,:) = CdU_u   (:,:)
+      zCdU_v  (:,:) = CdU_v   (:,:)
+
+!!gm ==>>> !!ts    ISSUe her on en discute    changer les cas ENS ENE  et triad ?
+      IF( kt == nit000 .AND. .NOT. ln_linssh )   CALL dyn_cor_2D_init( Kmm )   ! Set zwz, the barotropic Coriolis force coefficient
+      !                      ! recompute zwz = f/depth  at every time step for (.NOT.ln_linssh) as the water colomn height changes
+      !
+      
+#else
+      !                    !========================================!
+      !                    !==  Phase 1 for MLF time integration  ==!
+      !                    !========================================!
+      !
       !      
       !
       !                                   !=  zu_frc =  1/H e3*d/dt(Ua)  =!  (Vertical mean of Ua, the 3D trends)
       !                                   !  ---------------------------  !
-#if defined key_qco 
+# if defined key_qco 
       zu_frc(:,:) = SUM( e3u_0(:,:,:  ) * puu(:,:,:,Krhs) * umask(:,:,:), DIM=3 ) * r1_hu_0(:,:)
       zv_frc(:,:) = SUM( e3v_0(:,:,:  ) * pvv(:,:,:,Krhs) * vmask(:,:,:), DIM=3 ) * r1_hv_0(:,:)
-#else
+# else
       zu_frc(:,:) = SUM( e3u(:,:,:,Kmm) * puu(:,:,:,Krhs) * umask(:,:,:), DIM=3 ) * r1_hu(:,:,Kmm)
       zv_frc(:,:) = SUM( e3v(:,:,:,Kmm) * pvv(:,:,:,Krhs) * vmask(:,:,:), DIM=3 ) * r1_hv(:,:,Kmm)
-#endif 
+# endif 
       !
       !
@@ -256,27 +284,19 @@
       !                      ! recompute zwz = f/depth  at every time step for (.NOT.ln_linssh) as the water colomn height changes
       !
-      IF( .NOT. PRESENT(k_only_ADV) ) THEN   !* remove the 2D Coriolis trend  
-         zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:)        ! now fluxes 
-         zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:)        ! NB: FULL domain : put a value in last row and column
-         !
-         CALL dyn_cor_2d( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV,  &   ! <<== in
-            &                                                                          zu_trd, zv_trd   )   ! ==>> out
-         !
-         DO_2D( 0, 0, 0, 0 )                          ! Remove coriolis term (and possibly spg) from barotropic trend
-            zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
-            zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
-         END_2D
-      ENDIF
+      zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:)        ! now fluxes 
+      zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:)        ! NB: FULL domain : put a value in last row and column
+      !
+      CALL dyn_cor_2d( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV,  &   ! <<== in
+         &                                                                          zu_trd, zv_trd   )   ! ==>> out
+      !
+      DO_2D( 0, 0, 0, 0 )                          ! Remove coriolis term (and possibly spg) from barotropic trend
+         zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
+         zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
+      END_2D
       !
       !                                   !=  Add bottom stress contribution from baroclinic velocities  =!
       !                                   !  -----------------------------------------------------------  !
-      IF( PRESENT(k_only_ADV) ) THEN         !* only Advection in the RHS : provide the barotropic bottom drag coefficients
-         DO_2D( 0, 0, 0, 0 )
-            zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
-            zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
-         END_2D
-      ELSE				     !* remove baroclinic drag AND provide the barotropic drag coefficients
-         CALL dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b, pvv_b, zu_frc, zv_frc, zCdU_u, zCdU_v )
-      ENDIF
+      CALL dyn_drg_init( Kbb, Kmm, puu   , pvv   , puu_b , pvv_b  ,   &     !  <<= IN
+         &                         zu_frc, zv_frc, zCdU_u, zCdU_v )         !  =>> OUT
       !
       !                                   !=  Add atmospheric pressure forcing  =!
@@ -348,5 +368,5 @@
       END IF
       !
-#if defined key_asminc
+# if defined key_asminc
       !                                   !=  Add the IAU weighted SSH increment  =!
       !                                   !  ------------------------------------  !
@@ -354,5 +374,10 @@
          zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
       ENDIF
+# endif
+
+      !                    !==  END of  Phase 1 for MLF time integration  ==!
 #endif
+
+
       !                                   != Fill boundary data arrays for AGRIF
       !                                   ! ------------------------------------
@@ -701,32 +726,63 @@
          sshb_e (:,:) = sshn_e(:,:)
          sshn_e (:,:) = ssha_e(:,:)
-
-         !                                             !* Sum over whole bt loop
+         !
+         !                                             !* Sum over whole bt loop (except in weight average)
          !                                             !  ----------------------
-         za1 = wgtbtp1(jn)                                    
-         IF( ln_dynadv_vec .OR. ln_linssh ) THEN    ! Sum velocities
-            puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) 
-            pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) 
-         ELSE                                       ! Sum transports
-            IF ( .NOT.ln_wd_dl ) THEN  
-               puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:)
-               pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:)
-            ELSE 
-               puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:) * zuwdmask(:,:)
-               pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:) * zvwdmask(:,:)
-            END IF 
-         ENDIF
-         !                                          ! Sum sea level
-         pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:)
-
+         IF( ln_bt_av ) THEN
+            za1 = wgtbtp1(jn)                                    
+            IF( ln_dynadv_vec .OR. ln_linssh ) THEN    ! Sum velocities
+               puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) 
+               pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) 
+            ELSE                                       ! Sum transports
+               IF ( .NOT.ln_wd_dl ) THEN  
+                  puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:)
+                  pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:)
+               ELSE 
+                  puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:) * zuwdmask(:,:)
+                  pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:) * zvwdmask(:,:)
+               ENDIF
+            ENDIF
+            !                                          ! Sum sea level
+            pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:)
+         ENDIF
+         !
          !                                                 ! ==================== !
       END DO                                               !        end loop      !
       !                                                    ! ==================== !
+
+      
       ! -----------------------------------------------------------------------------
       ! Phase 3. update the general trend with the barotropic trend
       ! -----------------------------------------------------------------------------
       !
-      ! Set advection velocity correction:
-      IF (ln_bt_fw) THEN
+      IF(.NOT.ln_bt_av ) THEN                          !* Update Kaa barotropic external mode 
+         puu_b(:,:,Kaa) = ua_e  (:,:)
+         pvv_b(:,:,Kaa) = va_e  (:,:)
+         pssh (:,:,Kaa) = ssha_e(:,:)
+      ENDIF
+
+#if defined key_RK3
+      !                                                !*  RK3 case
+      !
+      ! advective transport from N to N+1 (i.e. Kbb to Kaa) 
+      ub2_b(:,:) = un_adv(:,:)      ! Save integrated transport for next computation (NOT USED)
+      vb2_b(:,:) = vn_adv(:,:)
+      ! 
+      IF( iom_use("ubar") ) THEN    ! RK3 single first: hu[N+1/2] = 1/2 ( hu[N] + hu[N+1] ) 
+         ALLOCATE( z2d(jpi,jpj) )
+         z2d(:,:) = 2._wp / ( hu_e(:,:) + hu(:,:,Kbb) + 1._wp - ssumask(:,:) ) 
+         CALL iom_put(  "ubar", un_adv(:,:)*z2d(:,:) )    ! barotropic i-current
+         z2d(:,:) = 2._wp / ( hv_e(:,:) + hv(:,:,Kbb) + 1._wp - ssvmask(:,:) )
+         CALL iom_put(  "vbar", vn_adv(:,:)*z2d(:,:) )    ! barotropic i-current
+         DEALLOCATE( z2d )
+      ENDIF
+      !
+      !                    !==  END Phase 3 for RK3 (forward mode) ==!
+
+#else
+      !                                                !*  MLF case
+      !
+      ! Set advective velocity correction:
+      IF( ln_bt_fw ) THEN
          IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN
             DO_2D( 1, 1, 1, 1 )
@@ -748,42 +804,51 @@
             ub2_b(:,:) = un_adv(:,:)      ! Save integrated transport for next computation
             vb2_b(:,:) = vn_adv(:,:)
-         END IF
-      ENDIF
-
-
+         ENDIF
+      ENDIF
       !
       ! Update barotropic trend:
       IF( ln_dynadv_vec .OR. ln_linssh ) THEN
          DO jk=1,jpkm1
-            puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_Dt_b
-            pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_Dt_b
+            puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_Dt
+            pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_Dt
          END DO
       ELSE
-         ! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points
-#if defined key_qcoTest_FluxForm
-         !                                ! 'key_qcoTest_FluxForm' : simple ssh average
-         DO_2D( 1, 0, 1, 0 )
-            zsshu_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji+1,jj  ,Kaa) ) * ssumask(ji,jj)
-            zsshv_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji  ,jj+1,Kaa) ) * ssvmask(ji,jj)
-         END_2D
-#else
-         DO_2D( 1, 0, 1, 0 )
-            zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji  ,jj) * pssh(ji  ,jj,Kaa)   &
-               &                                      + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) * ssumask(ji,jj)
-            zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj  ) * pssh(ji,jj  ,Kaa)   &
-               &                                      + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj)
-         END_2D
-#endif   
-         CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
-         !
-         DO jk=1,jpkm1
-            puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm)   &
-               &             * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b
-            pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm)   &
-               &             * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b
-         END DO
-         ! Save barotropic velocities not transport:
-         puu_b(:,:,Kaa) =  puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
-         pvv_b(:,:,Kaa) =  pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
+         IF(.NOT.ln_bt_av ) THEN   ! (puu_b,pvv_b)_Kaa is a velocity (hu,hv)_Kaa = (hu_e,hv_e)
+            ! 
+            DO jk=1,jpkm1
+               puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm)   &
+                  &             * ( puu_b(:,:,Kaa)*hu_e(:,:) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt
+               pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm)   &
+                  &             * ( pvv_b(:,:,Kaa)*hv_e(:,:) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt
+            END DO
+            !
+         ELSE                  ! at this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points
+            !
+# if defined key_qcoTest_FluxForm
+            !                                ! 'key_qcoTest_FluxForm' : simple ssh average
+            DO_2D( 1, 0, 1, 0 )
+               zsshu_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji+1,jj  ,Kaa) ) * ssumask(ji,jj)
+               zsshv_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji  ,jj+1,Kaa) ) * ssvmask(ji,jj)
+            END_2D
+# else
+            DO_2D( 1, 0, 1, 0 )
+               zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji  ,jj) * pssh(ji  ,jj,Kaa)   &
+                  &                                      + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) * ssumask(ji,jj)
+               zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj  ) * pssh(ji,jj  ,Kaa)   &
+                  &                                      + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj)
+            END_2D
+# endif   
+            CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
+            !
+            DO jk=1,jpkm1
+               puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm)   &
+                  &             * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt
+               pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm)   &
+                  &             * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt
+            END DO
+            ! Save barotropic velocities not transport:
+            puu_b(:,:,Kaa) =  puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
+            pvv_b(:,:,Kaa) =  pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
+         ENDIF
       ENDIF
 
@@ -795,16 +860,19 @@
       END DO
 
-      IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN 
+      IF( ln_wd_dl .AND. ln_wd_dl_bc) THEN 
          DO jk = 1, jpkm1
             puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) &
-                       & + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk) 
+               &            + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk) 
             pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) & 
-                       & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk)  
+               &            + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk) 
          END DO
-      END IF 
-
-      
+      ENDIF
+            
       CALL iom_put(  "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) )    ! barotropic i-current
       CALL iom_put(  "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) )    ! barotropic i-current
+
+      !                    !==  END Phase 3 for MLF time integration  ==!
+#endif
+
       !
 #if defined key_agrif
@@ -834,6 +902,6 @@
    END SUBROUTINE dyn_spg_ts
 
-
-   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
+   
+   SUBROUTINE ts_wgt( ll_av, ll_fw, Kpit, zwgt1, zwgt2)
       !!---------------------------------------------------------------------
       !!                   ***  ROUTINE ts_wgt  ***
@@ -841,71 +909,71 @@
       !! ** Purpose : Set time-splitting weights for temporal averaging (or not)
       !!----------------------------------------------------------------------
-      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
-      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
-      INTEGER, INTENT(inout) :: jpit      ! cycle length    
-      REAL(wp), DIMENSION(3*nn_e), INTENT(inout) ::   zwgt1, & ! Primary weights
-                                                         zwgt2    ! Secondary weights
-      
-      INTEGER ::  jic, jn, ji                      ! temporary integers
-      REAL(wp) :: za1, za2
-      !!----------------------------------------------------------------------
-
+      LOGICAL, INTENT(in   ) ::   ll_av     ! temporal averaging=.true.
+      LOGICAL, INTENT(in   ) ::   ll_fw     ! forward time splitting =.true.
+      INTEGER, INTENT(inout) ::   Kpit      ! cycle length
+      !!
+      INTEGER ::  jic, jn, ji   ! local integers
+      REAL(wp) :: za1, za2      ! loca scalars
+      REAL(wp), DIMENSION(3*nn_e), INTENT(inout) ::   zwgt1, zwgt2   ! Primary & Secondary weights
+      !!----------------------------------------------------------------------
+      !
       zwgt1(:) = 0._wp
       zwgt2(:) = 0._wp
-
-      ! Set time index when averaged value is requested
-      IF (ll_fw) THEN 
-         jic = nn_e
-      ELSE
-         jic = 2 * nn_e
-      ENDIF
-
-      ! Set primary weights:
-      IF (ll_av) THEN
-           ! Define simple boxcar window for primary weights 
-           ! (width = nn_e, centered around jic)     
-         SELECT CASE ( nn_bt_flt )
-              CASE( 0 )  ! No averaging
-                 zwgt1(jic) = 1._wp
-                 jpit = jic
-
-              CASE( 1 )  ! Boxcar, width = nn_e
-                 DO jn = 1, 3*nn_e
-                    za1 = ABS(float(jn-jic))/float(nn_e) 
-                    IF (za1 < 0.5_wp) THEN
-                      zwgt1(jn) = 1._wp
-                      jpit = jn
-                    ENDIF
-                 ENDDO
-
-              CASE( 2 )  ! Boxcar, width = 2 * nn_e
-                 DO jn = 1, 3*nn_e
-                    za1 = ABS(float(jn-jic))/float(nn_e) 
-                    IF (za1 < 1._wp) THEN
-                      zwgt1(jn) = 1._wp
-                      jpit = jn
-                    ENDIF
-                 ENDDO
-              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
+      !
+      !                          !==  Set time index when averaged value is requested  ==!
+      IF (ll_fw) THEN   ;   jic =     nn_e
+      ELSE              ;   jic = 2 * nn_e
+      ENDIF
+      !
+      !                          !==  Set primary weights  ==!
+      !
+      IF (ll_av) THEN               != Define simple boxcar window for primary weights 
+         !                                       ! (width = nn_e, centered around jic)     
+         SELECT CASE( nn_bt_flt )
+         !
+         CASE( 0 )                  ! No averaging
+            zwgt1(jic) = 1._wp
+            Kpit = jic
+            !
+         CASE( 1 )                  ! Boxcar, width = nn_e
+            DO jn = 1, 3*nn_e
+               za1 = ABS( REAL( jn-jic, wp) ) / REAL( nn_e, wp ) 
+               IF( za1 < 0.5_wp ) THEN
+                  zwgt1(jn) = 1._wp
+                  Kpit = jn
+               ENDIF
+            END DO
+            !
+         CASE( 2 )                  ! Boxcar, width = 2 * nn_e
+            DO jn = 1, 3*nn_e
+               za1 = ABS(REAL( jn-jic, wp) ) / REAL( nn_e, wp ) 
+               IF( za1 < 1._wp ) THEN
+                  zwgt1(jn) = 1._wp
+                  Kpit = jn
+               ENDIF
+            END DO
+            !
+         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
+         !
          END SELECT
 
-      ELSE ! No time averaging
+      ELSE                          !=  No time averaging
          zwgt1(jic) = 1._wp
-         jpit = jic
-      ENDIF
-    
-      ! Set secondary weights
-      DO jn = 1, jpit
-        DO ji = jn, jpit
-             zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
-        END DO
+         Kpit = jic
+      ENDIF
+      !
+      !                          !==  Set secondary weights  ==!
+      DO jn = 1, Kpit
+         DO ji = jn, Kpit
+            zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
+         END DO
       END DO
-
-      ! Normalize weigths:
-      za1 = 1._wp / SUM(zwgt1(1:jpit))
-      za2 = 1._wp / SUM(zwgt2(1:jpit))
-      DO jn = 1, jpit
-        zwgt1(jn) = zwgt1(jn) * za1
-        zwgt2(jn) = zwgt2(jn) * za2
+      !
+      !                          !==  Normalize weigths  ==!
+      za1 = 1._wp / SUM( zwgt1(1:Kpit) )
+      za2 = 1._wp / SUM( zwgt2(1:Kpit) )
+      DO jn = 1, Kpit
+         zwgt1(jn) = zwgt1(jn) * za1
+         zwgt2(jn) = zwgt2(jn) * za2
       END DO
       !
@@ -1023,5 +1091,5 @@
          IF(lwp) WRITE(numout,*) '     ln_ts_auto=.false.: Use nn_e in namelist   nn_e = ', nn_e
       ENDIF
-
+      !
       IF(ln_bt_av) THEN
          IF(lwp) WRITE(numout,*) '     ln_bt_av =.true.  ==> Time averaging over nn_e time steps is on '
@@ -1029,5 +1097,4 @@
          IF(lwp) WRITE(numout,*) '     ln_bt_av =.false. => No time averaging of barotropic variables '
       ENDIF
-      !
       !
       IF(ln_bt_fw) THEN
@@ -1088,5 +1155,6 @@
       !! and update depths at T- points (ht) at each barotropic time step
       !!
-      !! Compute zwz = f / ( height of the water colomn )
+      !! Compute zwz = f/h              (potential planetary voricity)
+      !! Compute ftne, ftnw, ftse, ftsw (triad of potential planetary voricity)
       !!----------------------------------------------------------------------
       INTEGER,  INTENT(in)         ::  Kmm  ! Time index
@@ -1138,5 +1206,5 @@
          !
       END SELECT
-      
+      !
    END SUBROUTINE dyn_cor_2d_init
 
@@ -1287,9 +1355,11 @@
       !! ** Purpose : 
       !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pshn
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
+      !
       INTEGER  ::   ji ,jj               ! dummy loop indices
       LOGICAL  ::   ll_tmp1, ll_tmp2
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pshn
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
-      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      !
       DO_2D( 0, 0, 0, 0 )
          ll_tmp1 = MIN(  pshn(ji,jj)               ,  pshn(ji+1,jj) ) >                &
@@ -1330,5 +1400,5 @@
          ENDIF
       END_2D
-            
+      !
    END SUBROUTINE wad_spg
      
@@ -1374,5 +1444,4 @@
       !
       IF( ln_bt_fw ) THEN                 ! FORWARD integration: use NOW bottom baroclinic velocities
-         
          DO_2D( 0, 0, 0, 0 )
             ikbu = mbku(ji,jj)       
@@ -1381,6 +1450,5 @@
             zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm)
          END_2D
-      ELSE                                ! CENTRED integration: use BEFORE bottom baroclinic velocities
-         
+      ELSE                                ! CENTRED integration: use BEFORE bottom baroclinic velocities        
          DO_2D( 0, 0, 0, 0 )
             ikbu = mbku(ji,jj)       
@@ -1400,5 +1468,4 @@
          END_2D
       ELSE                    ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
-         
          DO_2D( 0, 0, 0, 0 )
             pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj)
@@ -1412,5 +1479,4 @@
          !
          IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity
-            
             DO_2D( 0, 0, 0, 0 )
                iktu = miku(ji,jj)
@@ -1419,6 +1485,5 @@
                zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm)
             END_2D
-         ELSE                                ! CENTRED integration: use BEFORE top baroclinic velocity
-            
+         ELSE                               ! CENTRED integration: use BEFORE top baroclinic velocity
             DO_2D( 0, 0, 0, 0 )
                iktu = miku(ji,jj)
@@ -1428,7 +1493,5 @@
             END_2D
          ENDIF
-         !
-         !                    ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
-         
+         !                    ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)       
          DO_2D( 0, 0, 0, 0 )
             pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj)
@@ -1439,4 +1502,5 @@
       !
    END SUBROUTINE dyn_drg_init
+
 
    SUBROUTINE ts_bck_interp( jn, ll_init,       &   ! <== in
@@ -1477,5 +1541,4 @@
    END SUBROUTINE ts_bck_interp
 
-
    !!======================================================================
 END MODULE dynspg_ts
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynvor.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynvor.F90	(revision 14428)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynvor.F90	(revision 14547)
@@ -51,4 +51,8 @@
    IMPLICIT NONE
    PRIVATE
+   
+   INTERFACE dyn_vor
+      MODULE PROCEDURE dyn_vor_3D, dyn_vor_RK3
+   END INTERFACE
 
    PUBLIC   dyn_vor        ! routine called by step.F90
@@ -74,6 +78,7 @@
    INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme
 
-   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity
-   !                               ! associated indices:
+   !                                    !: choice of calculated vorticity
+   INTEGER, PUBLIC ::   ncor, nrvm, ntot   ! Coriolis, relative vorticity, total vorticity
+   !                                       ! associated indices:
    INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary)
    INTEGER, PUBLIC, PARAMETER ::   np_RVO = 2         ! relative vorticity
@@ -103,6 +108,6 @@
    !!----------------------------------------------------------------------
 CONTAINS
-
-   SUBROUTINE dyn_vor( kt, Kmm, puu, pvv, Krhs )
+   
+   SUBROUTINE dyn_vor_3D( kt, Kmm, puu, pvv, Krhs )
       !!----------------------------------------------------------------------
       !!
@@ -121,5 +126,5 @@
       !!----------------------------------------------------------------------
       !
-      IF( ln_timing )   CALL timing_start('dyn_vor')
+      IF( ln_timing )   CALL timing_start('dyn_vor_3D')
       !
       IF( l_trddyn ) THEN     !==  trend diagnostics case : split the added trend in two parts  ==!
@@ -209,7 +214,83 @@
          &                                tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
-      IF( ln_timing )   CALL timing_stop('dyn_vor')
-      !
-   END SUBROUTINE dyn_vor
+      IF( ln_timing )   CALL timing_stop('dyn_vor_3D')
+      !
+   END SUBROUTINE dyn_vor_3D
+
+
+   SUBROUTINE dyn_vor_RK3( kt, Kmm, puu, pvv, Krhs, knoco )
+      !!----------------------------------------------------------------------
+      !!
+      !! ** Purpose :   compute the lateral ocean tracer physics.
+      !!
+      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend
+      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
+      !!               and planetary vorticity trends) and send them to trd_dyn
+      !!               for futher diagnostics (l_trddyn=T)
+      !!----------------------------------------------------------------------
+      INTEGER                             , INTENT(in   ) ::   kt          ! ocean time-step index
+      INTEGER                             , INTENT(in   ) ::   Kmm, Krhs   ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv    ! ocean velocity field and RHS of momentum equation
+      INTEGER                             , INTENT(in   ) ::   knoco       ! specified vorticity trend (= np_MET or np_RVO)
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('dyn_vor_RK3')
+      !
+      !              !==  total vorticity trend added to the general trend  ==!
+      !!st WARNING 22/02 !!!!!!!! stoke drift or not stoke drift ? => bar to do later !!!
+      !! stoke drift a garder pas vortex force a priori !!
+      !! ATTENTION déja appelé avec Kbb !!
+      
+         !
+         SELECT CASE ( nvor_scheme )      !==  vorticity trend added to the general trend  ==!
+         CASE( np_ENT )                        !* energy conserving scheme  (T-pts)
+                             CALL vor_enT( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
+            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
+                             CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
+                             CALL vor_enT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force
+            ENDIF
+         CASE( np_EET )                        !* energy conserving scheme (een scheme using e3t)
+                             CALL vor_eeT( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
+            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
+                             CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
+                             CALL vor_eeT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force
+            ENDIF
+         CASE( np_ENE )                        !* energy conserving scheme
+                             CALL vor_ene( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
+            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
+                             CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
+                             CALL vor_ene( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force
+            ENDIF
+         CASE( np_ENS )                        !* enstrophy conserving scheme
+                             CALL vor_ens( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend
+
+            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
+                             CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend
+            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
+                             CALL vor_ens( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend and vortex force
+            ENDIF
+         CASE( np_MIX )                        !* mixed ene-ens scheme
+                             CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! relative vorticity or metric trend (ens)
+            IF( ln_stcor )        CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )        ! add the Stokes-Coriolis trend
+            IF( ln_vortex_force ) CALL vor_ens( kt, Kmm, nrvm, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add vortex force
+         CASE( np_EEN )                        !* energy and enstrophy conserving scheme
+                             CALL vor_een( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
+            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
+                             CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
+                             CALL vor_een( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force
+            ENDIF
+         END SELECT
+      !
+      !                       ! print sum trends (used for debugging)
+      IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' vor  - Ua: ', mask1=umask,               &
+         &                                tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      !
+      IF( ln_timing )   CALL timing_stop('dyn_vor_RK3')
+      !
+   END SUBROUTINE dyn_vor_RK3
 
 
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynzdf.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynzdf.F90	(revision 14428)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynzdf.F90	(revision 14547)
@@ -70,12 +70,12 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation
       !
-      INTEGER  ::   ji, jj, jk         ! dummy loop indices
-      INTEGER  ::   iku, ikv           ! local integers
-      REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars
-      REAL(wp) ::   zzws, ze3va        !   -      -
-      REAL(wp) ::   z1_e3ua, z1_e3va   !   -      -
-      REAL(wp) ::   zWu , zWv          !   -      -
-      REAL(wp) ::   zWui, zWvi         !   -      -
-      REAL(wp) ::   zWus, zWvs         !   -      -
+      INTEGER  ::   ji, jj, jk           ! dummy loop indices
+      INTEGER  ::   iku, ikv             ! local integers
+      REAL(wp) ::   zzwi, ze3ua, zDt_2   ! local scalars
+      REAL(wp) ::   zzws, ze3va          !   -      -
+      REAL(wp) ::   z1_e3ua, z1_e3va     !   -      -
+      REAL(wp) ::   zWu , zWv            !   -      -
+      REAL(wp) ::   zWui, zWvi           !   -      -
+      REAL(wp) ::   zWus, zWvs           !   -      -
       REAL(wp), DIMENSION(jpi,jpj,jpk)        ::  zwi, zwd, zws   ! 3D workspace 
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv   !  -      -
@@ -93,4 +93,7 @@
          ENDIF
       ENDIF
+      !
+      zDt_2 = rDt * 0.5_wp
+      !
       !                             !* explicit top/bottom drag case
       IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, Kmm, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add top/bottom friction trend to (puu(Kaa),pvv(Kaa))
@@ -138,6 +141,6 @@
             ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    &
                &             + r_vvl   * e3v(ji,jj,ikv,Kaa)
-            puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua
-            pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
+            puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua
+            pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 *( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
          END_2D
          IF( ln_isfcav.OR.ln_drgice_imp ) THEN    ! Ocean cavities (ISF)
@@ -149,6 +152,6 @@
                ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    &
                   &             + r_vvl   * e3v(ji,jj,ikv,Kaa)
-               puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua
-               pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
+               puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua
+               pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 *( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
             END_2D
          END IF
@@ -158,5 +161,4 @@
       !
       !                    !* Matrix construction
-      zdt = rDt * 0.5
       IF( ln_zad_Aimp ) THEN   !!
          SELECT CASE( nldf_dyn )
@@ -165,13 +167,13 @@
                ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &
                   &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
-               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
+               zzwi = - zDt_2 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
                   &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
-               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
+               zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
                   &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
                zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua
                zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
-               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 
-               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
-               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
+               zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp ) 
+               zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWus, 0._wp )
+               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
             END_3D
          CASE DEFAULT               ! iso-level lateral mixing
@@ -179,13 +181,13 @@
                ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &    ! after scale factor at U-point
                   &             + r_vvl   * e3u(ji,jj,jk,Kaa)
-               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) )   &
+               zzwi = - zDt_2 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) )   &
                   &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
-               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) )   &
+               zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) )   &
                   &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
                zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua
                zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
-               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )
-               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
-               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
+               zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp )
+               zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWus, 0._wp )
+               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
             END_3D
          END SELECT
@@ -194,9 +196,9 @@
             ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm)    &
                &             + r_vvl   * e3u(ji,jj,1,Kaa)
-            zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) )   &
+            zzws = - zDt_2 * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) )   &
                &         / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
             zWus = ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) ) / ze3ua
-            zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp )
-            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) )
+            zws(ji,jj,1 ) = zzws - zDt_2 * MAX( zWus, 0._wp )
+            zwd(ji,jj,1 ) = 1._wp - zzws - zDt_2 * ( MIN( zWus, 0._wp ) )
          END_2D
       ELSE
@@ -206,7 +208,7 @@
                ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &
                   &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
-               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
+               zzwi = - zDt_2 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
                   &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
-               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
+               zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
                   &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
                zwi(ji,jj,jk) = zzwi
@@ -218,7 +220,7 @@
                ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &
                   &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
-               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) )    &
+               zzwi = - zDt_2 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) )    &
                   &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
-               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) )    &
+               zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) )    &
                   &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
                zwi(ji,jj,jk) = zzwi
@@ -245,5 +247,5 @@
             ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    &
                &             + r_vvl   * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
-            zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
+            zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
          END_2D
          IF ( ln_isfcav.OR.ln_drgice_imp ) THEN   ! top friction (always implicit)
@@ -253,5 +255,5 @@
                ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    &
                   &             + r_vvl   * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
-               zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
+               zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
             END_2D
          END IF
@@ -280,5 +282,5 @@
          ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm)    &
             &             + r_vvl   * e3u(ji,jj,1,Kaa) 
-         puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
+         puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + zDt_2 * ( utau_b(ji,jj) + utau(ji,jj) )   &
             &                                      / ( ze3ua * rho0 ) * umask(ji,jj,1) 
       END_2D
@@ -297,5 +299,4 @@
       !
       !                       !* Matrix construction
-      zdt = rDt * 0.5
       IF( ln_zad_Aimp ) THEN   !!
          SELECT CASE( nldf_dyn )
@@ -304,13 +305,13 @@
                ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
                   &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
-               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
+               zzwi = - zDt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
                   &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
-               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
+               zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
                   &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
                zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va
                zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
-               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp )
-               zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp )
-               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
+               zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWvi, 0._wp )
+               zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWvs, 0._wp )
+               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
             END_3D
          CASE DEFAULT               ! iso-level lateral mixing
@@ -318,13 +319,13 @@
                ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
                   &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
-               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) )    &
+               zzwi = - zDt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) )    &
                   &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
-               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) )    &
+               zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) )    &
                   &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
                zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va
                zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
-               zwi(ji,jj,jk) = zzwi  + zdt * MIN( zWvi, 0._wp )
-               zws(ji,jj,jk) = zzws  - zdt * MAX( zWvs, 0._wp )
-               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
+               zwi(ji,jj,jk) = zzwi  + zDt_2 * MIN( zWvi, 0._wp )
+               zws(ji,jj,jk) = zzws  - zDt_2 * MAX( zWvs, 0._wp )
+               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
             END_3D
          END SELECT
@@ -333,9 +334,9 @@
             ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm)    &
                &             + r_vvl   * e3v(ji,jj,1,Kaa)
-            zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) )    &
+            zzws = - zDt_2 * ( avm(ji,jj+1,2) + avm(ji,jj,2) )    &
                &         / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
             zWvs = ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) ) / ze3va
-            zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp )
-            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) )
+            zws(ji,jj,1 ) = zzws - zDt_2 * MAX( zWvs, 0._wp )
+            zwd(ji,jj,1 ) = 1._wp - zzws - zDt_2 * ( MIN( zWvs, 0._wp ) )
          END_2D
       ELSE
@@ -345,7 +346,7 @@
                ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
                   &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
-               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
+               zzwi = - zDt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
                   &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
-               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
+               zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
                   &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
                zwi(ji,jj,jk) = zzwi
@@ -357,7 +358,7 @@
                ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
                   &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
-               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) )    &
+               zzwi = - zDt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) )    &
                   &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
-               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) )    &
+               zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) )    &
                   &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
                zwi(ji,jj,jk) = zzwi
@@ -383,5 +384,5 @@
             ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    &
                &             + r_vvl   * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
-            zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
+            zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - zDt_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
          END_2D
          IF ( ln_isfcav.OR.ln_drgice_imp ) THEN
@@ -390,5 +391,5 @@
                ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    &
                   &             + r_vvl   * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
-               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va
+               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - zDt_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va
             END_2D
          ENDIF
@@ -417,5 +418,5 @@
          ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm)    &
             &             + r_vvl   * e3v(ji,jj,1,Kaa) 
-         pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
+         pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + zDt_2*( vtau_b(ji,jj) + vtau(ji,jj) )   &
             &                                      / ( ze3va * rho0 ) * vmask(ji,jj,1) 
       END_2D
@@ -432,6 +433,6 @@
       !
       IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
-         ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / rDt - ztrdu(:,:,:)
-         ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / rDt - ztrdv(:,:,:)
+         ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) )*r1_Dt - ztrdu(:,:,:)
+         ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) )*r1_Dt - ztrdv(:,:,:)
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt, Kmm )
          DEALLOCATE( ztrdu, ztrdv ) 
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/sshwzv.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/sshwzv.F90	(revision 14428)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/sshwzv.F90	(revision 14547)
@@ -17,5 +17,7 @@
    !!   ssh_nxt       : after ssh
    !!   ssh_atf       : time filter the ssh arrays
-   !!   wzv           : compute now vertical velocity
+   !!   wzv           : generic interface of vertical velocity calculation
+   !!   wzv_MLF       : MLF: compute NOW vertical velocity
+   !!   wzv_RK3       : RK3: Compute a vertical velocity
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers variables
@@ -44,4 +46,8 @@
    IMPLICIT NONE
    PRIVATE
+   !                  !! * Interface
+   INTERFACE wzv
+      MODULE PROCEDURE wzv_MLF, wzv_RK3
+   END INTERFACE
 
    PUBLIC   ssh_nxt        ! called by step.F90
@@ -77,8 +83,9 @@
       INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! time level index
       REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! sea-surface height
-      ! 
+      !
       INTEGER  ::   jk      ! dummy loop index
       REAL(wp) ::   zcoef   ! local scalar
       REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace
       !!----------------------------------------------------------------------
       !
@@ -92,10 +99,9 @@
       !
       zcoef = 0.5_wp * r1_rho0
-
       !                                           !------------------------------!
       !                                           !   After Sea Surface Height   !
       !                                           !------------------------------!
       IF(ln_wd_il) THEN
-         CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv )
+         CALL wad_lmt( pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv )
       ENDIF
 
@@ -133,8 +139,8 @@
    END SUBROUTINE ssh_nxt
 
-   
-   SUBROUTINE wzv( kt, Kbb, Kmm, Kaa, pww )
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE wzv  ***
+
+   SUBROUTINE wzv_MLF( kt, Kbb, Kmm, Kaa, pww )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE wzv_MLF  ***
       !!                   
       !! ** Purpose :   compute the now vertical velocity
@@ -157,10 +163,10 @@
       !!----------------------------------------------------------------------
       !
-      IF( ln_timing )   CALL timing_start('wzv')
+      IF( ln_timing )   CALL timing_start('wzv_MLF')
       !
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
-         IF(lwp) WRITE(numout,*) '~~~~~ '
+         IF(lwp) WRITE(numout,*) 'wzv_MLF : now vertical velocity '
+         IF(lwp) WRITE(numout,*) '~~~~~~~'
          !
          pww(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
@@ -193,5 +199,5 @@
          END DO
          !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
-         DEALLOCATE( zhdiv ) 
+         DEALLOCATE( zhdiv )
          !                                            !=================================!
       ELSEIF( ln_linssh )   THEN                      !==  linear free surface cases  ==!
@@ -204,7 +210,15 @@
          !                                            !==========================================!
          DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence
-            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)                 &
-               &                            + r1_Dt * (  e3t(:,:,jk,Kaa)        &
-               &                                       - e3t(:,:,jk,Kbb)  )   ) * tmask(:,:,jk)
+            !                                               ! NB: [e3t[a] -e3t[b] ]=e3t_0*[r3t[a]-r3t[b]]
+!!gm old
+!            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)   &
+!               &                            + r1_Dt * (  e3t(:,:,jk,Kaa)      &
+!               &                                       - e3t(:,:,jk,Kbb)  )   ) * tmask(:,:,jk)
+!!gm slightly faster :
+            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)                            &
+               &                            + r1_Dt * e3t_0(:,:,jk) * ( r3t(:,:,Kaa) - r3t(:,:,Kbb) )  ) * tmask(:,:,jk)
+!!gm end
+
+
          END DO
       ENDIF
@@ -256,7 +270,137 @@
 #endif
       !
-      IF( ln_timing )   CALL timing_stop('wzv')
-      !
-   END SUBROUTINE wzv
+      IF( ln_timing )   CALL timing_stop('wzv_MLF')
+      !
+   END SUBROUTINE wzv_MLF
+
+
+   SUBROUTINE wzv_RK3( kt, Kbb, Kmm, Kaa, puu, pvv, pww )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE wzv_RK3  ***
+      !!                   
+      !! ** Purpose :   compute the now vertical velocity
+      !!
+      !! ** Method  : - Using the incompressibility hypothesis, the vertical 
+      !!      velocity is computed by integrating the horizontal divergence  
+      !!      from the bottom to the surface minus the scale factor evolution.
+      !!        The boundary conditions are w=0 at the bottom (no flux) and.
+      !!
+      !! ** action  :   pww      : now vertical velocity
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in   ) ::   kt             ! time step
+      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   puu, pvv       ! horizontal velocity at Kmm
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::             pww  ! vertical velocity at Kmm
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zhdiv
+      REAL(wp)             , DIMENSION(jpi,jpj,jpk) ::   ze3div
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('wzv_RK3')
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'wzv_RK3 : now vertical velocity '
+         IF(lwp) WRITE(numout,*) '~~~~~ '
+         !
+         pww(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
+      ENDIF
+      !
+      CALL div_hor( kt, Kbb, Kmm, puu, pvv, ze3div )
+      !                                           !------------------------------!
+      !                                           !     Now Vertical Velocity    !
+      !                                           !------------------------------!
+      !
+      !                                               !===============================!
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      !==  z_tilde and layer cases  ==!
+         !                                            !===============================!
+         ALLOCATE( zhdiv(jpi,jpj,jpk) ) 
+         !
+         DO jk = 1, jpkm1
+            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
+            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
+            DO_2D( 0, 0, 0, 0 )
+               zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
+            END_2D
+         END DO
+         CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
+         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
+         !                             ! Same question holds for hdiv. Perhaps just for security
+         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
+            ! computation of w
+            pww(:,:,jk) = pww(:,:,jk+1) - (   ze3div(:,:,jk) + zhdiv(:,:,jk)   &
+               &                            + r1_Dt * (  e3t(:,:,jk,Kaa)       &
+               &                                       - e3t(:,:,jk,Kbb) )   ) * tmask(:,:,jk)
+         END DO
+         !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
+         DEALLOCATE( zhdiv ) 
+         !                                            !=================================!
+      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) - ze3div(:,:,jk)
+         END DO
+         !                                            !==========================================!
+      ELSE                                            !==  Quasi-Eulerian vertical coordinate  ==!   ('key_qco')
+         !                                            !==========================================!
+         DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence
+            !                                               ! NB: [e3t[a] -e3t[b] ]=e3t_0*[r3t[a]-r3t[b]]
+            pww(:,:,jk) = pww(:,:,jk+1) - (  ze3div(:,:,jk)                            &
+               &                            + r1_Dt * e3t_0(:,:,jk) * ( r3t(:,:,Kaa) - r3t(:,:,Kbb) )  ) * tmask(:,:,jk)
+         END DO
+      ENDIF
+
+      IF( ln_bdy ) THEN
+         DO jk = 1, jpkm1
+            pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:)
+         END DO
+      ENDIF
+      !
+#if defined key_agrif
+      IF( .NOT. AGRIF_Root() ) THEN
+         !
+         ! Mask vertical velocity at first/last columns/row 
+         ! inside computational domain (cosmetic) 
+         DO jk = 1, jpkm1
+            IF( lk_west ) THEN                             ! --- West --- !
+               DO ji = mi0(2+nn_hls), mi1(2+nn_hls)
+                  DO jj = 1, jpj
+                     pww(ji,jj,jk) = 0._wp 
+                  END DO
+               END DO
+            ENDIF
+            IF( lk_east ) THEN                             ! --- East --- !
+               DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls)
+                  DO jj = 1, jpj
+                     pww(ji,jj,jk) = 0._wp
+                  END DO
+               END DO
+            ENDIF
+            IF( lk_south ) THEN                            ! --- South --- !
+               DO jj = mj0(2+nn_hls), mj1(2+nn_hls)
+                  DO ji = 1, jpi
+                     pww(ji,jj,jk) = 0._wp
+                  END DO
+               END DO
+            ENDIF
+            IF( lk_north ) THEN                            ! --- North --- !
+               DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls)
+                  DO ji = 1, jpi
+                     pww(ji,jj,jk) = 0._wp
+                  END DO
+               END DO
+            ENDIF
+            !
+         END DO
+         !
+      ENDIF 
+#endif
+      !
+      IF( ln_timing )   CALL timing_stop('wzv_RK3')
+      !
+   END SUBROUTINE wzv_RK3
 
 
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/eosbn2.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/eosbn2.F90	(revision 14428)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/eosbn2.F90	(revision 14547)
@@ -56,5 +56,5 @@
    !                  !! * Interface
    INTERFACE eos
-      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d, eos_insitu_pot_2d
+      MODULE PROCEDURE eos_insitu_New, eos_insitu, eos_insitu_pot, eos_insitu_2d, eos_insitu_pot_2d
    END INTERFACE
    !
@@ -188,4 +188,111 @@
    !!----------------------------------------------------------------------
 CONTAINS
+
+   SUBROUTINE eos_insitu_New( pts, Knn, prd )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE eos_insitu  ***
+      !!
+      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from
+      !!       potential temperature and salinity using an equation of state
+      !!       selected in the nameos namelist
+      !!
+      !! ** Method  :   prd(t,s,z) = ( rho(t,s,z) - rho0 ) / rho0
+      !!         with   prd    in situ density anomaly      no units
+      !!                t      TEOS10: CT or EOS80: PT      Celsius
+      !!                s      TEOS10: SA or EOS80: SP      TEOS10: g/kg or EOS80: psu
+      !!                z      depth                        meters
+      !!                rho    in situ density              kg/m^3
+      !!                rho0   reference density            kg/m^3
+      !!
+      !!     ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z).
+      !!         Check value: rho = 1028.21993233072 kg/m^3 for z=3000 dbar, ct=3 Celsius, sa=35.5 g/kg
+      !!
+      !!     ln_eos80 : polynomial EOS-80 equation of state is used for rho(t,s,z).
+      !!         Check value: rho = 1028.35011066567 kg/m^3 for z=3000 dbar, pt=3 Celsius, sp=35.5 psu
+      !!
+      !!     ln_seos : simplified equation of state
+      !!              prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rho0
+      !!              linear case function of T only: rn_alpha<>0, other coefficients = 0
+      !!              linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0
+      !!              Vallis like equation: use default values of coefficients
+      !!
+      !! ** Action  :   compute prd , the in situ density (no units)
+      !!
+      !! References :   Roquet et al, Ocean Modelling, in preparation (2014)
+      !!                Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006
+      !!                TEOS-10 Manual, 2010
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:,:,:), INTENT(in   ) ::   pts   ! T-S
+      INTEGER                     , INTENT(in   ) ::   Knn   ! time-level
+      REAL(wp), DIMENSION(:,:,:  ), INTENT(  out) ::   prd   ! in situ density
+      !
+      INTEGER  ::   ji, jj, jk                ! dummy loop indices
+      REAL(wp) ::   zt , zh , zs , ztm        ! local scalars
+      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('eos-insitu')
+      !
+      SELECT CASE( neos )
+      !
+      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
+         !
+         DO_3D(nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+            !
+            zh  = gdept(ji,jj,jk,Knn) * r1_Z0                                 ! depth
+            zt  = pts (ji,jj,jk,jp_tem,Knn) * r1_T0                           ! temperature
+            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal,Knn) + rdeltaS ) * r1_S0 )   ! square root salinity
+            ztm = tmask(ji,jj,jk)                                             ! tmask
+            !
+            zn3 = EOS013*zt   &
+               &   + EOS103*zs+EOS003
+               !
+            zn2 = (EOS022*zt   &
+               &   + EOS112*zs+EOS012)*zt   &
+               &   + (EOS202*zs+EOS102)*zs+EOS002
+               !
+            zn1 = (((EOS041*zt   &
+               &   + EOS131*zs+EOS031)*zt   &
+               &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
+               &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
+               &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
+               !
+            zn0 = (((((EOS060*zt   &
+               &   + EOS150*zs+EOS050)*zt   &
+               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
+               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
+               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
+               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
+               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+               !
+            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+            !
+            prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm  ! density anomaly (masked)
+            !
+         END_3D
+         !
+      CASE( np_seos )                !==  simplified EOS  ==!
+         !
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+            zt  = pts  (ji,jj,jk,jp_tem,Knn) - 10._wp
+            zs  = pts  (ji,jj,jk,jp_sal,Knn) - 35._wp
+            zh  = gdept(ji,jj,jk,Knn)
+            ztm = tmask(ji,jj,jk)
+            !
+            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   &
+               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
+               &  - rn_nu * zt * zs
+               !
+            prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked)
+         END_3D
+         !
+      END SELECT
+      !
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu  : ', kdim=jpk )
+      !
+      IF( ln_timing )   CALL timing_stop('eos-insitu')
+      !
+   END SUBROUTINE eos_insitu_New
+
 
    SUBROUTINE eos_insitu( pts, prd, pdep )
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/tranpc.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/tranpc.F90	(revision 14428)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/tranpc.F90	(revision 14547)
@@ -71,5 +71,5 @@
       LOGICAL  ::   l_bottom_reached, l_column_treated
       REAL(wp) ::   zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z
-      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_rDt
+      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw
       REAL(wp), PARAMETER ::   zn2_zero = 1.e-14_wp             ! acceptance criteria for neutrality (N2==0)
       REAL(wp), DIMENSION(    jpk     )   ::   zvn2             ! vertical profile of N2 at 1 given point...
@@ -311,7 +311,6 @@
          !
          IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic
-            z1_rDt = 1._wp / (2._wp * rn_Dt)
-            ztrdt(:,:,:) = ( pts(:,:,:,jp_tem,Kaa) - ztrdt(:,:,:) ) * z1_rDt
-            ztrds(:,:,:) = ( pts(:,:,:,jp_sal,Kaa) - ztrds(:,:,:) ) * z1_rDt
+            ztrdt(:,:,:) = ( pts(:,:,:,jp_tem,Kaa) - ztrdt(:,:,:) ) * r1_Dt
+            ztrds(:,:,:) = ( pts(:,:,:,jp_sal,Kaa) - ztrds(:,:,:) ) * r1_Dt
             CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_npc, ztrdt )
             CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_npc, ztrds )
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/nemogcm.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/nemogcm.F90	(revision 14428)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/nemogcm.F90	(revision 14547)
@@ -67,5 +67,9 @@
 #endif
 #if defined key_qco   ||   defined key_linssh
+# if defined key_RK3
+   USE stprk3
+# else
    USE stpmlf         ! NEMO time-stepping               (stp_MLF   routine)
+# endif
 #else
    USE step           ! NEMO time-stepping                 (stp     routine)
@@ -163,5 +167,9 @@
          !
 #  if defined key_qco   ||   defined key_linssh
+#   if defined key_RK3
+         CALL stp_RK3
+#   else
          CALL stp_MLF
+#   endif
 #  else
          CALL stp
@@ -184,5 +192,9 @@
             !
 #  if defined key_qco   ||   defined key_linssh
+#   if defined key_RK3
+            CALL stp_RK3( istp )
+#   else
             CALL stp_MLF( istp )
+#   endif
 #  else
             CALL stp    ( istp )
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stp2d.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stp2d.F90	(revision 14547)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stp2d.F90	(revision 14547)
@@ -0,0 +1,267 @@
+MODULE stp2d
+   !!======================================================================
+   !!                       ***  MODULE stp2d  ***
+   !! Time-stepping   : manager of the ocean, tracer and ice time stepping
+   !!                   using a 3rd order Rung Kuta  with fixed or quasi-eulerian coordinate
+   !!======================================================================
+   !! History :  4.5  !  2021-01  (S. Techene, G. Madec, N. Ducousso, F. Lemarie)  Original code
+   !!   NEMO     
+   !!----------------------------------------------------------------------
+#if defined key_qco   ||   defined key_linssh
+   !!----------------------------------------------------------------------
+   !!   'key_qco'                        Quasi-Eulerian vertical coordinate
+   !!                          OR
+   !!   'key_linssh                       Fixed in time vertical coordinate
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!    stp_2D   : RK3 case
+   !!----------------------------------------------------------------------
+   USE step_oce       ! time stepping used modules
+   USE domqco         ! quasi-eulerian coordinate      (dom_qco_r3c routine)
+   USE dynadv_cen2    ! centred flux form advection    (dyn_adv_cen2 routine)
+   USE dynadv_ubs     ! UBS flux form advection        (dyn_adv_ubs  routine)
+   USE dynkeg         ! kinetic energy gradient        (dyn_keg      routine)
+   USE dynspg_ts      ! 2D mode integration
+   USE sbc_ice , ONLY : snwice_mass, snwice_mass_b
+   USE sbcapr         ! surface boundary condition: atmospheric pressure
+   USE sbcwave,  ONLY : bhd_wave
+
+   PRIVATE
+
+   PUBLIC   stp_2D   ! called by nemogcm.F90
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: step.F90 12377 2020-02-12 14:39:06Z acc $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE stp_2D( kt, Kbb, Kmm, Kaa, Krhs )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE stp_2D  ***
+      !!
+      !! ** Purpose : - Compute sea-surface height and barotropic velocity at Kaa
+      !!                in single 1st RK3.
+      !!
+      !! ** Method  : -1- Compute the 3D to 2D forcing
+      !!                 * Momentum (Ue,Ve)_rhs :
+      !!                      3D to 2D dynamics, i.e. the vertical sum of :
+      !!                        - Hor. adv. : KEG   + RVO in vector form
+      !!                                    : ADV_h + MET in flux   form
+      !!                        - LDF Lateral mixing 
+      !!                        - HPG Hor. pressure gradient
+      !!                      External forcings
+      !!                        - baroclinic drag
+      !!                        - wind 
+      !!                        - atmospheric pressure
+      !!                        - snow+ice load
+      !!                        - surface wave load
+      !!                 * ssh (sshe_rhs) :
+      !!                      Net column average freshwater flux
+      !!
+      !!              -2- Solve the external mode Eqs. using sub-time step
+      !!                  by a call to dyn_spg_ts (will be renamed dyn_2D or stp_2D)
+      !!
+      !! ** action  :   ssh            : N+1 sea surface height (Kaa=N+1)
+      !!                (uu_b,vv_b)    : N+1 barotropic velocity 
+      !!                (un_adv,vn_adv): barotropic transport from N to N+1 
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt, Kbb, Kmm, Kaa, Krhs    ! ocean time-step and time-level indices
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zg_2, zintp, zgrho0r, zld, zztmp     ! local scalars
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zpice   ! 2D workspace
+      !! ---------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('stp_2D')
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'stp_2D : barotropic field in single first '
+         IF(lwp) WRITE(numout,*) '~~~~~~'
+      ENDIF
+      !
+      IF( ln_linssh ) THEN    !==  Compute ww(:,:,1)  ==!   (needed for momentum advection)
+!!gm  only in Flux Form, Vector Form  dzU_z=0 assumed to be zero 
+!!gm  ww(k=1) = div_h(uu_b) ==> modif dans dynadv                        <<<=== TO BE DONE
+      ENDIF
+
+      ALLOCATE( sshe_rhs(jpi,jpj) , Ue_rhs(jpi,jpj) , Ve_rhs(jpi,jpj) , CdU_u(jpi,jpj) , CdU_v(jpi,jpj) )
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      !   RHS of barotropic momentum  Eq.
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+      !                       !======================================!
+      !                       !==  Dynamics 2D RHS from 3D trends  ==!   (HADV + LDF + HPG) (No Coriolis trend)
+      !                       !======================================!
+
+      uu(:,:,:,Krhs) = 0._wp        ! set dynamics trends to zero
+      vv(:,:,:,Krhs) = 0._wp
+
+      SELECT CASE( n_dynadv )       !*  compute horizontal advection  only  *!
+      CASE( np_VEC_c2  )                                    !- vector form  ( HADV = KEG + VOR )
+         CALL dyn_keg( kt, nn_dynkeg, Kbb, uu, vv, Krhs )                  ! grad_h(KE)
+         CALL dyn_vor( kt,            Kbb, uu, vv, Krhs, np_RVO )          ! relative vorticity
+      CASE( np_FLX_c2  )                                    !-  flux form   ( HADV = ADV_h + MET )
+         CALL dyn_adv_cen2( kt      , Kbb, uu, vv, Krhs, no_zad = 1 )      ! 2nd order centered scheme
+         CALL dyn_vor     ( kt      , Kbb, uu, vv, Krhs, np_MET )          ! metric term
+      CASE( np_FLX_ubs )  
+         CALL dyn_adv_ubs ( kt      , Kbb, Kbb, uu, vv, Krhs, no_zad = 1)  ! 3rd order UBS scheme (UP3)
+         CALL dyn_vor     ( kt      , Kbb, uu, vv, Krhs, np_MET )          ! metric term
+      END SELECT
+      !
+      !                             !*  lateral viscosity  *!
+      CALL dyn_ldf( kt,   Kbb, Kbb, uu, vv, Krhs )
+      !
+      !                             !*  hydrostatic pressure gradient  *!  
+      CALL eos    ( ts , Kbb, rhd )                          ! in situ density anomaly at Kbb
+      CALL dyn_hpg( kt  , Kbb     , uu, vv, Krhs )           ! horizontal gradient of Hydrostatic pressure
+      !
+      !                             !*  vertical averaging  *!
+      Ue_rhs(:,:) = SUM( e3u_0(:,:,:) * uu(:,:,:,Krhs) * umask(:,:,:), DIM=3 ) * r1_hu_0(:,:)
+      Ve_rhs(:,:) = SUM( e3v_0(:,:,:) * vv(:,:,:,Krhs) * vmask(:,:,:), DIM=3 ) * r1_hv_0(:,:)
+      !
+      !
+      !                       !===========================!
+      !                       !==  external 2D forcing  ==!
+      !                       !===========================!
+      !
+      ! 			    !* baroclinic drag forcing *!   (also provide the barotropic drag coeff.)
+      !
+      CALL dyn_drg_init( Kbb, Kbb, uu, vv, uu_b, vv_b, Ue_rhs, Ve_rhs, CdU_u, CdU_v )
+      !
+      !                             !* wind forcing *!
+!!st ATTENTION stoke drift !!
+      IF( ln_bt_fw ) THEN
+         DO_2D( 0, 0, 0, 0 )
+            Ue_rhs(ji,jj) =  Ue_rhs(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kbb)
+            Ve_rhs(ji,jj) =  Ve_rhs(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kbb)
+         END_2D
+      ELSE
+         zztmp = r1_rho0 * r1_2
+         DO_2D( 0, 0, 0, 0 )
+            Ue_rhs(ji,jj) =  Ue_rhs(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kbb)
+            Ve_rhs(ji,jj) =  Ve_rhs(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kbb)
+         END_2D
+      ENDIF
+      !
+      !                             !* atmospheric pressure forcing *!
+      IF( ln_apr_dyn ) THEN
+         IF( ln_bt_fw ) THEN                          ! FORWARD integration: use kt+1/2 pressure (NOW+1/2)
+            DO_2D( 0, 0, 0, 0 )
+               Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
+               Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
+            END_2D
+         ELSE                                         ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW)
+            zztmp = grav * r1_2
+            DO_2D( 0, 0, 0, 0 )
+               Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + zztmp * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)  &
+                    &                                   + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
+               Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + zztmp * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)  &
+                    &                                   + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj)
+            END_2D
+         ENDIF
+      ENDIF
+      !
+      !                             !* snow+ice load *!   (embedded sea ice)
+      IF( ln_ice_embd ) THEN
+         ALLOCATE( zpice(jpi,jpj) )
+         zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc )
+         zgrho0r     = - grav * r1_rho0
+         zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrho0r
+         DO_2D( 0, 0, 0, 0 )
+            Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
+            Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
+         END_2D
+         DEALLOCATE( zpice )
+      ENDIF
+      !
+      !                             !* surface wave load *!   (Bernoulli head)
+      !
+      IF( ln_wave .AND. ln_bern_srfc ) THEN
+         DO_2D( 0, 0, 0, 0 )
+            Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) * r1_e1u(ji,jj)   !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3]
+            Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) * r1_e1u(ji,jj)
+         END_2D
+      ENDIF
+      !
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      !   RHS of see surface height  Eq.
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      ! 		!==  Net water flux forcing  ==!  (applied to a water column)
+      !
+      IF (ln_bt_fw) THEN                          ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2)
+         sshe_rhs(:,:) = r1_rho0 * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) )
+      ELSE                                        ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW)
+         zztmp = r1_rho0 * r1_2
+         sshe_rhs(:,:) = zztmp * (   emp(:,:)        + emp_b(:,:)          &
+            &                      - rnf(:,:)        - rnf_b(:,:)          &
+            &                      + fwfisf_cav(:,:) + fwfisf_cav_b(:,:)   &
+            &                      + fwfisf_par(:,:) + fwfisf_par_b(:,:)   )
+      ENDIF
+      !
+      ! 		!==  Stokes drift divergence  ==!   (if exist)
+      !
+      IF( ln_sdw )    sshe_rhs(:,:) = sshe_rhs(:,:) + div_sd(:,:)
+      !
+      !
+      ! 		!==  ice sheet coupling  ==!
+      !
+      IF( ln_isf .AND. ln_isfcpl ) THEN
+         IF( ln_rstart .AND. kt == nit000 )   sshe_rhs(:,:) = sshe_rhs(:,:) + risfcpl_ssh(:,:)
+         IF( ln_isfcpl_cons               )   sshe_rhs(:,:) = sshe_rhs(:,:) + risfcpl_cons_ssh(:,:)
+      ENDIF
+      !
+#if defined key_asminc
+      !                 !==  Add the IAU weighted SSH increment  ==!
+      !
+      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau )   zsshe_rhs(:,:) = zsshe_rhs(:,:) - ssh_iau(:,:)
+#endif
+      !
+#if defined key_agrif
+      !                 !==  AGRIF : fill boundary data arrays (on both )
+         IF( .NOT.Agrif_Root() )   CALL agrif_dta_ts( kt )
+#endif
+      !
+      !
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      !             Compute ssh and (uu_b,vv_b)  at N+1  (Kaa)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+      !    using a split-explicit time integration in forward mode 
+      !    ( ABM3-AM4 time-integration Shchepetkin et al. OM2005) with temporal diffusion (Demange et al. JCP2019) )
+
+      CALL dyn_spg_ts( kt, Kbb, Kbb, Krhs, uu, vv, ssh, uu_b, vv_b, Kaa ) ! time-splitting
+                   
+
+      DEALLOCATE( sshe_rhs , Ue_rhs , Ve_rhs , CdU_u , CdU_v )
+
+!!gm  this is useless I guess : RK3,  done in each stages
+!
+!      IF( ln_dynspg_ts ) THEN      ! With split-explicit free surface, since now transports have been updated and ssh(:,:,Krhs)
+!                                   ! as well as vertical scale factors and vertical velocity need to be updated
+!                            CALL div_hor    ( kstp, Kbb, Kmm )                ! Horizontal divergence  (2nd call in time-split case)
+!         IF(.NOT.lk_linssh) CALL dom_qco_r3c( ssh(:,:,Kaa), r3t(:,:,Kaa), r3u(:,:,Kaa), r3v(:,:,Kaa), r3f(:,:) )   ! update ssh/h_0 ratio at t,u,v,f pts 
+!      ENDIF
+!!gm   
+      !
+      IF( ln_timing )   CALL timing_stop('stp_2D')
+      !
+   END SUBROUTINE stp_2D
+
+
+#else
+   !!----------------------------------------------------------------------
+   !!   default option             EMPTY MODULE           qco not activated
+   !!----------------------------------------------------------------------
+#endif
+   
+   !!======================================================================
+END MODULE stp2d
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stpmlf.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stpmlf.F90	(revision 14428)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stpmlf.F90	(revision 14547)
@@ -62,5 +62,4 @@
 #  include "do_loop_substitute.h90"
 #  include "domzgr_substitute.h90"
-#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stprk3.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stprk3.F90	(revision 14547)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stprk3.F90	(revision 14547)
@@ -0,0 +1,394 @@
+MODULE stprk3
+   !!======================================================================
+   !!                       ***  MODULE stpRK3  ***
+   !! Time-stepping   : manager of the ocean, tracer and ice time stepping
+   !!                   using a 3rd order Rung Kuta  with fixed or quasi-eulerian coordinate
+   !!======================================================================
+   !! History :  4.5  !  2021-01  (S. Techene, G. Madec, N. Ducousso, F. Lemarie)  Original code
+   !!   NEMO     
+   !!----------------------------------------------------------------------
+#if defined key_qco   ||   defined key_linssh
+   !!----------------------------------------------------------------------
+   !!   'key_qco'                        Quasi-Eulerian vertical coordinate
+   !!                          OR
+   !!   'key_linssh                       Fixed in time vertical coordinate
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   stp_RK3       : NEMO 3rd order Runge-Kutta time-stepping 
+   !!----------------------------------------------------------------------
+   USE step_oce       ! time stepping used modules
+   USE domqco         ! quasi-eulerian coordinate      (dom_qco_r3c routine)
+   USE rk3stg         ! RK3 stages
+   USE stp2d          ! external mode solver
+   USE dynspg_ts, ONLY: un_adv, vn_adv   ! updated Kmm barotropic transport
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   stp_RK3   ! called by nemogcm.F90
+
+   !                                          !**  time level indices  **!
+   INTEGER, PUBLIC ::   Nbb, Nnn, Naa, Nrhs   !: used by nemo_init
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: step.F90 12377 2020-02-12 14:39:06Z acc $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+#if defined key_agrif
+   RECURSIVE SUBROUTINE stp_RK3( )
+      INTEGER             ::   kstp   ! ocean time-step index
+#else
+   SUBROUTINE stp_RK3( kstp )
+      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
+#endif
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE stp_RK3  ***
+      !!
+      !! ** Purpose : - Time stepping of OCE  (momentum and active tracer Eqs.) (RK3)
+      !!              - Time stepping of SI3 (dynamic and thermodynamic Eqs.)   (FBS)
+      !!              - Time stepping of TRC  (passive tracer Eqs.)
+      !!
+      !! ** Method  : -1- Update forcings and data
+      !!              -2- Update ocean physics
+      !!              -3- Compute the after (Naa) ssh and velocity 
+      !!              -4- diagnostics and output at Now (Nnn)
+      !!              -4- Compute the after (Naa) T-S
+      !!              -5- Update now 
+      !!              -6- Update the horizontal velocity
+      !!              -7- Compute the diagnostics variables (rd,N2, hdiv,w)
+      !!              -8- Outputs and diagnostics
+      !!----------------------------------------------------------------------
+      INTEGER ::   ji, jj, jk, jtile   ! dummy loop indice
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zgdept
+      !! ---------------------------------------------------------------------
+#if defined key_agrif
+      IF( nstop > 0 )   RETURN   ! avoid to go further if an error was detected during previous time step (child grid)
+      kstp = nit000 + Agrif_Nb_Step()
+      Kbb_a = Nbb   ; Kmm_a = Nnn   ;   Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
+      IF( lk_agrif_debug ) THEN
+         IF( Agrif_Root() .AND. lwp)   WRITE(*,*) '---'
+         IF(lwp)   WRITE(*,*) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint()
+      ENDIF
+      IF( kstp == nit000 + 1 )   lk_agrif_fstep = .FALSE.
+# if defined key_xios
+      IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( cxios_context )
+# endif
+#endif
+      !
+      IF( ln_timing )   CALL timing_start('stp_RK3')
+      !
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! update I/O and calendar
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
+                             CALL iom_init( cxios_context, ld_closedef=.FALSE. )   ! for model grid (including possible AGRIF zoom)
+         IF( lk_diamlr   )   CALL dia_mlr_iom_init    ! with additional setup for multiple-linear-regression analysis
+                             CALL iom_init_closedef
+         IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid
+      ENDIF
+      IF( kstp == nitrst .AND. lwxios ) THEN
+                             CALL iom_swap(                     cw_ocerst_cxt )
+                             CALL iom_init_closedef(            cw_ocerst_cxt )
+                             CALL iom_setkt( kstp - nit000 + 1, cw_ocerst_cxt )
+#if defined key_top
+                             CALL iom_swap(                     cw_toprst_cxt )
+                             CALL iom_init_closedef(            cw_toprst_cxt )
+                             CALL iom_setkt( kstp - nit000 + 1, cw_toprst_cxt )
+#endif
+      ENDIF
+#if defined key_si3
+      IF( kstp + nn_fsbc - 1 == nitrst .AND. lwxios ) THEN
+                             CALL iom_swap(                     cw_icerst_cxt )
+                             CALL iom_init_closedef(            cw_icerst_cxt )
+                             CALL iom_setkt( kstp - nit000 + 1, cw_icerst_cxt )
+      ENDIF
+#endif
+      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
+                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell IOM we are at time step kstp
+      IF( ln_crs         )   CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell IOM we are at time step kstp
+
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Update external forcing (tides, open boundaries, ice shelf interaction and surface boundary condition (including sea-ice)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( ln_tide    )   CALL tide_update( kstp )                          ! update tide potential
+      IF( ln_apr_dyn )   CALL sbc_apr    ( kstp )                          ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib)
+      IF( ln_bdy     )   CALL bdy_dta    ( kstp, Nnn )                     ! update dynamic & tracer data at open boundaries
+      IF( ln_isf     )   CALL isf_stp    ( kstp, Nnn )                     ! update iceshelf geometry
+                         CALL sbc        ( kstp, Nbb, Nnn )                ! Sea Boundary Condition (including sea-ice)
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Update stochastic parameters and random T/S fluctuations
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      IF( ln_sto_eos )   CALL sto_par( kstp )                              ! Stochastic parameters
+      IF( ln_sto_eos )   CALL sto_pts( ts(:,:,:,:,Nnn)  )                  ! Random T/S fluctuations
+
+!!gm  ocean physic computed at stage 3 ?
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Ocean physics update
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !  THERMODYNAMICS
+!!gm only Before is needed for bn2 and rab  (except at stage 3 for n2)
+!!gm issue with   Nnn  used in rad(Nbb) 
+                         CALL eos_rab( ts(:,:,:,:,Nbb), rab_b, Nnn )       ! before local thermal/haline expension ratio at T-points
+!!st                         CALL eos_rab( ts(:,:,:,:,Nnn), rab_n, Nnn )       ! now    local thermal/haline expension ratio at T-points
+                         CALL bn2    ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency
+!!st                         CALL bn2    ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn  ) ! now    Brunt-Vaisala frequency
+!!gm
+                         rab_n = rab_b
+                         rn2   = rn2b
+!!gm sh2 computed at the end of the time-step
+!!gm       or  call zdf_phy at the end !
+      !  VERTICAL PHYSICS
+!!st                         CALL zdf_phy( kstp, Nbb, Nnn, Nrhs )   ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
+                         CALL zdf_phy( kstp, Nbb, Nbb, Nrhs )   ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
+!!gm
+      !  LATERAL  PHYSICS
+      !
+      IF( l_ldfslp ) THEN                             ! slope of lateral mixing
+!!gm gdep
+                         CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) )               ! before in situ density
+
+         IF( ln_zps .AND. .NOT. ln_isfcav)                                    &
+            &            CALL zps_hde    ( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
+            &                                          rhd, gru , grv    )       ! of t, s, rd at the last ocean level
+
+         IF( ln_zps .AND.       ln_isfcav)                                                &
+            &            CALL zps_hde_isf( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+            &                                          rhd, gru , grv , grui, grvi   )       ! of t, s, rd at the first ocean level
+         IF( ln_traldf_triad ) THEN
+                         CALL ldf_slp_triad( kstp, Nbb, Nnn )             ! before slope for triad operator
+         ELSE
+                         CALL ldf_slp     ( kstp, rhd, rn2b, Nbb, Nnn )   ! before slope for standard operator
+         ENDIF
+      ENDIF
+      !                                                                        ! eddy diffusivity coeff.
+      IF( l_ldftra_time .OR. l_ldfeiv_time )   CALL ldf_tra( kstp, Nbb, Nnn )  !       and/or eiv coeff.
+      IF( l_ldfdyn_time                    )   CALL ldf_dyn( kstp, Nbb )       ! eddy viscosity coeff.
+
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      !  RK3 : single first external mode computation
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+      CALL stp_2D( kstp, Nbb, Nbb, Naa, Nrhs )   ! out: ssh, (uu_b,vv_b) and (un_adv,vn_adv) at Naa
+
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      !  RK3 time integration
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+      ! Stage 1 :
+      CALL stp_RK3_stg( 1, kstp, Nbb, Nbb, Nrhs, Naa )
+      !
+      Nrhs = Nnn   ;   Nnn  = Naa   ;   Naa  = Nrhs    ! Swap: Nbb unchanged, Nnn <==> Naa
+      !
+      ! Stage 2 :
+      CALL stp_RK3_stg( 2, kstp, Nbb, Nnn, Nrhs, Naa )
+      !
+      Nrhs = Nnn   ;   Nnn  = Naa   ;   Naa  = Nrhs    ! Swap: Nbb unchanged, Nnn <==> Naa
+      !
+      ! Stage 3 :
+      CALL stp_RK3_stg( 3, kstp, Nbb, Nnn, Nrhs, Naa )
+      !
+      Nrhs = Nbb   ;   Nbb  = Naa   ;   Naa  = Nrhs    ! Swap: Nnn unchanged, Nbb <==> Naa
+    
+      ! Swap time levels 
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! diagnostics and outputs
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      
+!==>>>  at Nbb  no more Nnn 
+      
+      IF( ln_diacfl  )   CALL dia_cfl   ( kstp,      Nbb )      ! Courant number diagnostics
+                         CALL dia_hth   ( kstp,      Nbb )      ! Thermocline depth (20 degres isotherm depth)
+      IF( ln_diadct  )   CALL dia_dct   ( kstp,      Nbb )      ! Transports
+                         CALL dia_ar5   ( kstp,      Nbb )      ! ar5 diag
+                         CALL dia_ptr   ( kstp,      Nbb )      ! Poleward adv/ldf TRansports diagnostics
+                         CALL dia_wri   ( kstp,      Nbb )      ! ocean model: outputs
+      IF( ln_crs     )   CALL crs_fld   ( kstp,      Nbb )      ! ocean model: online field coarsening & output
+      IF( lk_diadetide ) CALL dia_detide( kstp )                ! Weights computation for daily detiding of model diagnostics
+      IF( lk_diamlr  )   CALL dia_mlr                           ! Update time used in multiple-linear-regression analysis
+
+!!====>>>> to be modified for RK3
+!      IF( ln_floats  )   CALL flo_stp   ( kstp, Nbb, Nnn )      ! drifting Floats
+!      IF( ln_diahsb  )   CALL dia_hsb   ( kstp, Nbb, Nnn )  ! - ML - global conservation diagnostics
+
+!!gm : This does not only concern the dynamics ==>>> add a new title
+!!gm2: why ouput restart before AGRIF update?
+!!
+!!jc: That would be better, but see comment above
+!!
+!!====>>>> to be modified for RK3
+      IF( lrst_oce   )   CALL rst_write    ( kstp, Nbb, Nnn )   ! write output ocean restart file
+      IF( ln_sto_eos )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters
+
+#if defined key_agrif
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! AGRIF recursive integration
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                         Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs      ! agrif_oce module copies of time level indices
+                         CALL Agrif_Integrate_ChildGrids( stp_RK3 )       ! allows to finish all the Child Grids before updating
+
+#endif
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Control
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                         CALL stp_ctl      ( kstp, Nnn )
+
+#if defined key_agrif
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! AGRIF update
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( Agrif_NbStepint() == 0 .AND. nstop == 0 )   &
+         &               CALL Agrif_update_all( )                  ! Update all components
+
+#endif
+      IF( ln_diaobs .AND. nstop == 0 )   &
+         &               CALL dia_obs( kstp, Nnn )  ! obs-minus-model (assimilation) diags (after dynamics update)
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! File manipulation at the end of the first time step
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( kstp == nit000 ) THEN                          ! 1st time step only
+                                        CALL iom_close( numror )   ! close input  ocean restart file
+         IF( lrxios )                   CALL iom_context_finalize( cr_ocerst_cxt )
+         IF(lwm)                        CALL FLUSH    ( numond )   ! flush output namelist oce
+         IF(lwm .AND. numoni /= -1 )    CALL FLUSH    ( numoni )   ! flush output namelist ice (if exist)
+      ENDIF
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Coupled mode
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( lk_oasis .AND. nstop == 0 )   CALL sbc_cpl_snd( kstp, Nbb, Nnn )     ! coupled mode : field exchanges
+      !
+#if defined key_xios
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Finalize contextes if end of simulation or error detected
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( kstp == nitend .OR. nstop > 0 ) THEN
+                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
+         IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !
+      ENDIF
+#endif
+      !
+      IF( l_1st_euler ) THEN         ! recover Leap-frog timestep
+         rDt   = 2._wp * rn_Dt
+         r1_Dt = 1._wp / rDt
+         l_1st_euler = .FALSE.
+      ENDIF
+      !
+      IF( ln_timing )   CALL timing_stop('stp_RK3')
+      !
+   END SUBROUTINE stp_RK3
+
+
+   SUBROUTINE mlf_baro_corr( Kmm, Kaa, puu, pvv )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE mlf_baro_corr  ***
+      !!
+      !! ** Purpose :   Finalize after horizontal velocity.
+      !!
+      !! ** Method  : * Ensure after velocities transport matches time splitting
+      !!             estimate (ln_dynspg_ts=T)
+      !!
+      !! ** Action :   puu(Kmm),pvv(Kmm)   updated now horizontal velocity (ln_bt_fw=F)
+      !!               puu(Kaa),pvv(Kaa)   after horizontal velocity
+      !!----------------------------------------------------------------------
+      USE dynspg_ts, ONLY : un_adv, vn_adv   ! updated Kmm barotropic transport 
+      !!
+      INTEGER                             , INTENT(in   ) ::   Kmm, Kaa   ! before and after time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv   ! velocities
+      !
+      INTEGER  ::   jk   ! dummy loop indices
+      REAL(wp), DIMENSION(jpi,jpj) ::   zue, zve
+      !!----------------------------------------------------------------------
+
+      ! Ensure below that barotropic velocities match time splitting estimate
+      ! Compute actual transport and replace it with ts estimate at "after" time step
+      zue(:,:) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
+      zve(:,:) = e3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
+      DO jk = 2, jpkm1
+         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
+         puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk)
+         pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk)
+      END DO
+      !
+!!st      IF( .NOT.ln_bt_fw ) THEN
+!!st         ! Remove advective velocity from "now velocities"
+!!st         ! prior to asselin filtering
+!!st         ! In the forward case, this is done below after asselin filtering
+!!st         ! so that asselin contribution is removed at the same time
+!!st         DO jk = 1, jpkm1
+!!st            puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk)
+!!st            pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk)
+!!st         END DO
+!!st      ENDIF
+      !
+   END SUBROUTINE mlf_baro_corr
+
+
+   SUBROUTINE finalize_lbc( kt, Kbb, Kaa, puu, pvv, pts )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE finalize_lbc  ***
+      !!
+      !! ** Purpose :   Apply the boundary condition on the after velocity
+      !!
+      !! ** Method  : * Apply lateral boundary conditions on after velocity
+      !!             at the local domain boundaries through lbc_lnk call,
+      !!             at the one-way open boundaries (ln_bdy=T),
+      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
+      !!
+      !! ** Action :   puu(Kaa),pvv(Kaa)   after horizontal velocity and tracers
+      !!----------------------------------------------------------------------
+#if defined key_agrif
+      USE agrif_oce_interp
+#endif
+      USE bdydyn         ! ocean open boundary conditions (define bdy_dyn)
+      !!
+      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) ::   Kbb, 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,jpts,jpt), INTENT(inout) ::   pts        ! active tracers
+      !!----------------------------------------------------------------------
+      !
+      ! Update after tracer and velocity on domain lateral boundaries
+      !
+# if defined key_agrif
+            CALL Agrif_tra                     !* AGRIF zoom boundaries
+            CALL Agrif_dyn( kt )
+# endif
+      !                                        ! local domain boundaries  (T-point, unchanged sign)
+      CALL lbc_lnk_multi( 'finalize_lbc', puu(:,:,:,       Kaa), 'U', -1., pvv(:,:,:       ,Kaa), 'V', -1.   &
+                       &                , pts(:,:,:,jp_tem,Kaa), 'T',  1., pts(:,:,:,jp_sal,Kaa), 'T',  1. )
+      !
+      !                                        !* BDY open boundaries
+      IF( ln_bdy )   THEN
+                               CALL bdy_tra( kt, Kbb, pts,      Kaa )
+         IF( ln_dynspg_exp )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa )
+         IF( ln_dynspg_ts  )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only=.true. )
+      ENDIF
+      !
+   END SUBROUTINE finalize_lbc
+
+#else
+   !!----------------------------------------------------------------------
+   !!   default option             EMPTY MODULE           qco not activated
+   !!----------------------------------------------------------------------
+#endif
+   
+   !!======================================================================
+END MODULE stprk3
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stprk3_stg.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stprk3_stg.F90	(revision 14547)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stprk3_stg.F90	(revision 14547)
@@ -0,0 +1,362 @@
+MODULE stprk3_stg
+   !!======================================================================
+   !!                       ***  MODULE stprk3_stg  ***
+   !! Time-stepping   : manager of the ocean, tracer and ice time stepping
+   !!                   using a 3rd order Runge-Kutta  with fixed or quasi-eulerian coordinate
+   !!======================================================================
+   !! History :  4.5  !  2021-01  (S. Techene, G. Madec, N. Ducousso, F. Lemarie)  Original code
+   !!   NEMO     
+   !!----------------------------------------------------------------------
+#if defined key_qco   ||   defined key_linssh
+   !!----------------------------------------------------------------------
+   !!   'key_qco'                        Quasi-Eulerian vertical coordinate
+   !!                          OR
+   !!   'key_linssh                       Fixed in time vertical coordinate
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!    stp_RK3_stg  : NEMO 3rd order Runge-Kutta stage with qco or linssh
+   !!----------------------------------------------------------------------
+   USE step_oce       ! time stepping used modules
+   USE domqco         ! quasi-eulerian coordinate      (dom_qco_r3c routine)
+   USE dynspg_ts, ONLY:   un_adv , vn_adv   ! advective transport from N to N+1
+   USE bdydyn         ! ocean open boundary conditions (define bdy_dyn)
+#if defined key_agrif
+   USE agrif_oce_interp
+#endif
+   
+   !
+   USE prtctl         ! print control
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   stp_RK3_stg   ! called by nemogcm.F90
+
+   REAL(wp) ::   r1_2 = 1._wp / 2._wp
+   REAL(wp) ::   r1_3 = 1._wp / 3._wp
+   REAL(wp) ::   r2_3 = 2._wp / 3._wp
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssha         ! sea-surface height  at N+1
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ua_b, va_b   ! barotropic velocity at N+1
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: step.F90 12377 2020-02-12 14:39:06Z acc $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE stp_RK3_stg( kstg, kstp, Kbb, Kmm, Krhs, Kaa )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE stp_RK3_stg  ***
+      !!
+      !! ** Purpose : - stage of RK3 time stepping of OCE and TOP 
+      !!
+      !! ** Method  :   input: computed in dynspg_ts
+      !!              ssh             shea surface height at N+1           (oce.F90)
+      !!              (uu_b,vv_b)     barotropic velocity at N, N+1        (oce.F90)
+      !!              (un_adv,vn_adv) barotropic transport from N to N+1   (dynspg_ts.F90)
+      !!              , 
+      !!              -1- set ssh(Naa) (Naa=N+1/3, N+1/2, or N)
+      !!              -2- set the advective velocity (zadU,zaV)
+      !!              -4- Compute the after (Naa) T-S
+      !!              -5- Update now 
+      !!              -6- Update the horizontal velocity
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kstg                        ! RK3 stage
+      INTEGER, INTENT(in) ::   kstp, Kbb, Kmm, Krhs, Kaa   ! ocean time-step and time-level indices
+      !
+      INTEGER  ::   ji, jj, jk, jtile                      ! dummy loop indices
+      REAL(wp) ::   ze3Ub, ze3Vb, ze3Tb, ze3Sb, z1_e3t     ! local scalars
+      REAL(wp) ::   ze3Ur, ze3Vr, ze3Tr, ze3Sr             !   -      -
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zaU, zaV       ! advective horizontal velocity
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zub, zvb       ! advective transport 
+      !! ---------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('stp_RK3_stg')
+      !
+      IF( kstp == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'stp_RK3_stg : Runge Kutta 3rd order at stage ', kstg
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      ENDIF
+      !
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      !  ssh, uu_b, vv_b, and  ssh/h0 at Kaa
+      !  3D advective velocity at Kmm
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !     
+      SELECT CASE( kstg )
+      !                    !---------------!
+      CASE ( 1 )           !==  Stage 1  ==!   Kbb = Kmm = N  ;  Kaa = N+1/3
+         !                 !---------------!
+         !
+         ALLOCATE( ssha(jpi,jpj) , ua_b(jpi,jpj) , va_b(jpi,jpj) )
+         !
+         rDt = r1_3 * rn_Dt            ! set time-step : rn_Dt/3
+         r1_Dt = 1._wp / rDt
+         !
+         ssha(:,:) = ssh (:,:,Kaa)     ! save ssh, uu_b, vv_b at N+1  (computed in dynspg_ts)
+         ua_b(:,:) = uu_b(:,:,Kaa)
+         va_b(:,:) = vv_b(:,:,Kaa)
+         !                             ! interpolated ssh and (uu_b,vv_b) at Kaa (N+1/3)
+         ssh (:,:,Kaa) = r2_3 * ssh (:,:,Kbb) + r1_3 * ssha(:,:)
+         uu_b(:,:,Kaa) = r2_3 * uu_b(:,:,Kbb) + r1_3 * ua_b(:,:)
+         vv_b(:,:,Kaa) = r2_3 * vv_b(:,:,Kbb) + r1_3 * va_b(:,:)
+         !
+         !                 !---------------!
+      CASE ( 2 )           !==  Stage 2  ==!   Kbb = N   ;   Kmm = N+1/3   ;   Kaa = N+1/2
+         !                 !---------------!
+         !
+         rDt = r1_2 * rn_Dt            ! set time-step : rn_Dt/2
+         r1_Dt = 1._wp / rDt
+         !
+         !                             ! set ssh and (uu_b,vv_b) at N+1/2  (Kaa)
+         ssh (:,:,Kaa) = r1_2 * ( ssh (:,:,Kbb) + ssha(:,:) )
+         uu_b(:,:,Kaa) = r1_2 * ( uu_b(:,:,Kbb) + ua_b(:,:) )
+         vv_b(:,:,Kaa) = r1_2 * ( vv_b(:,:,Kbb) + va_b(:,:) )
+         !
+         !                 !---------------!
+      CASE ( 3 )           !==  Stage 3  ==!   Kbb = N   ;   Kmm = N+1/2   ;   Kaa = N+1
+         !                 !---------------!
+         !
+         rDt = rn_Dt                   ! set time-step : rn_Dt
+         r1_Dt = 1._wp / rDt
+         !
+         ssh (:,:,Kaa) = ssha(:,:)     ! recover ssh and (uu_b,vv_b) at N + 1
+         uu_b(:,:,Kaa) = ua_b(:,:)                     
+         vv_b(:,:,Kaa) = va_b(:,:)
+         !
+         DEALLOCATE( ssha , ua_b , va_b )
+         !
+      END SELECT
+      !
+      !                     !==  ssh/h0 ratio at Kaa  ==! 
+      !
+      CALL dom_qco_r3c( ssh(:,:,Kaa), r3t(:,:,Kaa), r3u(:,:,Kaa), r3v(:,:,Kaa), r3f(:,:) )   ! "after" ssh/h_0 ratio at t,u,v-column
+      !
+      !
+      !                     !==  advective velocity at Kmm  ==!
+      !
+      !                                            !- horizontal components -!   (zaU,zaV) 
+      zub(:,:) = un_adv(:,:)*r1_hu(:,:,Kmm) - uu_b(:,:,Kmm)    ! barotropic velocity correction
+      zvb(:,:) = vn_adv(:,:)*r1_hv(:,:,Kmm) - vv_b(:,:,Kmm)
+      DO jk = 1, jpkm1                                         ! horizontal advective velocity
+         zaU(:,:,jk) = uu(:,:,jk,Kmm) + zub(:,:)*umask(:,:,jk)
+         zaV(:,:,jk) = vv(:,:,jk,Kmm) + zvb(:,:)*vmask(:,:,jk)
+      END DO
+      !                                            !- vertical components -!   ww
+                         CALL wzv  ( kstp, Kbb, Kmm, Kaa, zaU, zaV, ww )     ! ww cross-level velocity
+
+!!st      IF( ln_zad_Aimp )  CALL wAimp( kstp, ww, wi       )     ! Adaptive-implicit vertical advection partitioning
+      !                                                            !==>>>  implicite for stages 1 & 2 ????
+      !
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      !  RHS of ocean dynamics : ADV + VOR/COR + HPG (+ ASM )      <<<===  Question:  Stokes drift ?
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      uu(:,:,:,Krhs) = 0._wp           ! set dynamics trends to zero
+      vv(:,:,:,Krhs) = 0._wp
+      !
+!===>>>>>> Modify dyn_adv_... dyn_keg routines so that Krhs to zero useless
+      !                                         ! advection (VF or FF)	==> RHS
+      CALL dyn_adv( kstp, Kbb, Kmm      , uu, vv, Krhs, zaU, zaV, ww )
+      !                                         ! Coriolis / vorticity  ==> RHS
+      CALL dyn_vor( kstp,      Kmm      , uu, vv, Krhs )
+      !
+!!gm à appeler que pour ln_zad_Aimp=T   et en ne faisant que wi  par zdf
+!!                                              ! ZAD (implicit part)   ==> RHS
+!!    CALL dyn_zdf    ( kstp, Kbb, Kmm, Krhs, uu, vv, Kaa  )               
+
+!===>>>>>> Modify dyn_hpg & dyn_hpg_...  routines : rhd computed in dyn_hpg and pass in argument to dyn_hpg_...
+
+      CALL eos    ( ts, Kmm, rhd )              ! Kmm in situ density anomaly for hpg computation
+
+!!gm end
+      CALL dyn_hpg( kstp,      Kmm      , uu, vv, Krhs )
+      !
+!!gm ===>>>>>> Probably useless since uu_b(Kaa) will be imposed at the end of stage 1 and 2
+!                   but may be necessary in stage 3 due to implicite in dynzdf.
+!                   except if my idea for the matrice construction is OK ! 
+!      !                                         ! grad_h of ps          ==> RHS       
+!      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+!         uu(ji,jj,jk,Krhs) = uu(ji,jj,jk,Krhs) - grav * ( ssh(ji+1,jj  ,Kmm) - ssh(ji,jj,Kmm) )
+!         vv(ji,jj,jk,Krhs) = vv(ji,jj,jk,Krhs) - grav * ( ssh(ji  ,jj+1,Kmm) - ssh(ji,jj,Kmm) )
+!      END_3D
+!!gm      
+      
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! RHS of tracers : ADV only using (zaU,zaV,ww)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+#if defined key_top
+      !                       !==  Passive Tracer  ==!
+!===>>>>>> create a specific trc_RK3_adv routine with Advective velocity 
+!                                                and with Kstg as input argument
+!                         CALL trc_stp    ( kstp, Kbb, Kmm, Krhs, Kaa )            ! time-stepping
+#endif
+
+      !                       !==  T-S Tracers  ==!
+
+!===>>>>>> Modify tra_adv_...  routines so that Krhs to zero useless
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         ts(ji,jj,jk,:,Krhs) = 0._wp                                    ! set tracer trends to zero
+      END_3D
+
+!===>>> CAUTION here sans GM velocity  but with stokes drift ! 
+      CALL tra_adv( kstp, Kbb, Kmm, ts, Krhs, zaU, zaV, ww )       ! hor. + vert. advection	==> RHS
+
+
+!===>>>>>> stg1&2:  Verify the necessity of these trends (we may need it as there are in the RHS of dynspg_ts ?)
+!!gm ====>>>>   needed for heat and salt fluxes associated with mass/volume flux
+!                     CALL tra_sbc( kstp,      Kmm, ts, Krhs )   ! surface boundary condition
+!      IF( ln_isf )   CALL tra_isf( kstp,      Kmm, ts, Krhs )   ! ice shelf heat flux
+!!gm
+
+      !                                
+!!gm ===>>>>>>  Verify the necessity of these trends  at stages 1 and 2 
+!           (we may need it as they are in the RHS of dynspg_ts ?)
+      IF(  lk_asminc .AND. ln_asmiau ) THEN               ! apply assimilation increment
+         IF( ln_dyninc )   CALL dyn_asm_inc( kstp, Kbb, Kmm, uu, vv, Krhs )   ! dynamics   ==> RHS
+         IF( ln_trainc )   CALL tra_asm_inc( kstp, Kbb, Kmm, ts    , Krhs )   ! tracers    ==> RHS
+      ENDIF
+!!gm  end Verif
+
+      !
+      SELECT CASE( kstg )
+      !                    !-------------------!
+      CASE ( 1 , 2 )       !==  Stage 1 & 2  ==!   stg1:  Kbb = N  ;  Kaa = N+1/3
+         !                 !-------------------!   stg2:  Kbb = N  ;  Kmm = N+1/3  ;  Kaa = N+1/2
+         !
+         !                                      !==  time integration  ==!   ∆t = rn_Dt/3 (stg1) or rn_Dt/2 (stg2)
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            ze3Ub = e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb )
+            ze3Vb = e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb )
+            ze3Ur = e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Krhs)
+            ze3Vr = e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Krhs)
+            uu(ji,jj,jk,Kaa) = ( ze3Ub + rDt * ze3Ur*umask(ji,jj,jk) ) / e3u(ji,jj,jk, Kaa)
+            vv(ji,jj,jk,Kaa) = ( ze3Vb + rDt * ze3Vr*vmask(ji,jj,jk) ) / e3v(ji,jj,jk, Kaa)
+            !
+            ze3Tb = e3t(ji,jj,jk,Kbb) * ts(ji,jj,jk,jp_tem,Kbb )
+            ze3Sb = e3t(ji,jj,jk,Kbb) * ts(ji,jj,jk,jp_sal,Kbb )
+            ze3Tr = e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Krhs)
+            ze3Sr = e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Krhs)
+            z1_e3t= 1._wp / e3t(ji,jj,jk, Kaa)
+            ts(ji,jj,jk,jp_tem,Kaa) = ( ze3Tb + rDt * ze3Tr*tmask(ji,jj,jk) ) * z1_e3t
+            ts(ji,jj,jk,jp_sal,Kaa) = ( ze3Sb + rDt * ze3Sr*tmask(ji,jj,jk) ) * z1_e3t
+         END_3D
+         !
+         !
+         IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=uu(:,:,:,Kaa), clinfo1='stp stg   - Ua: ', mask1=umask,   &
+            &                                  tab3d_2=vv(:,:,:,Kaa), clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
+         !
+         IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Kaa), clinfo1='stp stg   - Ta: ', mask1=tmask,   &
+            &                                  tab3d_2=ts(:,:,:,jp_sal,Kaa), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
+         !
+         !
+         !                 !---------------!
+      CASE ( 3 )           !==  Stage 3  ==!   add all remaining RHS terms
+         !                 !---------------!
+         !
+         !                                      !==  complete the momentum RHS ==!   except ZDF (implicit)
+         !                                                   ! lateral mixing                    ==> RHS
+                            CALL dyn_ldf( kstp, Kbb, Kmm      , uu, vv, Krhs )
+         !                                                   ! OSMOSIS non-local velocity fluxes ==> RHS
+         IF( ln_zdfosm  )   CALL dyn_osm( kstp,      Kmm      , uu, vv, Krhs )
+         !                    
+         IF( ln_bdy     ) THEN                               ! bdy damping trends     ==> RHS
+                            CALL bdy_dyn3d_dmp ( kstp, Kbb, uu, vv, Krhs )
+                            CALL bdy_tra_dmp   ( kstp, Kbb, ts    , Krhs )
+         ENDIF
+
+#if defined key_agrif
+         IF(.NOT. Agrif_Root() ) THEN     ! AGRIF:   sponge ==> momentum and tracer RHS
+            CALL Agrif_Sponge_dyn
+            CALL Agrif_Sponge_tra
+         ENDIF
+#endif
+         !                                      !==  complete the tracers RHS  ==!   except ZDF (implicit)
+#if defined key_top
+         !                                            !*  Passive Tracer  *!
+
+!===>>>>>> create a specific trc_RK3_NOadv routine with other trend than ADV 
+!                                                and with Kstg as input argument
+!                         CALL trc_stp    ( kstp, Kbb, Kmm, Krhs, Kaa )            ! time-stepping
+#endif
+         !                                            !*  T-S Tracer  *!
+         !
+                            CALL tra_ldf( kstp, Kbb, Kmm, ts, Krhs )  ! lateral mixing
+!!gm ====>>>>   needed for heat and salt fluxes associated with mass/volume flux
+                            CALL tra_sbc( kstp,      Kmm, ts, Krhs )   ! surface boundary condition
+         IF( ln_isf )       CALL tra_isf( kstp,      Kmm, ts, Krhs )   ! ice shelf heat flux
+!!gm
+         !
+         IF( ln_traqsr  )   CALL tra_qsr( kstp,      Kmm, ts, Krhs )  ! penetrative solar radiation qsr
+         IF( ln_trabbc  )   CALL tra_bbc( kstp,      Kmm, ts, Krhs )  ! bottom heat flux
+         IF( ln_trabbl  )   CALL tra_bbl( kstp, Kbb, Kmm, ts, Krhs )  ! advective (and/or diffusive) bottom boundary layer scheme
+         IF( ln_tradmp  )   CALL tra_dmp( kstp, Kbb, Kmm, ts, Krhs )  ! internal damping trends
+
+         IF( ln_zdfmfc  )   CALL tra_mfc( kstp, Kbb,      ts, Krhs )  ! Mass Flux Convection
+         IF( ln_zdfosm  ) THEN
+                            CALL tra_osm( kstp,      Kmm, ts, Krhs )  ! OSMOSIS non-local tracer fluxes ==> RHS
+            IF( lrst_oce )  CALL osm_rst( kstp,      Kmm, 'WRITE'  )  ! write OSMOSIS outputs + ww (so must do here) to restarts
+         ENDIF
+
+         !
+         !                                      !==  DYN & TRA time integration + ZDF  ==!   ∆t = rDt
+         !
+                            CALL dyn_zdf( kstp, Kbb, Kmm, Krhs, uu, vv, Kaa  )  ! vertical diffusion and time integration 
+                            CALL tra_zdf( kstp, Kbb, Kmm, Krhs, ts    , Kaa  )  ! vertical mixing and after tracer fields
+         !
+         IF( ln_zdfnpc  )   CALL tra_npc( kstp,      Kmm, Krhs, ts    , Kaa  )  ! update after fields by non-penetrative convection
+
+
+!==>>>>>   Missing obc CALL....
+
+         !         
+      END SELECT      
+
+      !                                         !==  correction of the barotropic (all stages)  ==!    at Kaa = N+1/3, N+1/2 or N+1
+      !                                                           ! barotropic velocity correction
+      zub(:,:) = uu_b(:,:,Kaa) - SUM( e3u_0(:,:,:)*uu(:,:,:,Kaa), 3 ) * r1_hu_0(:,:)
+      zvb(:,:) = vv_b(:,:,Kaa) - SUM( e3v_0(:,:,:)*vv(:,:,:,Kaa), 3 ) * r1_hv_0(:,:)
+      DO jk = 1, jpkm1                                            ! corrected horizontal velocity
+         uu(:,:,jk,Kaa) = uu(:,:,jk,Kaa) + zub(:,:)*umask(:,:,jk)
+         vv(:,:,jk,Kaa) = vv(:,:,jk,Kaa) + zvb(:,:)*vmask(:,:,jk)
+      END DO
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Set boundary conditions
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+# if defined key_agrif
+            CALL Agrif_tra                     !* AGRIF zoom boundaries
+            CALL Agrif_dyn( kt )
+# endif
+      !                                        ! local domain boundaries  (T-point, unchanged sign)
+      CALL lbc_lnk_multi( 'stp_RK3_stg', uu(:,:,:,       Kaa), 'U', -1., vv(:,:,:       ,Kaa), 'V', -1.   &
+         &                             , ts(:,:,:,jp_tem,Kaa), 'T',  1., ts(:,:,:,jp_sal,Kaa), 'T',  1. )
+      !
+      !                                        !* BDY open boundaries
+      IF( ln_bdy )   THEN
+                               CALL bdy_tra( kstp, Kbb, ts,     Kaa )
+         IF( ln_dynspg_exp )   CALL bdy_dyn( kstp, Kbb, uu, vv, Kaa )
+         IF( ln_dynspg_ts  )   CALL bdy_dyn( kstp, Kbb, uu, vv, Kaa, dyn3d_only=.true. )
+      ENDIF
+      !
+      IF( ln_timing )   CALL timing_stop('stp_RK3_stg')
+      !
+   END SUBROUTINE stp_RK3_stg
+
+#else
+   !!----------------------------------------------------------------------
+   !!   default option             EMPTY MODULE           qco not activated
+   !!----------------------------------------------------------------------
+#endif
+   
+   !!======================================================================
+END MODULE stprk3_stg
