Index: NEMO/branches/2021/dev_r14447_HPC-07_Irrmann_try_new_pt2pt/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/2021/dev_r14447_HPC-07_Irrmann_try_new_pt2pt/src/OCE/DYN/dynspg_ts.F90	(revision 15044)
+++ NEMO/branches/2021/dev_r14447_HPC-07_Irrmann_try_new_pt2pt/src/OCE/DYN/dynspg_ts.F90	(revision 15054)
@@ -517,4 +517,5 @@
          ENDIF    
          !
+         !                    In loop_ssha_e
          !
          !     Compute Sea Level at step jit+1
@@ -522,13 +523,9 @@
          !--        ssh    =  ssh   - delta_t' * [ frc + div( flux      ) ]      --!
          !-------------------------------------------------------------------------!
-
          !
          IF( ln_async )  THEN
-            CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T', 1._wp,  zhU, 'U', -1._wp,  zhV, 'V', -1._wp, loop_fct=loop_fct1 )
+            CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T', 1._wp,  zhU, 'U', -1._wp,  zhV, 'V', -1._wp, loop_fct=loop_ssha_e )
          ELSE
-            DO_2D( 0, 0, 0, 0 )
-                zhdiv = (   zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1)   ) * r1_e1e2t(ji,jj)
-                ssha_e(ji,jj) = (  sshn_e(ji,jj) - rDt_e * ( zssh_frc(ji,jj) + zhdiv )  ) * ssmask(ji,jj)
-            END_2D
+            CALL loop_ssha_e( 2, jpi-1, 2, jpj-1, 1, jpkm1 )   ! arguments are useless in that case
             CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T', 1._wp,  zhU, 'U', -1._wp,  zhV, 'V', -1._wp )
          END IF
@@ -614,4 +611,6 @@
          ENDIF
          !
+         !                               In loop_velocity
+         !
          ! Set next velocities:
          !     Compute barotropic speeds at step jit+1    (h : total height of the water colomn)
@@ -620,11 +619,221 @@
          !--  u     =             u   + delta_t' * \         (1-r)*g * grad_x( ssh') -         f * k vect u      +     frc /    --!
          !--                                                                                                                    --!
-         !--                             FLUX FORM                                                                              --!
+         !--                                FLUX FORM                                                                           --!
          !--  m+1   __1__  /  m    m               /  m+1/2                             m+1/2              m+1/2    n      \ \  --!
          !-- u    =   m+1 |  h  * u   + delta_t' * \ h     * (1-r)*g * grad_x( ssh') - h     * f * k vect u      + h * frc /  | --!
          !--         h     \                                                                                                 /  --!
          !------------------------------------------------------------------------------------------------------------------------!
