Index: /NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DOM/domqco.F90
===================================================================
--- /NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DOM/domqco.F90	(revision 15626)
+++ /NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DOM/domqco.F90	(revision 15627)
@@ -39,5 +39,6 @@
    PUBLIC  dom_qco_init       ! called by domain.F90
    PUBLIC  dom_qco_zgr        ! called by isfcpl.F90
-   PUBLIC  dom_qco_r3c        ! called by steplf.F90
+   PUBLIC  dom_qco_r3c        ! called by stpmlf.F90
+   PUBLIC  dom_qco_r3c_RK3    ! called by stprk3_stg.F90
 
    !                                                      !!* Namelist nam_vvl
@@ -118,5 +119,5 @@
 #if defined key_RK3
       CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb), r3f(:,:) )
-      CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm)           )
+      CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm)           )  !! needed for agrif ???
 #else
       CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb)           )
@@ -156,7 +157,4 @@
       !                                      !==  ratio at u-,v-point  ==!
       !
-!!st      IF( ln_dynadv_vec ) THEN                     !- Vector Form   (thickness weighted averaging)
-#if ! defined key_qcoTest_FluxForm
-      !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
       DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
          pr3u(ji,jj) = 0.5_wp * (  e1e2t(ji  ,jj) * pssh(ji  ,jj)  &
@@ -164,8 +162,67 @@
          pr3v(ji,jj) = 0.5_wp * (  e1e2t(ji,jj  ) * pssh(ji,jj  )  &
             &                    + e1e2t(ji,jj+1) * pssh(ji,jj+1)  ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj)
+      END_2D      
+      !
+      IF( .NOT.PRESENT( pr3f ) ) THEN              !- lbc on ratio at u-, v-points only
+         IF (nn_hls==1) CALL lbc_lnk( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp )
+         !
+         !
+      ELSE                                   !==  ratio at f-point  ==!
+         !
+         DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
+            ! round brackets added to fix the order of floating point operations
+            ! needed to ensure halo 1 - halo 2 compatibility
+            pr3f(ji,jj) = 0.25_wp * (   (  e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )      &
+               &                         + e1e2t(ji+1,jj  ) * pssh(ji+1,jj  )  )   & ! bracket for halo 1 - halo 2 compatibility
+               &                     +  (  e1e2t(ji  ,jj+1) * pssh(ji  ,jj+1)      &
+               &                         + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  )   & ! bracket for halo 1 - halo 2 compatibility
+               &                    ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj)
+         END_2D
+         !                                                 ! lbc on ratio at u-,v-,f-points
+         IF (nn_hls==1) CALL lbc_lnk( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp, pr3f, 'F', 1._wp )
+         !
+      ENDIF
+      !
+   END SUBROUTINE dom_qco_r3c
+
+   
+   SUBROUTINE dom_qco_r3c_RK3( pssh, pr3t, pr3u, pr3v, pr3f )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE r3c  ***
+      !!
+      !! ** Purpose :   compute the filtered ratio ssh/h_0 at t-,u-,v-,f-points
+      !!
+      !! ** Method  : - compute the ssh at u- and v-points (f-point optional)
+      !!                   Vector Form : surface weighted averaging
+      !!                   Flux   Form : simple           averaging
+      !!              - compute the ratio ssh/h_0 at t-,u-,v-pts, (f-pt optional)
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:)          , INTENT(in   )  ::   pssh               ! sea surface height   [m]
+      REAL(wp), DIMENSION(:,:)          , INTENT(  out)  ::   pr3t, pr3u, pr3v   ! ssh/h0 ratio at t-, u-, v-,points  [-]
+      REAL(wp), DIMENSION(:,:), OPTIONAL, INTENT(  out)  ::   pr3f               ! ssh/h0 ratio at f-point   [-]
+      !
+      INTEGER ::   ji, jj   ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      !
+      DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
+         pr3t(ji,jj) = pssh(ji,jj) * r1_ht_0(ji,jj)   !==  ratio at t-point  ==!
+      END_2D
+      !
+      !
+      !                                      !==  ratio at u-,v-point  ==!
+      !
+!!st      IF( ln_dynadv_vec ) THEN                     !- Vector Form   (thickness weighted averaging)
+#if ! defined key_qcoTest_FluxForm
+      !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
+      DO_2D( 0, 0, 0, 0 )
+         pr3u(ji,jj) = 0.5_wp * (  e1e2t(ji  ,jj) * pssh(ji  ,jj)  &
+            &                    + e1e2t(ji+1,jj) * pssh(ji+1,jj)  ) * r1_hu_0(ji,jj) * r1_e1e2u(ji,jj)
+         pr3v(ji,jj) = 0.5_wp * (  e1e2t(ji,jj  ) * pssh(ji,jj  )  &
+            &                    + e1e2t(ji,jj+1) * pssh(ji,jj+1)  ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj)
       END_2D
 !!st      ELSE                                         !- Flux Form   (simple averaging)
 #else
