Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/stprk3.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/stpRK3.F90	(revision 13604)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/stprk3.F90	(revision 13769)
@@ -1,5 +1,5 @@
-MODULE stpRK3
+MODULE stprk3
    !!======================================================================
-   !!                       ***  MODULE stpRK3  ***
+   !!                       ***  MODULE stprk3  ***
    !! Time-stepping   : manager of the shallow water equation time stepping
    !!                   3rd order Runge-Kutta scheme
@@ -69,6 +69,6 @@
       INTEGER ::   indic        ! error indicator if < 0
       REAL(wp)::   z1_2rho0,  z5_6,  z3_4  ! local scalars
-      REAL(wp)::   zue3a, zue3b, zua    ! local scalars
-      REAL(wp)::   zve3a, zve3b, zva    !   -      -
+      REAL(wp)::   zue3a, zue3b, zua, zrhs_u    ! local scalars
+      REAL(wp)::   zve3a, zve3b, zva, zrhs_v    !   -      -
       !! ---------------------------------------------------------------------
       !
@@ -129,45 +129,44 @@
       CALL dyn_ldf( kstp, Nbb, Nbb, uu, vv, Nrhs )  ! lateral mixing
 #endif
-      !
-      DO_3D( 0,0, 0,0, 1,jpkm1 )
+!!st       !
+!!st       DO_3D( 0,0, 0,0, 1,jpkm1 )
+!!st          !                                          ! horizontal pressure gradient
+!!st          uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nbb) - ssh(ji,jj,Nbb) ) * r1_e1u(ji,jj)
+!!st          vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nbb) - ssh(ji,jj,Nbb) ) * r1_e2v(ji,jj)
+!!st       END_3D
+!!st       !
+!!st #if defined key_RK3all 
+!!st       !                                             ! wind stress and layer friction
+!!st       z5_6 = 5._wp/6._wp
+!!st       DO_3D( 0, 0, 0, 0,1,jpkm1)
+!!st          uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + r1_rho0 * ( z5_6*utau_b(ji,jj) + (1._wp - z5_6)*utau(ji,jj) ) / e3u(ji,jj,jk,Nbb)   &
+!!st             &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
+!!st          vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + r1_rho0 * ( z5_6*vtau_b(ji,jj) + (1._wp - z5_6)*vtau(ji,jj) ) / e3v(ji,jj,jk,Nbb)   &
+!!st             &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
+!!st       END_3D
+!!st #endif
+!!st why not ?
+      z5_6 = 5._wp/6._wp
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          !                                          ! horizontal pressure gradient
-         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nbb) - ssh(ji,jj,Nbb) ) * r1_e1u(ji,jj)
-         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nbb) - ssh(ji,jj,Nbb) ) * r1_e2v(ji,jj)
+         zrhs_u =        - grav    * ( ssh(ji+1,jj,Nbb) - ssh(ji,jj,Nbb) ) * r1_e1u(ji,jj)
+         zrhs_v =        - grav    * ( ssh(ji,jj+1,Nbb) - ssh(ji,jj,Nbb) ) * r1_e2v(ji,jj)
+#if defined key_RK3all
+         !                                          ! wind stress and layer friction
+         zrhs_u = zrhs_u + r1_rho0 * ( z5_6*utau_b(ji,jj) + (1._wp - z5_6)*utau(ji,jj) ) / e3u(ji,jj,jk,Nbb)   &
+            &            - rn_rfr  * uu(ji,jj,jk,Nbb)
+         zrhs_v = zrhs_v + r1_rho0 * ( z5_6*vtau_b(ji,jj) + (1._wp - z5_6)*vtau(ji,jj) ) / e3v(ji,jj,jk,Nbb)   &
+            &            - rn_rfr  * vv(ji,jj,jk,Nbb)
+#endif
+         !                                          ! ==> RHS
+         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + zrhs_u
+         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + zrhs_v
       END_3D