-         IF( ln_dynadv_vec .OR. ln_linssh ) THEN      !* Vector form
-            DO_2D( 0, 0, 0, 0 )
+
+         IF( ln_async )  THEN
+            IF( .NOT.ln_linssh ) THEN   !* Update ocean depth (variable volume case only)
+               CALL lbc_lnk( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  &
+                    &                   , hu_e , 'U',  1._wp, hv_e , 'V',  1._wp  &
+                    &                   , hur_e, 'U',  1._wp, hvr_e, 'V',  1._wp, loop_fct=loop_velocity  )
+            ELSE
+               CALL lbc_lnk( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp, loop_fct=loop_velocity  )
+            ENDIF
+         ELSE
+            CALL loop_velocity( 2, jpi-1, 2, jpj-1, 1, jpkm1 )   ! arguments are useless in that case
+            IF( .NOT.ln_linssh ) THEN   !* Update ocean depth (variable volume case only)
+               CALL lbc_lnk( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  &
+                    &                   , hu_e , 'U',  1._wp, hv_e , 'V',  1._wp  &
+                    &                   , hur_e, 'U',  1._wp, hvr_e, 'V',  1._wp  )
+            ELSE
+               CALL lbc_lnk( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  )
+            ENDIF
+         END IF
+         
+         !                                                 ! open boundaries
+         IF( ln_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
+#if defined key_agrif                                                           
+         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
+#endif
+         !                                             !* Swap
+         !                                             !  ----
+         ubb_e  (:,:) = ub_e  (:,:)
+         ub_e   (:,:) = un_e  (:,:)
+         un_e   (:,:) = ua_e  (:,:)
+         !
+         vbb_e  (:,:) = vb_e  (:,:)
+         vb_e   (:,:) = vn_e  (:,:)
+         vn_e   (:,:) = va_e  (:,:)
+         !
+         sshbb_e(:,:) = sshb_e(:,:)
+         sshb_e (:,:) = sshn_e(:,:)
+         sshn_e (:,:) = ssha_e(:,:)
+
+         !                                             !* Sum over whole bt loop
+         !                                             !  ----------------------
+         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(:,:)
+
+         !                                                 ! ==================== !
+      END DO                                               !        end loop      !
+      !                                                    ! ==================== !
+      lints = .FALSE.
+      ! -----------------------------------------------------------------------------
+      ! Phase 3. update the general trend with the barotropic trend
+      ! -----------------------------------------------------------------------------
+      !
+      ! Set advection velocity correction:
+      IF (ln_bt_fw) THEN
+         IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN
+            DO_2D( 1, 1, 1, 1 )
+               zun_save = un_adv(ji,jj)
+               zvn_save = vn_adv(ji,jj)
+               !                          ! apply the previously computed correction 
+               un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - rn_atfp * un_bf(ji,jj) )
+               vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - rn_atfp * vn_bf(ji,jj) )
+               !                          ! Update corrective fluxes for next time step
+               un_bf(ji,jj)  = rn_atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )
+               vn_bf(ji,jj)  = rn_atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )
+               !                          ! Save integrated transport for next computation
+               ub2_b(ji,jj) = zun_save
+               vb2_b(ji,jj) = zvn_save
+            END_2D
+         ELSE
+            un_bf(:,:) = 0._wp            ! corrective fluxes for next time step set to zero
+            vn_bf(:,:) = 0._wp
+            ub2_b(:,:) = un_adv(:,:)      ! Save integrated transport for next computation
+            vb2_b(:,:) = vn_adv(:,:)
+         END IF
+      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
+         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( '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(:,:) )
+      ENDIF
+
+
+      ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)  
+      DO jk = 1, jpkm1
+         puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)*r1_hu(:,:,Kmm) - puu_b(:,:,Kmm) ) * umask(:,:,jk)
+         pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) + vn_adv(:,:)*r1_hv(:,:,Kmm) - pvv_b(:,:,Kmm) ) * vmask(:,:,jk)
+      END DO
+
+      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) 
+            pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) & 
+                       & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk)  
+         END DO
+      END IF 
+
+      
+      CALL iom_put(  "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) )    ! barotropic i-current
+      CALL iom_put(  "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) )    ! barotropic i-current
+      !
+#if defined key_agrif
+      ! Save time integrated fluxes during child grid integration
+      ! (used to update coarse grid transports at next time step)
+      !
+      IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
+         IF( Agrif_NbStepint() == 0 ) THEN
+            ub2_i_b(:,:) = 0._wp
+            vb2_i_b(:,:) = 0._wp
+         END IF
+         !
+         za1 = 1._wp / REAL(Agrif_rhot(), wp)
+         ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
+         vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
+      ENDIF
+#endif      
+      !                                   !* write time-spliting arrays in the restart
+      IF( lrst_oce .AND.ln_bt_fw )   CALL ts_rst( kt, 'WRITE' )
+      !
+      IF( ln_wd_il )   DEALLOCATE( zcpx, zcpy )
+      IF( ln_wd_dl )   DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
+      !
+      CALL iom_put( "baro_u" , puu_b(:,:,Kmm) )  ! Barotropic  U Velocity
+      CALL iom_put( "baro_v" , pvv_b(:,:,Kmm) )  ! Barotropic  V Velocity
+      !
+      !
+    CONTAINS
+      SUBROUTINE loop_ssha_e(i0, i1, j0, j1, k0, k1, buf)
+        !!---------------------------------------------------------------------
+        !!                   ***  SUBROUTINE loop_ssha_e  ***
+        !!
+        !! ** Purpose : Set ssha_e for next sub time step
+        !!----------------------------------------------------------------------
+        INTEGER, INTENT(in) :: i0, i1, j0, j1, k0, k1
+        REAL*8, DIMENSION(:,:,:,:,:,:), OPTIONAL, INTENT(out) :: buf
+        !     Compute Sea Level at step jit+1
+        !--           m+1        m                               m+1/2          --!
+        !--        ssh    =  ssh   - delta_t' * [ frc + div( flux      ) ]      --!
+        !-------------------------------------------------------------------------!
+        DO_2D( 0, 0, 0, 0 )
+            zhdiv = (   zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1)   ) * r1_e1e2t(ji,jj)
+            ssha_e(ji,jj) = (  sshn_e(ji,jj) - rDt_e * ( zssh_frc(ji,jj) + zhdiv )  ) * ssmask(ji,jj)
+        END_2D
+      END SUBROUTINE loop_ssha_e
+      !
+      SUBROUTINE loop_velocity(i0, i1, j0, j1, k0, k1, buf)
+        !!---------------------------------------------------------------------
+        !!                   ***  SUBROUTINE loop_velocity  ***
+        !!
+        !! ** Purpose : Set velocities for next sub time step
+        !!----------------------------------------------------------------------
+        INTEGER, INTENT(in) :: i0, i1, j0, j1, k0, k1
+        REAL*8, DIMENSION(:,:,:,:,:,:), OPTIONAL, INTENT(out) :: buf
+        !
+        ! Set next velocities:
+        !     Compute barotropic speeds at step jit+1    (h : total height of the water colomn)
+        !--                              VECTOR FORM
+        !--   m+1                 m               /                                                       m+1/2           \    --!
+        !--  u     =             u   + delta_t' * \         (1-r)*g * grad_x( ssh') -         f * k vect u      +     frc /    --!
+        !--                                                                                                                    --!
+        !--                             FLUX FORM                                                                              --!
+        !--  m+1   __1__  /  m    m               /  m+1/2                             m+1/2              m+1/2    n      \ \  --!
+        !-- u    =   m+1 |  h  * u   + delta_t' * \ h     * (1-r)*g * grad_x( ssh') - h     * f * k vect u      + h * frc /  | --!
+        !--         h     \                                                                                                 /  --!
+        !------------------------------------------------------------------------------------------------------------------------!
+        IF( ln_dynadv_vec .OR. ln_linssh ) THEN      !* Vector form
+           DO_2D( 0, 0, 0, 0 )
                ua_e(ji,jj) = (                                 un_e(ji,jj)   & 
                          &     + rDt_e * (                   zu_spg(ji,jj)   &
@@ -638,8 +847,8 @@
                          &                                 + zv_frc(ji,jj) ) &
                          &   ) * ssvmask(ji,jj)
