Index: NEMO/branches/UKMO/dev_r10448_bdyvol/src/OCE/BDY/bdyini.F90
===================================================================
--- NEMO/branches/UKMO/dev_r10448_bdyvol/src/OCE/BDY/bdyini.F90	(revision 10454)
+++ NEMO/branches/UKMO/dev_r10448_bdyvol/src/OCE/BDY/bdyini.F90	(revision 10455)
@@ -133,5 +133,4 @@
       INTEGER  ::   i_offset, j_offset                     !   -       -
       INTEGER , POINTER  ::  nbi, nbj, nbr                 ! short cuts
-      REAL(wp), POINTER  ::  flagu, flagv                  !    -   -
       REAL(wp), POINTER, DIMENSION(:,:)       ::   pmask    ! pointer to 2D mask fields
       REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars
@@ -1244,38 +1243,4 @@
          !
       END DO
-
-      ! Compute total lateral surface for volume correction:
-      ! ----------------------------------------------------
-      ! JC: this must be done at each time step with non-linear free surface
-      bdysurftot = 0._wp 
-      IF( ln_vol ) THEN  
-         igrd = 2      ! Lateral surface at U-points
-         DO ib_bdy = 1, nb_bdy
-            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
-               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
-               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
-               flagu => idx_bdy(ib_bdy)%flagu(ib,igrd)
-               bdysurftot = bdysurftot + hu_n   (nbi  , nbj)                           &
-                  &                    * e2u    (nbi  , nbj) * ABS( flagu ) &
-                  &                    * tmask_i(nbi  , nbj)                           &
-                  &                    * tmask_i(nbi+1, nbj)                   
-            END DO
-         END DO
-
-         igrd=3 ! Add lateral surface at V-points
-         DO ib_bdy = 1, nb_bdy
-            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
-               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
-               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
-               flagv => idx_bdy(ib_bdy)%flagv(ib,igrd)
-               bdysurftot = bdysurftot + hv_n   (nbi, nbj  )                           &
-                  &                    * e1v    (nbi, nbj  ) * ABS( flagv ) &
-                  &                    * tmask_i(nbi, nbj  )                           &
-                  &                    * tmask_i(nbi, nbj+1)
-            END DO
-         END DO
-         !
-         CALL mpp_sum( 'bdyini', bdysurftot )      ! sum over the global domain
-      END IF   
       !
       ! Tidy up
@@ -1284,5 +1249,4 @@
       !
    END SUBROUTINE bdy_segs
-
 
    SUBROUTINE bdy_ctl_seg
Index: NEMO/branches/UKMO/dev_r10448_bdyvol/src/OCE/BDY/bdyvol.F90
===================================================================
--- NEMO/branches/UKMO/dev_r10448_bdyvol/src/OCE/BDY/bdyvol.F90	(revision 10454)
+++ NEMO/branches/UKMO/dev_r10448_bdyvol/src/OCE/BDY/bdyvol.F90	(revision 10455)
@@ -9,4 +9,5 @@
    !!            3.0  !  2008-04  (NEMO team)  add in the reference version
    !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
+   !!            4.0  !  2019-01  (P. Mathiot) adapted to time splitting     
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers 
@@ -24,5 +25,5 @@
    PRIVATE
 
-   PUBLIC   bdy_vol    ! called by ???
+   PUBLIC   bdy_vol2d    ! called by dynspg_ts
 
    !!----------------------------------------------------------------------
@@ -33,5 +34,5 @@
 CONTAINS
 
-   SUBROUTINE bdy_vol( kt )
+   SUBROUTINE bdy_vol2d( kt, kc, pua2d, pva2d, phu, phv )
       !!----------------------------------------------------------------------
       !!                      ***  ROUTINE bdyvol  ***
@@ -40,6 +41,4 @@
       !!      A correction velocity is calculated to correct the total transport 
       !!      through the unstructured OBC. 