-      DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
+      DO_2D( 0, 0, 0, 0 )
          pr3u(ji,jj) = 0.5_wp * (  pssh(ji,jj) + pssh(ji+1,jj  )  ) * r1_hu_0(ji,jj)
          pr3v(ji,jj) = 0.5_wp * (  pssh(ji,jj) + pssh(ji  ,jj+1)  ) * r1_hv_0(ji,jj)
@@ -174,9 +231,5 @@
 #endif         
       !
-      IF( .NOT.PRESENT( pr3f ) ) THEN              !- lbc on ratio at u-, v-points only
-         IF (nn_hls==1) CALL lbc_lnk( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp )
-         !
-         !
-      ELSE                                   !==  ratio at f-point  ==!
+      IF( PRESENT( pr3f ) ) THEN             !==  ratio at f-point  ==!
          !
 !!st         IF( ln_dynadv_vec )   THEN                !- Vector Form   (thickness weighted averaging)
@@ -184,35 +237,30 @@
          !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
 
-      DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
+      DO_2D( 0, 0, 0, 0 )
          ! round brackets added to fix the order of floating point operations
          ! needed to ensure halo 1 - halo 2 compatibility
-         pr3f(ji,jj) = 0.25_wp * ( ( e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )   &
-            &                      + e1e2t(ji+1,jj  ) * pssh(ji+1,jj  )   &
-            &                      )                                      & ! bracket for halo 1 - halo 2 compatibility
-            &                     + ( e1e2t(ji  ,jj+1) * pssh(ji  ,jj+1)  &
-            &                       + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  &
-            &                       )                                     & ! bracket for halo 1 - halo 2 compatibility
+         pr3f(ji,jj) = 0.25_wp * (   (  e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )      &
+            &                         + e1e2t(ji+1,jj  ) * pssh(ji+1,jj  )  )   & ! bracket for halo 1 - halo 2 compatibility
+            &                     +  (  e1e2t(ji  ,jj+1) * pssh(ji  ,jj+1)      &
+            &                         + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  )   & ! bracket for halo 1 - halo 2 compatibility
             &                    ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj)
       END_2D
 !!st         ELSE                                      !- Flux Form   (simple averaging)
 #else
-      DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
+      DO_2D( 0, 0, 0, 0 )
          ! round brackets added to fix the order of floating point operations
          ! needed to ensure halo 1 - halo 2 compatibility
-         pr3f(ji,jj) = 0.25_wp * ( ( pssh(ji,jj  ) + pssh(ji+1,jj  ) ) &
-            &                    + ( pssh(ji,jj+1) + pssh(ji+1,jj+1)   &
-            &                       )                                  & ! bracket for halo 1 - halo 2 compatibility
+         pr3f(ji,jj) = 0.25_wp * (   (  pssh(ji,jj  ) + pssh(ji+1,jj  )  )   &
+            &                     +  (  pssh(ji,jj+1) + pssh(ji+1,jj+1)  )   & ! bracket for halo 1 - halo 2 compatibility
             &                    ) * r1_hf_0(ji,jj)
       END_2D
 !!st         ENDIF
 #endif
-         !                                                 ! lbc on ratio at u-,v-,f-points
-         IF (nn_hls==1) CALL lbc_lnk( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp, pr3f, 'F', 1._wp )
          !
       ENDIF
       !
-   END SUBROUTINE dom_qco_r3c
-
-
+   END SUBROUTINE dom_qco_r3c_RK3
+
+    
    SUBROUTINE qco_ctl
       !!---------------------------------------------------------------------
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 15626)
+++ /NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stprk3_stg.F90	(revision 15627)
@@ -49,6 +49,7 @@
    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
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssha                    ! sea-surface height  at N+1
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   r3ta, r3ua, r3va, r3fa, r3fb   ! ssh/h_0 ratio at t,u,v-column
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ua_b, va_b              ! barotropic velocity at N+1
 
    !! * Substitutions