-            END_2D
-            !
-         ELSE                           !* Flux form
-            DO_2D( 0, 0, 0, 0 )
+           END_2D
+           !
+        ELSE                           !* Flux form
+           DO_2D( 0, 0, 0, 0 )
                !                    ! hu_e, hv_e hold depth at jn,  zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2
                !                    ! backward interpolated depth used in spg terms at jn+1/2
@@ -667,187 +876,23 @@
                     &                       + zhvp2_e(ji,jj) * zv_trd (ji,jj)  &   !
                     &                       +  hv(ji,jj,Kmm) * zv_frc (ji,jj)  )   ) * z1_hv
-            END_2D
-         ENDIF
+           END_2D
+        ENDIF
 !jth implicit bottom friction:
-         IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
-            DO_2D( 0, 0, 0, 0 )
-               ua_e(ji,jj) =  ua_e(ji,jj) / ( 1._wp - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj) )
-               va_e(ji,jj) =  va_e(ji,jj) / ( 1._wp - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj) )
-            END_2D
-         ENDIF
+        IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
+           DO_2D( 0, 0, 0, 0 )
+              ua_e(ji,jj) =  ua_e(ji,jj) / ( 1._wp - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj) )
+              va_e(ji,jj) =  va_e(ji,jj) / ( 1._wp - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj) )
+           END_2D
+        ENDIF
        