-      !!      The total depth used is constant (H0) to be consistent with the 
-      !!      linear free surface coded in OPA 8.2    <<<=== !!gm  ???? true ????
       !!
       !! ** Method  :   The correction velocity (zubtpecor here) is defined calculating
@@ -65,28 +64,28 @@
       !!            (set nn_volctl to 1 in tne namelist for this option)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   kt, kc   ! ocean time-step index, cycle time step
       !
       INTEGER  ::   ji, jj, jk, jb, jgrd
       INTEGER  ::   ib_bdy, ii, ij
-      REAL(wp) ::   zubtpecor, z_cflxemp, ztranst
+      REAL(wp) ::   zubtpecor, ztranst, zeh
+      REAL(wp), SAVE :: z_cflxemp   
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: phu, phv
       TYPE(OBC_INDEX), POINTER :: idx
       !!-----------------------------------------------------------------------------
       !
-      IF( ln_vol ) THEN
-      !
-      IF( kt == nit000 ) THEN 
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*)'bdy_vol : Correction of velocities along unstructured OBC'
-         IF(lwp) WRITE(numout,*)'~~~~~~~'
-      END IF 
-
       ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain
       ! -----------------------------------------------------------------------
-!!gm replace these lines :
-      z_cflxemp = SUM ( ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / rau0
-      CALL mpp_sum( 'bdyvol', z_cflxemp )     ! sum over the global domain
-!!gm   by :
-!!gm      z_cflxemp = glob_sum(  'bdyvol', ( emp(:,:)-rnf(:,:)+fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rau0
-!!gm
+      IF ( kc == 1 ) z_cflxemp = glob_sum( 'bdyvol', ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rau0
+
+      ! Compute bdy surface each cycle if non linear free surface
+      ! ---------------------------------------------------------
+      IF ( .NOT. ln_linssh ) THEN
+         ! compute area each time step
+         bdysurftot = bdy_segs_surf( phu, phv )
+      ELSE
+         ! compute area only the first time
+         IF ( ( kt == nit000 ) .AND. ( kc == 1 ) ) bdysurftot = bdy_segs_surf( phu, phv )
+      END IF
 
       ! Transport through the unstructured open boundary
@@ -98,59 +97,70 @@
          jgrd = 2                               ! cumulate u component contribution first 
          DO jb = 1, idx%nblenrim(jgrd)
-            DO jk = 1, jpkm1
+            ii = idx%nbi(jb,jgrd)
+            ij = idx%nbj(jb,jgrd)
+            zubtpecor = zubtpecor + idx%flagu(jb,jgrd) * pua2d(ii,ij) * e2u(ii,ij) * phu(ii,ij) * tmask_i(ii,ij) * tmask_i(ii+1,ij)
+         END DO
+         jgrd = 3                               ! then add v component contribution
+         DO jb = 1, idx%nblenrim(jgrd)
+            ii = idx%nbi(jb,jgrd)
+            ij = idx%nbj(jb,jgrd)
+            zubtpecor = zubtpecor + idx%flagv(jb,jgrd) * pva2d(ii,ij) * e1v(ii,ij) * phv(ii,ij) * tmask_i(ii,ij) * tmask_i(ii,ij+1)
+         END DO
+         !
+      END DO
+      IF( lk_mpp )   CALL mpp_sum( 'bdyvol', zubtpecor )   ! sum over the global domain
+
+      ! The normal velocity correction
+      ! ------------------------------
+      IF( nn_volctl==1 ) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp ) / bdysurftot  ! maybe should be apply only once at the end
+      ELSE                      ;   zubtpecor =   zubtpecor               / bdysurftot
+      END IF
+
+      ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation
+      ! -------------------------------------------------------------
+      DO ib_bdy = 1, nb_bdy
+         idx => idx_bdy(ib_bdy)
+         !
+         jgrd = 2                               ! correct u component
+         DO jb = 1, idx%nblenrim(jgrd)
                ii = idx%nbi(jb,jgrd)
                ij = idx%nbj(jb,jgrd)
-               zubtpecor = zubtpecor + idx%flagu(jb,jgrd) * ua(ii,ij, jk) * e2u(ii,ij) * e3u_n(ii,ij,jk)
-            END DO
-         END DO
-         jgrd = 3                               ! then add v component contribution
-         DO jb = 1, idx%nblenrim(jgrd)
-            DO jk = 1, jpkm1
+               pua2d(ii,ij) = pua2d(ii,ij) - idx%flagu(jb,jgrd) * zubtpecor * tmask_i(ii,ij) * tmask_i(ii+1,ij)
+         END DO
+         jgrd = 3                              ! correct v component
+         DO jb = 1, idx%nblenrim(jgrd)
                ii = idx%nbi(jb,jgrd)
                ij = idx%nbj(jb,jgrd)
-               zubtpecor = zubtpecor + idx%flagv(jb,jgrd) * va(ii,ij, jk) * e1v(ii,ij) * e3v_n(ii,ij,jk) 
-            END DO
-         END DO
-         !
-      END DO
-      CALL mpp_sum( 'bdyvol', zubtpecor )   ! sum over the global domain
-
-      ! The normal velocity correction
-      ! ------------------------------
-      IF( nn_volctl==1 ) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp ) / bdysurftot 
-      ELSE                      ;   zubtpecor =   zubtpecor               / bdysurftot
-      END IF
-
-      ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation
-      ! -------------------------------------------------------------
-      ztranst = 0._wp
-      DO ib_bdy = 1, nb_bdy
-         idx => idx_bdy(ib_bdy)
-         !
-         jgrd = 2                               ! correct u component
-         DO jb = 1, idx%nblenrim(jgrd)
-            DO jk = 1, jpkm1
-               ii = idx%nbi(jb,jgrd)
-               ij = idx%nbj(jb,jgrd)
-               ua(ii,ij,jk) = ua(ii,ij,jk) - idx%flagu(jb,jgrd) * zubtpecor * umask(ii,ij,jk)
-               ztranst = ztranst + idx%flagu(jb,jgrd) * ua(ii,ij,jk) * e2u(ii,ij) * e3u_n(ii,ij,jk)
-            END DO
-         END DO
-         jgrd = 3                              ! correct v component
-         DO jb = 1, idx%nblenrim(jgrd)
-            DO jk = 1, jpkm1
-               ii = idx%nbi(jb,jgrd)
-               ij = idx%nbj(jb,jgrd)
-               va(ii,ij,jk) = va(ii,ij,jk) -idx%flagv(jb,jgrd) * zubtpecor * vmask(ii,ij,jk)
-               ztranst = ztranst + idx%flagv(jb,jgrd) * va(ii,ij,jk) * e1v(ii,ij) * e3v_n(ii,ij,jk)
-            END DO
-         END DO
-         !
-      END DO
-      CALL mpp_sum( 'bdyvol', ztranst )   ! sum over the global domain
- 
+               pva2d(ii,ij) = pva2d(ii,ij) - idx%flagv(jb,jgrd) * zubtpecor * tmask_i(ii,ij) * tmask_i(ii,ij+1)
+         END DO
+         !
+      END DO
+      ! 
       ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected
       ! ------------------------------------------------------
-      IF( lwp .AND. MOD( kt, nwrite ) == 0 ) THEN
+      IF( MOD( kt, nwrite ) == 0 .AND. ( kc == 1 ) ) THEN
+         !
+         ! compute residual transport across boundary
+         ztranst = 0._wp
+         DO ib_bdy = 1, nb_bdy
+            idx => idx_bdy(ib_bdy)
+            !
+            jgrd = 2                               ! correct u component
+            DO jb = 1, idx%nblenrim(jgrd)
+                  ii = idx%nbi(jb,jgrd)
+                  ij = idx%nbj(jb,jgrd)
+                  ztranst = ztranst + idx%flagu(jb,jgrd) * pua2d(ii,ij) * e2u(ii,ij) * phu(ii,ij) * tmask_i(ii,ij) * tmask_i(ii+1,ij)
+            END DO
+            jgrd = 3                              ! correct v component
+            DO jb = 1, idx%nblenrim(jgrd)
+                  ii = idx%nbi(jb,jgrd)
+                  ij = idx%nbj(jb,jgrd)
+                  ztranst = ztranst + idx%flagv(jb,jgrd) * pva2d(ii,ij) * e1v(ii,ij) * phv(ii,ij) * tmask_i(ii,ij) * tmask_i(ii,ij+1)
+            END DO
+            !
+         END DO
+         IF( lk_mpp )   CALL mpp_sum('bdyvol', ztranst )   ! sum over the global domain
+
+
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*)'bdy_vol : time step :', kt
@@ -162,8 +172,50 @@
       END IF 
       !
-      END IF ! ln_vol
-      !
-   END SUBROUTINE bdy_vol
-
+   END SUBROUTINE bdy_vol2d
+   !
+   REAL(wp) FUNCTION bdy_segs_surf(phu, phv)
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE bdy_ctl_seg  ***
+      !!
+      !! ** Purpose :   Compute total lateral surface for volume correction
+      !!
+      !!----------------------------------------------------------------------
+
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: phu, phv ! water column thickness at U and V points
+      INTEGER            ::  igrd, ib_bdy, ib                ! loop indexes
+      INTEGER , POINTER  ::  nbi, nbj                        ! short cuts
+      REAL(wp), POINTER  ::  zflagu, zflagv                  !    -   -
+
+      ! Compute total lateral surface for volume correction:
+      ! ----------------------------------------------------
+      bdy_segs_surf = 0._wp
+      igrd = 2      ! Lateral surface at U-points
+      DO ib_bdy = 1, nb_bdy
+         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
+            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
+            zflagu => idx_bdy(ib_bdy)%flagu(ib,igrd)
+            bdy_segs_surf = bdy_segs_surf + phu(nbi, nbj)                              &
+               &                            * e2u(nbi, nbj) * ABS( zflagu )            &
+               &                            * tmask_i(nbi, nbj) * tmask_i(nbi+1, nbj)
+         END DO
+      END DO
+
+      igrd=3 ! Add lateral surface at V-points
+      DO ib_bdy = 1, nb_bdy
+         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
+            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
+            zflagv => idx_bdy(ib_bdy)%flagv(ib,igrd)
+            bdy_segs_surf = bdy_segs_surf + phv(nbi, nbj)                              &
+               &                            * e1v(nbi, nbj) * ABS( zflagv )            &
+               &                            * tmask_i(nbi, nbj) * tmask_i(nbi, nbj+1)
+         END DO
+      END DO
+      !
+      ! redirect the time to bdyvol as this variable is only used by bdyvol
+      IF( lk_mpp )   CALL mpp_sum( 'bdyvol', bdy_segs_surf ) ! sum over the global domain
+      !
+   END FUNCTION bdy_segs_surf
    !!======================================================================
 END MODULE bdyvol
Index: NEMO/branches/UKMO/dev_r10448_bdyvol/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/UKMO/dev_r10448_bdyvol/src/OCE/DYN/dynspg_ts.F90	(revision 10454)
+++ NEMO/branches/UKMO/dev_r10448_bdyvol/src/OCE/DYN/dynspg_ts.F90	(revision 10455)
@@ -41,4 +41,5 @@
    USE wet_dry         ! wetting/drying flux limter
    USE bdy_oce         ! open boundary
+   USE bdyvol          ! open boundary volume conservation
    USE bdytides        ! open boundary condition data
    USE bdydyn2d        ! open boundary conditions on barotropic variables
@@ -792,6 +793,7 @@
          !                                                !* after ssh
          !                                                !  -----------
-         ! One should enforce volume conservation at open boundaries here
-         ! considering fluxes below:
+         !
+         ! Enforce volume conservation at open boundaries:
+         IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e )
          !
          zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5