-      !
-#if defined key_RK3all 
-      !                                             ! wind stress and layer friction
-      z5_6 = 5._wp/6._wp
-      DO_3D( 0, 0, 0, 0,1,jpkm1)
-         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + r1_rho0 * ( z5_6*utau_b(ji,jj) + (1._wp - z5_6)*utau(ji,jj) ) / e3u(ji,jj,jk,Nbb)   &
-            &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
-         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + r1_rho0 * ( z5_6*vtau_b(ji,jj) + (1._wp - z5_6)*vtau(ji,jj) ) / e3v(ji,jj,jk,Nbb)   &
-            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
-      END_3D
-#endif
-!!!!st why not ?
-!!      z5_6 = 5._wp/6._wp
-!!      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-!!         !                                                              ! horizontal pressure gradient
-!!         zrhs_u =        - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
-!!         zrhs_v =        - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
-!!         !                                                              ! wind stress and layer friction
-!!#if defined key_RK3all
-!!         zrhs_u = zrhs_u + r1_rho0 * ( z5_6*utau_b(ji,jj) + (1._wp - z5_6)*utau(ji,jj) ) / e3u(ji,jj,jk,Nbb)   &
-!!            &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
-!!         zrhs_v = zrhs_v + r1_rho0 * ( z5_6*vtau_b(ji,jj) + (1._wp - z5_6)*vtau(ji,jj) ) / e3v(ji,jj,jk,Nbb)   &
-!!            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
-!!         !                                                              ! ==> RHS
-!!#endif
-
-!!         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + zrhs_u
-!!         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + zrhs_v
-!!      END_3D
-!!!!st end
+!!st end
       !
       !                                 !==  Time stepping of ssh Eq.  ==!   (and update r3_Naa)
       !
       CALL ssh_nxt( kstp, Nbb, Nbb, ssh, Naa )      ! after ssh
-      !                                                                 ! after ssh/h_0 ratio explicit
+      !                                             ! after ssh/h_0 ratio
       CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )
       !
@@ -175,10 +174,10 @@
       !
       IF( ln_dynadv_vec ) THEN                      ! vector invariant form : applied on velocity
-         DO_3D( 0, 0, 0, 0,1,jpkm1)
+         DO_3D( 0, 0, 0, 0, 1,jpkm1)
             uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
             vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
          END_3D          
       ELSE
-         DO_3D( 0, 0, 0, 0,1,jpkm1)                 ! flux form : applied on thickness weighted velocity
+         DO_3D( 0, 0, 0, 0, 1,jpkm1)                 ! flux form : applied on thickness weighted velocity
             zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
             zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
@@ -191,6 +190,7 @@
       ENDIF
       !
+      CALL lbc_lnk_multi( 'stp_RK3', uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1. )
+      !
       !                                 !==  Swap time levels  ==!
-      !
       Nrhs= Nnn
       Nnn = Naa
@@ -213,25 +213,42 @@
 #endif
       !
+      z3_4 = 3._wp/4._wp
       DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          !                                          ! horizontal pressure gradient
-         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
-         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
+         zrhs_u =        - grav    * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
+         zrhs_v =        - grav    * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
+#if defined key_RK3all
+         !                                          ! wind stress and layer friction
+         zrhs_u = zrhs_u + r1_rho0 * ( z3_4*utau_b(ji,jj) + (1._wp - z3_4)*utau(ji,jj) ) / e3u(ji,jj,jk,Nnn)   &
+            &            - rn_rfr  * uu(ji,jj,jk,Nbb)
+         zrhs_v = zrhs_v + r1_rho0 * ( z3_4*vtau_b(ji,jj) + (1._wp - z3_4)*vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   &
+            &            - rn_rfr  * vv(ji,jj,jk,Nbb)
+#endif
+         !                                          ! ==> RHS
+         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + zrhs_u
+         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + zrhs_v
       END_3D
-      !
-#if defined key_RK3all
-      !                                             ! wind stress and layer friction
-      z3_4 = 3._wp/4._wp
-      DO_3D( 0, 0, 0, 0,1,jpkm1)
-         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + r1_rho0 * ( z3_4*utau_b(ji,jj) + (1._wp - z3_4)*utau(ji,jj) ) / e3u(ji,jj,jk,Nbb)   &
-            &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
-         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + r1_rho0 * ( z3_4*vtau_b(ji,jj) + (1._wp - z3_4)*vtau(ji,jj) ) / e3v(ji,jj,jk,Nbb)   &
-            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
-      END_3D
-#endif
+!!st       !
+!!st       DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+!!st          !                                          ! horizontal pressure gradient
+!!st          uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
+!!st          vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
+!!st       END_3D
+!!st       !
+!!st #if defined key_RK3all
+!!st       !                                             ! wind stress and layer friction
+!!st       z3_4 = 3._wp/4._wp
+!!st       DO_3D( 0, 0, 0, 0,1,jpkm1)
+!!st          uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + r1_rho0 * ( z3_4*utau_b(ji,jj) + (1._wp - z3_4)*utau(ji,jj) ) / e3u(ji,jj,jk,Nbb)   &
+!!st             &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
+!!st          vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + r1_rho0 * ( z3_4*vtau_b(ji,jj) + (1._wp - z3_4)*vtau(ji,jj) ) / e3v(ji,jj,jk,Nbb)   &
+!!st             &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
+!!st       END_3D
+!!st #endif
       !
       !                                 !==  Time stepping of ssh Eq.  ==!   (and update r3_Naa)
       !
       CALL ssh_nxt( kstp, Nbb, Nnn, ssh, Naa )      ! after ssh
-      !                                                                 ! after ssh/h_0 ratio explicit
+      !                                             ! after ssh/h_0 ratio
       CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )
       !