-         IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
-            hu_e (2:jpim1,2:jpjm1) =    hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
-            hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / (  hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1)  )
-            hv_e (2:jpim1,2:jpjm1) =    hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)
-            hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / (  hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1)  )
-         ENDIF
-         !
-         IF( .NOT.ln_linssh ) THEN   !* Update ocean depth (variable volume case only)
-            CALL lbc_lnk( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  &
-                 &                   , hu_e , 'U',  1._wp, hv_e , 'V',  1._wp  &
-                 &                   , hur_e, 'U',  1._wp, hvr_e, 'V',  1._wp  )
-         ELSE
-            CALL lbc_lnk( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  )
-         ENDIF
-         !                                                 ! open boundaries
-         IF( ln_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
-#if defined key_agrif                                                           
-         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
-#endif
-         !                                             !* Swap
-         !                                             !  ----
-         ubb_e  (:,:) = ub_e  (:,:)
-         ub_e   (:,:) = un_e  (:,:)
-         un_e   (:,:) = ua_e  (:,:)
-         !
-         vbb_e  (:,:) = vb_e  (:,:)
-         vb_e   (:,:) = vn_e  (:,:)
-         vn_e   (:,:) = va_e  (:,:)
-         !
-         sshbb_e(:,:) = sshb_e(:,:)
-         sshb_e (:,:) = sshn_e(:,:)
-         sshn_e (:,:) = ssha_e(:,:)
-
-         !                                             !* Sum over whole bt loop
-         !                                             !  ----------------------
-         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(:,:)
-
-         !                                                 ! ==================== !
-      END DO                                               !        end loop      !
-      !                                                    ! ==================== !
-      lints = .FALSE.
-      ! -----------------------------------------------------------------------------
-      ! Phase 3. update the general trend with the barotropic trend
-      ! -----------------------------------------------------------------------------
-      !
-      ! Set advection velocity correction:
-      IF (ln_bt_fw) THEN
-         IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN
-            DO_2D( 1, 1, 1, 1 )
-               zun_save = un_adv(ji,jj)
-               zvn_save = vn_adv(ji,jj)
-               !                          ! apply the previously computed correction 
-               un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - rn_atfp * un_bf(ji,jj) )
-               vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - rn_atfp * vn_bf(ji,jj) )
-               !                          ! Update corrective fluxes for next time step
-               un_bf(ji,jj)  = rn_atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )
-               vn_bf(ji,jj)  = rn_atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )
-               !                          ! Save integrated transport for next computation
-               ub2_b(ji,jj) = zun_save
-               vb2_b(ji,jj) = zvn_save
-            END_2D
-         ELSE
-            un_bf(:,:) = 0._wp            ! corrective fluxes for next time step set to zero
-            vn_bf(:,:) = 0._wp
-            ub2_b(:,:) = un_adv(:,:)      ! Save integrated transport for next computation
-            vb2_b(:,:) = vn_adv(:,:)
-         END IF
-      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
-         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( '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(:,:) )
-      ENDIF
-
-
-      ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)  
-      DO jk = 1, jpkm1
-         puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)*r1_hu(:,:,Kmm) - puu_b(:,:,Kmm) ) * umask(:,:,jk)
-         pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) + vn_adv(:,:)*r1_hv(:,:,Kmm) - pvv_b(:,:,Kmm) ) * vmask(:,:,jk)
-      END DO
-
-      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) 
-            pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) & 
-                       & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk)  
-         END DO
-      END IF 
-
-      
-      CALL iom_put(  "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) )    ! barotropic i-current
-      CALL iom_put(  "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) )    ! barotropic i-current
-      !
-#if defined key_agrif
-      ! Save time integrated fluxes during child grid integration
-      ! (used to update coarse grid transports at next time step)
-      !
-      IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
-         IF( Agrif_NbStepint() == 0 ) THEN
-            ub2_i_b(:,:) = 0._wp
-            vb2_i_b(:,:) = 0._wp
-         END IF
-         !
-         za1 = 1._wp / REAL(Agrif_rhot(), wp)
-         ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
-         vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
-      ENDIF
-#endif      
-      !                                   !* write time-spliting arrays in the restart
-      IF( lrst_oce .AND.ln_bt_fw )   CALL ts_rst( kt, 'WRITE' )
-      !
-      IF( ln_wd_il )   DEALLOCATE( zcpx, zcpy )
-      IF( ln_wd_dl )   DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
-      !
-      CALL iom_put( "baro_u" , puu_b(:,:,Kmm) )  ! Barotropic  U Velocity
-      CALL iom_put( "baro_v" , pvv_b(:,:,Kmm) )  ! Barotropic  V Velocity
-      !
-    CONTAINS
-      subroutine loop_fct1(i0, i1, j0, j1, k0, k1, buf)
-        integer, intent(in) :: i0, i1, j0, j1, k0, k1
-        REAL*8, dimension(:,:,:,:,:,:), optional, intent(out) :: buf
-        DO_2D( 0, 0, 0, 0 )
-            zhdiv = (   zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1)   ) * r1_e1e2t(ji,jj)
-            ssha_e(ji,jj) = (  sshn_e(ji,jj) - rDt_e * ( zssh_frc(ji,jj) + zhdiv )  ) * ssmask(ji,jj)
-        END_2D
-      end subroutine loop_fct1
+        IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
+           hu_e (2:jpim1,2:jpjm1) =    hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
+           hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / (  hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1)  )
+           hv_e (2:jpim1,2:jpjm1) =    hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)
+           hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / (  hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1)  )
+        ENDIF
+        !
+      END SUBROUTINE loop_velocity
+      !
    END SUBROUTINE dyn_spg_ts
 