@@ -120,4 +121,22 @@
          vv_b(:,:,Kaa) = r2_3 * vv_b(:,:,Kbb) + r1_3 * va_b(:,:)
          !
+         !
+         !                     !==  ssh/h0 ratio at Kaa  ==! 
+         !
+         IF( .NOT.lk_linssh ) THEN     ! "after" ssh/h_0 ratio at t,u,v-column
+            !
+            ALLOCATE( r3ta(jpi,jpj) , r3ua(jpi,jpj) , r3va(jpi,jpj) , r3fa(jpi,jpj) , r3fb(jpi,jpj) )
+            !
+            r3fb(:,:) = r3f(:,:)       !!st dirty fix check with gm
+            CALL dom_qco_r3c_RK3( ssha, r3ta, r3ua, r3va, r3fa )
+            !
+            CALL lbc_lnk( 'stprk3_stg', r3ua, 'U', 1._wp, r3va, 'V', 1._wp, r3fa, 'F', 1._wp )
+            !
+            r3t(:,:,Kaa) = r2_3 * r3t(:,:,Kbb) + r1_3 * r3ta(:,:)
+            r3u(:,:,Kaa) = r2_3 * r3u(:,:,Kbb) + r1_3 * r3ua(:,:)
+            r3v(:,:,Kaa) = r2_3 * r3v(:,:,Kbb) + r1_3 * r3va(:,:)
+            r3f(:,:) = r2_3 * r3fb(:,:) + r1_3 * r3fa(:,:)
+         ENDIF
+         !
          !                 !---------------!
       CASE ( 2 )           !==  Stage 2  ==!   Kbb = N   ;   Kmm = N+1/3   ;   Kaa = N+1/2
@@ -132,4 +151,11 @@
          vv_b(:,:,Kaa) = r1_2 * ( vv_b(:,:,Kbb) + va_b(:,:) )
          !
+         IF( .NOT.lk_linssh ) THEN
+            r3t(:,:,Kaa) = r1_2 * ( r3t(:,:,Kbb) + r3ta(:,:) )
+            r3u(:,:,Kaa) = r1_2 * ( r3u(:,:,Kbb) + r3ua(:,:) )
+            r3v(:,:,Kaa) = r1_2 * ( r3v(:,:,Kbb) + r3va(:,:) )
+            r3f(:,:) = r1_2 * ( r3fb(:,:) + r3fa(:,:) )
+         ENDIF
+         !
          !                 !---------------!
       CASE ( 3 )           !==  Stage 3  ==!   Kbb = N   ;   Kmm = N+1/2   ;   Kaa = N+1
@@ -145,16 +171,14 @@
          DEALLOCATE( ssha , ua_b , va_b )
          !
+         IF( .NOT.lk_linssh ) THEN
+            r3t(:,:,Kaa) = r3ta(:,:)
+            r3u(:,:,Kaa) = r3ua(:,:)
+            r3v(:,:,Kaa) = r3va(:,:)
+            r3f(:,:    ) = r3fa(:,:)
+            DEALLOCATE( r3ta, r3ua, r3va, r3fa, r3fb )
+            !
+         ENDIF
+         !
       END SELECT
-      !
-      !                     !==  ssh/h0 ratio at Kaa  ==! 
-      !
-      IF( .NOT.lk_linssh )   CALL dom_qco_r3c( ssh(:,:,Kaa), r3t(:,:,Kaa), r3u(:,:,Kaa), r3v(:,:,Kaa), r3f(:,:) )   ! "after" ssh/h_0 ratio at t,u,v-column
-!!      SELECT CASE( kstg )
-         !
-!!      CASE ( 3 )
-         !!st required at each stage for div hor loops
-         CALL lbc_lnk( 'stprk3_stg', r3u(:,:,Kaa), 'U', 1._wp, r3v(:,:,Kaa), 'V', 1._wp, r3f(:,:), 'F', 1._wp )
-         !
-!!      END SELECT
       !
       !