@@ -239,10 +256,10 @@
       !
       IF( ln_dynadv_vec ) THEN                      ! vector invariant form : applied on velocity
-         DO_3D( 0, 0, 0, 0,1,jpkm1)
+         DO_3D( 0, 0, 0, 0, 1,jpkm1)
             uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
             vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
          END_3D          
       ELSE
-         DO_3D( 0, 0, 0, 0,1,jpkm1)                 ! flux form : applied on thickness weighted velocity
+         DO_3D( 0, 0, 0, 0, 1,jpkm1)                 ! flux form : applied on thickness weighted velocity
             zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
             zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
@@ -255,6 +272,7 @@
       ENDIF
       !
+      CALL lbc_lnk_multi( 'stp_RK3', uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1. )
+      !
       !                                 !==  Swap time levels  ==!
-      !
       Nrhs= Nnn
       Nnn = Naa
@@ -276,22 +294,23 @@
       CALL dyn_ldf( kstp, Nbb, Nnn, uu, vv, Nrhs )  ! lateral mixing
 
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+      z1_2rho0 = 0.5_wp * r1_rho0
+      DO_3D( 0, 0, 0, 0, 1,jpkm1 )
          !                                          ! horizontal pressure gradient
-         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
-         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
+         zrhs_u =        - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
+         zrhs_v =        - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
+         !                                          ! wind stress and layer friction
+         zrhs_u = zrhs_u + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn)   &
+            &            - rn_rfr   * uu(ji,jj,jk,Nbb)
+         zrhs_v = zrhs_v + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   &
+            &            - rn_rfr   * vv(ji,jj,jk,Nbb)
+         !                                          ! ==> RHS
+         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + zrhs_u
+         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + zrhs_v
       END_3D
-      !                                             ! wind stress and layer friction
-      z1_2rho0 = 0.5_wp * r1_rho0
-      DO_3D( 0, 0, 0, 0,1,jpkm1)
-         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn)   &
-            &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
-         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   &
-            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
-      END_3D
       !
       !                                 !==  Time stepping of ssh Eq.  ==!   (and update r3_Naa)
       !
       CALL ssh_nxt( kstp, Nbb, Nnn, ssh, Naa )      ! after ssh
-      !                                                                 ! after ssh/h_0 ratio explicit
+      !                                             ! after ssh/h_0 ratio
       CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )
       !
@@ -299,17 +318,11 @@
       !
       IF( ln_dynadv_vec ) THEN                      ! vector invariant form : applied on velocity
-         DO_3D( 1, 1, 1, 1,1,jpkm1)
-            zua = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
-            zva = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
-            !                                                                  ! Asselin time filter on u,v (Nnn)
-            uu(ji,jj,jk,Nnn) = uu(ji,jj,jk,Nnn) + rn_atfp * (uu(ji,jj,jk,Nbb) - 2._wp * uu(ji,jj,jk,Nnn) + zua)
-            vv(ji,jj,jk,Nnn) = vv(ji,jj,jk,Nnn) + rn_atfp * (vv(ji,jj,jk,Nbb) - 2._wp * vv(ji,jj,jk,Nnn) + zva)
-            !              
-            uu(ji,jj,jk,Naa) = zua
-            vv(ji,jj,jk,Naa) = zva
+         DO_3D( 0, 0, 0, 0, 1,jpkm1)
+            uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
+            vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
          END_3D
          !
       ELSE                                          ! flux form : applied on thickness weighted velocity
-         DO_3D( 1, 1, 1, 1,1,jpkm1)
+         DO_3D( 0, 0, 0, 0, 1,jpkm1)
             zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
             zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
@@ -321,7 +334,6 @@
          END_3D
       ENDIF
-
-      CALL lbc_lnk_multi( 'stp_RK3', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1.,   &   !* local domain boundaries
-         &                           uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1.    )     
+      !
+      CALL lbc_lnk_multi( 'stp_RK3', uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1. )
       !
       !                                 !==  Swap time levels  ==!
@@ -363,3 +375,3 @@
 
    !!======================================================================
-END MODULE stpRK3
+END MODULE stprk3