Index: NEMO/branches/2021/dev_r14447_HPC-07_Irrmann_try_new_pt2pt/src/OCE/LBC/lbc_lnk_persistent.h90
===================================================================
--- NEMO/branches/2021/dev_r14447_HPC-07_Irrmann_try_new_pt2pt/src/OCE/LBC/lbc_lnk_persistent.h90	(revision 15044)
+++ NEMO/branches/2021/dev_r14447_HPC-07_Irrmann_try_new_pt2pt/src/OCE/LBC/lbc_lnk_persistent.h90	(revision 15054)
@@ -143,6 +143,6 @@
       ! Size of region
       ifldmax = 6   ! 6 arrays updated max in a single call in dynspg_ts, hypothesis : 1 halo in time splitting
-      icount(1:2) = ifldmax * (jpi-2)   ! west  - east
-      icount(3:4) = ifldmax * (jpj-2)   ! south - north
+      icount(1:2) = ifldmax * (jpj-2)   ! west  - east
+      icount(3:4) = ifldmax * (jpi-2)   ! south - north
       icount(5:8) = ifldmax             ! diagonals
       !
@@ -185,5 +185,5 @@
       ! ----------------------- !
       !
-      CALL MPI_WAITALL(16, nreq_pers, MPI_STATUSES_IGNORE, ierr)
+      CALL MPI_WAITALL( isnd+ircv, nreq_pers, MPI_STATUSES_IGNORE, ierr)
       DO jn = 1, 8
 #define MPI_FILL    
Index: NEMO/branches/2021/dev_r14447_HPC-07_Irrmann_try_new_pt2pt/src/OCE/LBC/lbclnk.F90
===================================================================
--- NEMO/branches/2021/dev_r14447_HPC-07_Irrmann_try_new_pt2pt/src/OCE/LBC/lbclnk.F90	(revision 15044)
+++ NEMO/branches/2021/dev_r14447_HPC-07_Irrmann_try_new_pt2pt/src/OCE/LBC/lbclnk.F90	(revision 15054)
@@ -66,6 +66,8 @@
    PUBLIC   lbc_lnk_icb        ! iceberg lateral boundary conditions
 
-   REAL(dp), DIMENSION(:), ALLOCATABLE ::   buffsnd_dp, buffrcv_dp   ! MPI send/recv buffers
-   REAL(sp), DIMENSION(:), ALLOCATABLE ::   buffsnd_sp, buffrcv_sp   !
+   REAL(sp), DIMENSION(:), ALLOCATABLE ::   buffsnd_sp      , buffrcv_sp         ! MPI send/recv buffers
+   REAL(sp), DIMENSION(:), ALLOCATABLE ::   buffsnd_async_sp, buffrcv_async_sp   ! MPI send/recv buffers
+   REAL(dp), DIMENSION(:), ALLOCATABLE ::   buffsnd_dp      , buffrcv_dp         ! MPI send/recv buffers
+   REAL(dp), DIMENSION(:), ALLOCATABLE ::   buffsnd_async_dp, buffrcv_async_dp   ! MPI send/recv buffers
    INTEGER,  DIMENSION(8)              ::   nreq_p2p                 ! request id for MPI_Isend in point-2-point communication
 
@@ -153,9 +155,13 @@
 #  include "lbc_lnk_oldpt2pt_generic.h90"
 #  include "lbc_lnk_neicoll_generic.h90"
+#  undef BUFFSND
+#  undef BUFFRCV
+#  define BUFFSND buffsnd_async_sp
+#  define BUFFRCV buffrcv_async_sp
+#  include "lbc_lnk_pt2pt_async.h90"
+#  undef BUFFSND
+#  undef BUFFRCV
 #  include "lbc_lnk_persistent.h90"
-#  include "lbc_lnk_pt2pt_async.h90"
 #  undef MPI_TYPE
-#  undef BUFFSND
-#  undef BUFFRCV
 #undef PRECISION
    !!
@@ -170,9 +176,13 @@
 #  include "lbc_lnk_oldpt2pt_generic.h90"
 #  include "lbc_lnk_neicoll_generic.h90"
+#  undef BUFFSND
+#  undef BUFFRCV
+#  define BUFFSND buffsnd_async_dp
+#  define BUFFRCV buffrcv_async_dp
+#  include "lbc_lnk_pt2pt_async.h90"
+#  undef BUFFSND
+#  undef BUFFRCV
 #  include "lbc_lnk_persistent.h90"
-#  include "lbc_lnk_pt2pt_async.h90"
 #  undef MPI_TYPE
-#  undef BUFFSND
-#  undef BUFFRCV
 #undef PRECISION
 
Index: NEMO/branches/2021/dev_r14447_HPC-07_Irrmann_try_new_pt2pt/src/OCE/LBC/lib_mpp.F90
===================================================================
--- NEMO/branches/2021/dev_r14447_HPC-07_Irrmann_try_new_pt2pt/src/OCE/LBC/lib_mpp.F90	(revision 15044)
+++ NEMO/branches/2021/dev_r14447_HPC-07_Irrmann_try_new_pt2pt/src/OCE/LBC/lib_mpp.F90	(revision 15054)
@@ -1179,7 +1179,7 @@
       !!----------------------------------------------------------------------
 
-      if( wp == dp ) then
+      if     ( wp == dp ) then
          MPI_TYPE = MPI_DOUBLE_PRECISION
-      else if ( wp == sp ) then
+      else if( wp == sp ) then
          MPI_TYPE = MPI_REAL
       else
@@ -1190,6 +1190,6 @@
       ! Size of region
       ifldmax = 6   ! 6 arrays updated max in a single call in dynspg_ts
-      icount(1:2) = ifldmax * (jpi-2)   ! west  - east
-      icount(3:4) = ifldmax * (jpj-2)   ! south - north
+      icount(1:2) = ifldmax * (jpj-2)   ! west  - east
+      icount(3:4) = ifldmax * (jpi-2)   ! south - north
       icount(5:8) = ifldmax             ! diagonals
       !
