Changeset 3970 for branches/2013/dev_r3867_MERCATOR1_DYN

Timestamp:

2013-07-11T15:59:14+02:00 (11 years ago)

Author:

cbricaud

Message:

Time splitting update, see ticket #1079

Location:

branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC

Files:

• BDY/bdydta.F90 (modified) (6 diffs)
• BDY/bdydyn.F90 (modified) (2 diffs)
• BDY/bdydyn2d.F90 (modified) (7 diffs)
• BDY/bdyini.F90 (modified) (15 diffs)
• BDY/bdytides.F90 (modified) (9 diffs)
• DOM/dom_oce.F90 (modified) (3 diffs)
• DOM/domain.F90 (modified) (3 diffs)
• DOM/domvvl.F90 (modified) (4 diffs)
• DYN/dynadv_ubs.F90 (modified) (3 diffs)
• DYN/dynnxt.F90 (modified) (10 diffs)
• DYN/dynspg.F90 (modified) (4 diffs)
• DYN/dynspg_ts.F90 (modified) (15 diffs)
• DYN/dynzdf_imp.F90 (modified) (19 diffs)
• DYN/sshwzv.F90 (modified) (10 diffs)
• SBC/updtide.F90 (modified) (5 diffs)
• nemogcm.F90 (modified) (1 diff)
• oce.F90 (modified) (2 diffs)
• step.F90 (modified) (4 diffs)

branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

@@ -29 +29 @@
    USE iom             ! IOM library
    USE in_out_manager  ! I/O logical units
+   ! bg jchanut tschanges
+   USE dynspg_oce, ONLY: lk_dynspg_ts ! Split-explicit free surface flag
+   ! end jchanut tschanges
+
 #if defined key_lim2
    USE ice_2
@@ -314 +318 @@
                      END DO
                   ENDIF
-                  IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing
-                     CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy),  &
-                                        & td=tides(ib_bdy), time_offset=time_offset )
-                  ENDIF
+                  ! bg jchanut tschanges
+                  !IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing
+                  !   CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy),  &
+                  !                      & td=tides(ib_bdy), time_offset=time_offset )
+                  !ENDIF
+                  ! end jchanut tschanges
                ENDIF
             ENDIF
@@ -323 +329 @@
          END IF ! nn_dta(ib_bdy) = 1
       END DO  ! ib_bdy
+
+      ! bg jchanut tschanges
+#if defined key_tide
+      ! Add tides if not split-explicit free surface else this is done in ts loop
+      IF (.NOT.lk_dynspg_ts) CALL bdy_dta_tides( kt=kt, time_offset=time_offset )
+#endif
+      ! end jchanut tschanges
 
       IF ( ln_apr_obc ) THEN
@@ -476 +489 @@
             IF( nn_dyn2d(ib_bdy) .gt. 0 .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN 
 
-               IF( nn_tra(ib_bdy) .ne. 4 ) THEN ! runoff condition : no ssh reading
+               IF( nn_dyn2d(ib_bdy) .ne. jp_frs .and. nn_tra(ib_bdy) .ne. 4 ) THEN ! runoff condition : no ssh reading
                   jfld = jfld + 1
                   blf_i(jfld) = bn_ssh
@@ -572 +585 @@
             ! Recalculate field counts
             !-------------------------
-            nb_bdy_fld_sum = 0
             IF( ib_bdy .eq. 1 ) THEN 
+               nb_bdy_fld_sum = 0
                nb_bdy_fld(ib_bdy) = jfld
                nb_bdy_fld_sum     = jfld              
@@ -616 +629 @@
                ALLOCATE( dta_bdy(ib_bdy)%u2d(ilen0(2)) )
                ALLOCATE( dta_bdy(ib_bdy)%v2d(ilen0(3)) )
-               IF (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3) THEN
+               IF ( nn_dyn2d(ib_bdy) .ne. jp_frs .and. (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3) ) THEN
                   jfld = jfld + 1
                   dta_bdy(ib_bdy)%ssh => bf(jfld)%fnow(:,1,1)

branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90

@@ -30 +30 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  !
+   USE domvvl          ! variable volume
 
    IMPLICIT NONE
@@ -84 +85 @@
       pu2d(:,:) = 0.e0
       pv2d(:,:) = 0.e0
-      DO jk = 1, jpkm1   !! Vertically integrated momentum trends
-          pu2d(:,:) = pu2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
-          pv2d(:,:) = pv2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
-      END DO
-      pu2d(:,:) = pu2d(:,:) * phur(:,:)
-      pv2d(:,:) = pv2d(:,:) * phvr(:,:)
+      ! bg jchanut tschanges (not specifically related to ts; this is a bug)
+      IF (lk_vvl) THEN
+         DO jk = 1, jpkm1   !! Vertically integrated momentum trends
+            pu2d(:,:) = pu2d(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
+            pv2d(:,:) = pv2d(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
+         END DO
+         pu2d(:,:) = pu2d(:,:) / ( hu_0(:,:) + sshu_a(:,:) + 1._wp - umask(:,:,1) )
+         pv2d(:,:) = pv2d(:,:) / ( hv_0(:,:) + sshv_a(:,:) + 1._wp - vmask(:,:,1) )
+      ! end jchanut tschanges
+      ELSE
+         DO jk = 1, jpkm1   !! Vertically integrated momentum trends
+            pu2d(:,:) = pu2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
+            pv2d(:,:) = pv2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
+         END DO
+         pu2d(:,:) = pu2d(:,:) * phur(:,:)
+         pv2d(:,:) = pv2d(:,:) * phvr(:,:)
+      ENDIF
+
       DO jk = 1 , jpkm1
-         ua(:,:,jk) = ua(:,:,jk) - pu2d(:,:)
-         va(:,:,jk) = va(:,:,jk) - pv2d(:,:)
+         ua(:,:,jk) = ua(:,:,jk) - pu2d(:,:) * umask(:,:,jk)
+         va(:,:,jk) = va(:,:,jk) - pv2d(:,:) * vmask(:,:,jk)
       END DO
 

branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90

@@ -6 +6 @@
    !! History :  3.4  !  2011     (D. Storkey) new module as part of BDY rewrite
    !!            3.5  !  2012     (S. Mocavero, I. Epicoco) Optimization of BDY communications
+   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
    !!----------------------------------------------------------------------
 #if defined key_bdy 
@@ -27 +28 @@
 
    PUBLIC   bdy_dyn2d     ! routine called in dynspg_ts and bdy_dyn
+   PUBLIC   bdy_ssh       ! routine called in dynspg_ts or sshwzv
 
    !!----------------------------------------------------------------------
@@ -135 +137 @@
       REAL(wp) ::   zcorr                            ! Flather correction
       REAL(wp) ::   zforc                            ! temporary scalar
+      REAL(wp) ::   zflag, z1_2                      !    "        "
       !!----------------------------------------------------------------------
 
       IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_fla')
+
+      z1_2 = 0.5_wp
 
       ! ---------------------------------!
@@ -164 +169 @@
          !
          zcorr = - idx%flagu(jb) * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
-         zforc = dta%u2d(jb)
-         pu2d(ii,ij) = zforc + zcorr * umask(ii,ij,1) 
+         ! bg jchanut tschanges: Set zflag to 0 below to revert to Flather scheme
+!!         zforc = dta%u2d(jb)
+         zflag = ABS(idx%flagu(jb))
+         iim1 = ii + idx%flagu(jb)
+         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pu2d(iim1,ij)
+         pu2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1) 
+         ! end jchanut tschanges
       END DO
       !
@@ -177 +187 @@
          !
          zcorr = - idx%flagv(jb) * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
-         zforc = dta%v2d(jb)
-         pv2d(ii,ij) = zforc + zcorr * vmask(ii,ij,1)
+         ! bg jchanut tschanges: Set zflag to 0 below to revert to std Flather scheme
+!!         zforc = dta%v2d(jb)
+         zflag = ABS(idx%flagv(jb))
+         ijm1 = ij + idx%flagv(jb)
+         zforc = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pv2d(ii,ijm1)
+         pv2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1)
+         ! end jchanut tschanges
       END DO
       CALL lbc_bdy_lnk( pu2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
@@ -186 +201 @@
       !
    END SUBROUTINE bdy_dyn2d_fla
+
+   SUBROUTINE bdy_ssh( zssh )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE bdy_ssh  ***
+      !!
+      !! ** Purpose : Duplicate sea level across open boundaries
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zssh ! Sea level
+      !!
+      INTEGER  ::   ib_bdy, ib, igrd                        ! local integers
+      INTEGER  ::   ii, ij, zcoef, zcoef1, zcoef2, ip, jp   !   "       "
+
+      igrd = 1                       ! Everything is at T-points here
+
+      DO ib_bdy = 1, nb_bdy
+         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+            ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+            ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+            ! Set gradient direction:
+            zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  )
+            zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1)
+            IF ( zcoef1+zcoef2 == 0 ) THEN
+               ! corner
+!               zcoef = tmask(ii-1,ij,1) + tmask(ii+1,ij,1) +  tmask(ii,ij-1,1) +  tmask(ii,ij+1,1)
+!               zssh(ii,ij) = zssh(ii-1,ij  ) * tmask(ii-1,ij  ,1) + &
+!                 &           zssh(ii+1,ij  ) * tmask(ii+1,ij  ,1) + &
+!                 &           zssh(ii  ,ij-1) * tmask(ii  ,ij-1,1) + &
+!                 &           zssh(ii  ,ij+1) * tmask(ii  ,ij+1,1)
+               zcoef = bdytmask(ii-1,ij) + bdytmask(ii+1,ij) +  bdytmask(ii,ij-1) +  bdytmask(ii,ij+1)
+               zssh(ii,ij) = zssh(ii-1,ij  ) * bdytmask(ii-1,ij  ) + &
+                 &           zssh(ii+1,ij  ) * bdytmask(ii+1,ij  ) + &
+                 &           zssh(ii  ,ij-1) * bdytmask(ii  ,ij-1) + &
+                 &           zssh(ii  ,ij+1) * bdytmask(ii  ,ij+1)
+               zssh(ii,ij) = ( zssh(ii,ij) / MAX( 1, zcoef) ) * tmask(ii,ij,1)
+            ELSE
+               ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )
+               jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)
+               zssh(ii,ij) = zssh(ii+ip,ij+jp) * tmask(ii+ip,ij+jp,1)
+            ENDIF
+         END DO
+
+         ! Boundary points should be updated
+         CALL lbc_bdy_lnk( zssh(:,:), 'T', 1., ib_bdy )
+      END DO
+
+   END SUBROUTINE bdy_ssh
 #else
    !!----------------------------------------------------------------------
@@ -192 +254 @@
 CONTAINS
    SUBROUTINE bdy_dyn2d( kt )      ! Empty routine
-      WRITE(*,*) 'bdy_dyn_frs: You should not have seen this print! error?', kt
+      INTEGER, intent(in) :: kt
+      WRITE(*,*) 'bdy_dyn2: You should not have seen this print! error?', kt
    END SUBROUTINE bdy_dyn2d
+
 #endif
 
    !!======================================================================
 END MODULE bdydyn2d
+

branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

@@ -161 +161 @@
 
       DO ib_bdy = 1,nb_bdy
+        icount = 0 ! flag to set max rimwidth to 1 if no relaxation
         IF(lwp) WRITE(numout,*) ' ' 
         IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------' 
@@ -175 +176 @@
           CASE(jp_none)         ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
           CASE(jp_frs)          ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+          icount = icount + 1
           CASE(jp_flather)      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition'
           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' )
@@ -196 +198 @@
           CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
           CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+          icount = icount + 1
           CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value'
           CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)'
@@ -215 +218 @@
               CALL ctl_stop( 'Use FRS OR relaxation' )
            ELSE
+              icount = icount + 1
               IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone'
               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
@@ -228 +232 @@
           CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
           CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+          icount = icount + 1
           CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value'
           CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Neumann conditions'
@@ -248 +253 @@
               CALL ctl_stop( 'Use FRS OR relaxation' )
            ELSE
+              icount = icount + 1
               IF(lwp) WRITE(numout,*) '      + T/S relaxation zone'
               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
@@ -262 +268 @@
           CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
           CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+          icount = icount + 1
           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
         END SELECT
@@ -273 +280 @@
         IF(lwp) WRITE(numout,*)
 #endif
-
-        IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy)
-        IF(lwp) WRITE(numout,*)
+        IF ( icount>0 ) THEN
+           IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy)
+           IF(lwp) WRITE(numout,*)
+        ELSE
+           nn_rimwidth(ib_bdy) = 1 ! no relaxation
+        ENDIF
 
       ENDDO
@@ -391 +401 @@
             ENDDO
             CALL iom_close( inum )
-
          ENDIF 
 
@@ -398 +407 @@
       IF (nb_bdy>0) THEN
          jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy))
-
          ! Allocate arrays
          !---------------
@@ -446 +454 @@
          ENDIF 
 
-      ENDDO      
+      ENDDO     
     
       ! 2. Now fill indices corresponding to straight open boundary arrays:
@@ -752 +760 @@
                ! check if point is in local domain
                IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   &
-                  & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in       ) THEN
+                  & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in .AND.   &
+                  & nbrdta(ib,igrd,ib_bdy) <= nn_rimwidth(ib_bdy)     ) THEN      
                   !
                   icount = icount  + 1
@@ -765 +774 @@
          ! Allocate index arrays for this boundary set
          !--------------------------------------------
-         ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:))
+
+         ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(1:jpbgrd))
+         ilen1 = MAX(1,ilen1)
          ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) )
          ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) )
@@ -773 +784 @@
          ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) )
          ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1) )
-         ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) )
+         ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) )      
 
          ! Dispatch mapping indices and discrete distances on each processor
@@ -1008 +1019 @@
       ! bdytmask = 1  on the computational domain AND on open boundaries
       !          = 0  elsewhere   
- 
+      bdytmask(:,:) = 1.e0
+      bdyumask(:,:) = 1.e0
+      bdyvmask(:,:) = 1.e0
+
       IF( ln_mask_file ) THEN
          CALL iom_open( cn_mask_file, inum )

branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90

@@ -9 +9 @@
    !!            3.3  !  2010-09  (D.Storkey and E.O'Dea)  bug fixes
    !!            3.4  !  2012-09  (G. Reffray and J. Chanut) New inputs + mods
+   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
    !!----------------------------------------------------------------------
 #if defined key_bdy
@@ -32 +33 @@
 !   USE tide_mod       ! Useless ??
    USE fldread, ONLY: fld_map
+   USE dynspg_oce, ONLY: lk_dynspg_ts
 
    IMPLICIT NONE
@@ -38 +40 @@
    PUBLIC   bdytide_init     ! routine called in bdy_init
    PUBLIC   bdytide_update   ! routine called in bdy_dta
+   PUBLIC   bdy_dta_tides    ! routine called in dyn_spg_ts
 
    TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data
@@ -49 +52 @@
 
    TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data
+   TYPE(OBC_DATA)  , PRIVATE, DIMENSION(jp_bdy) :: dta_bdy_s  !: bdy external data (slow component)
 
    !!----------------------------------------------------------------------
@@ -252 +256 @@
             ENDIF
             !
+            IF ( lk_dynspg_ts ) THEN ! Allocate arrays to save slowly varying boundary data during
+                                     ! time splitting integration
+               ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) )
+               ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) )
+               ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) )
+               dta_bdy_s(ib_bdy)%ssh(:) = 0.e0
+               dta_bdy_s(ib_bdy)%u2d(:) = 0.e0
+               dta_bdy_s(ib_bdy)%v2d(:) = 0.e0
+            ENDIF
+            !
          ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2
          !
@@ -318 +332 @@
          
       IF( PRESENT(jit) ) THEN  
-         z_arg = ( ((kt-kt_tide)-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) )
+         z_arg = ((kt-kt_tide) * rdt + (jit+0.5*(time_add-1)) * rdt / REAL(nn_baro,wp) )
       ELSE                              
          z_arg = ((kt-kt_tide)+time_add) * rdt
@@ -324 +338 @@
 
       ! Linear ramp on tidal component at open boundaries 
-      zramp = 1.
-      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0.),1.)
+      zramp = 1._wp
+      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0._wp),1._wp)
 
       DO itide = 1, nb_harmo
@@ -351 +365 @@
       !
    END SUBROUTINE bdytide_update
+
+   SUBROUTINE bdy_dta_tides( kt, kit, time_offset )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_dta_tides  ***
+      !!                
+      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays. 
+      !!                
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in )            ::   kt          ! Main timestep counter
+      INTEGER, INTENT( in ),OPTIONAL   ::   kit         ! Barotropic timestep counter (for timesplitting option)
+      INTEGER, INTENT( in ),OPTIONAL   ::   time_offset ! time offset in units of timesteps. NB. if kit
+                                                        ! is present then units = subcycle timesteps.
+                                                        ! time_offset = 0  => get data at "now"    time level
+                                                        ! time_offset = -1 => get data at "before" time level
+                                                        ! time_offset = +1 => get data at "after"  time level
+                                                        ! etc.
+      !!
+      LOGICAL  :: lk_first_btstp  ! =.TRUE. if time splitting and first barotropic step
+      INTEGER,          DIMENSION(jpbgrd) :: ilen0 
+      INTEGER, DIMENSION(1:jpbgrd) :: nblen, nblenrim  ! short cuts
+      INTEGER  :: itide, ib_bdy, ib, igrd                     ! loop indices
+      INTEGER  :: time_add                                    ! time offset in units of timesteps
+      REAL(wp) :: z_arg, z_sarg, zramp, zoff, z_cost, z_sist      
+      !!----------------------------------------------------------------------
+
+      IF( nn_timing == 1 ) CALL timing_start('bdy_dta_tides')
+
+      lk_first_btstp=.TRUE.
+      IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF
+
+      time_add = 0
+      IF( PRESENT(time_offset) ) THEN
+         time_add = time_offset
+      ENDIF
+      
+      ! Absolute time from model initialization:   
+      IF( PRESENT(kit) ) THEN  
+         z_arg = ( kt + (kit+0.5*(time_add-1)) / REAL(nn_baro,wp) ) * rdt
+      ELSE                              
+         z_arg = ( kt + time_add ) * rdt
+      ENDIF
+
+      ! Linear ramp on tidal component at open boundaries 
+      zramp = 1.
+      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rdttideramp*rday),0.),1.)
+
+      DO ib_bdy = 1,nb_bdy
+
+         ! line below should be simplified (runoff case)
+         IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(nn_tra(ib_bdy).NE.4)) THEN
+
+            nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd)
+            nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd)
+
+            IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN
+               ilen0(:)=nblen(:)
+            ELSE
+               ilen0(:)=nblenrim(:)
+            ENDIF     
+
+            ! We refresh nodal factors every day below
+            ! This should be done somewhere else
+            IF ( nsec_day == NINT(0.5 * rdttra(1)) .AND. lk_first_btstp ) THEN
+               !
+               kt_tide = kt               
+               !
+               IF(lwp) THEN
+               WRITE(numout,*)
+               WRITE(numout,*) 'bdy_tide_dta : Refresh nodal factors for tidal open bdy data at kt=',kt
+               WRITE(numout,*) '~~~~~~~~~~~~~~ '
+               ENDIF
+               !
+               CALL tide_init_elevation ( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
+               CALL tide_init_velocities( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
+               !
+            ENDIF
+            zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time
+            !
+            ! If time splitting, save data at first barotropic iteration
+            IF ( PRESENT(kit) ) THEN
+               IF ( lk_first_btstp ) THEN ! Save slow varying open boundary data:
+                  dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy(ib_bdy)%ssh(1:ilen0(1))
+                  dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy(ib_bdy)%u2d(1:ilen0(2))
+                  dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy(ib_bdy)%v2d(1:ilen0(3))
+
+               ELSE ! Initialize arrays from slow varying open boundary data:            
+                  dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1))
+                  dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2))
+                  dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3))
+               ENDIF
+            ENDIF
+            !
+            ! Update open boundary data arrays:
+            DO itide = 1, nb_harmo
+               !
+               z_sarg = (z_arg + zoff) * omega_tide(itide)
+               z_cost = zramp * COS( z_sarg )
+               z_sist = zramp * SIN( z_sarg )
+               !
+               igrd=1                              ! SSH on tracer grid
+               DO ib = 1, ilen0(igrd)
+                  dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + &
+                     &                      ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + &
+                     &                        tides(ib_bdy)%ssh(ib,itide,2)*z_sist )
+               END DO
+               !
+               igrd=2                              ! U grid
+               DO ib = 1, ilen0(igrd)
+                  dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + &
+                     &                      ( tides(ib_bdy)%u(ib,itide,1)*z_cost + &
+                     &                        tides(ib_bdy)%u(ib,itide,2)*z_sist )
+               END DO
+               !
+               igrd=3                              ! V grid
+               DO ib = 1, ilen0(igrd) 
+                  dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + &
+                     &                      ( tides(ib_bdy)%v(ib,itide,1)*z_cost + &
+                     &                        tides(ib_bdy)%v(ib,itide,2)*z_sist )
+               END DO
+            END DO
+         END IF
+      END DO
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_dta_tides')
+      !
+   END SUBROUTINE bdy_dta_tides
 
    SUBROUTINE tide_init_elevation( idx, td )
@@ -457 +597 @@
       WRITE(*,*) 'bdytide_update: You should not have seen this print! error?', kt, jit
    END SUBROUTINE bdytide_update
+   SUBROUTINE bdy_dta_tides( kt, jit )   ! Empty routine
+      WRITE(*,*) 'bdy_dta_tides: You should not have seen this print! error?', kt, jit
+   END SUBROUTINE bdy_dta_tides
 #endif
 
    !!======================================================================
 END MODULE bdytides
+

branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

@@ -158 +158 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   hur  , hvr    !: inverse of u and v-points ocean depth (1/m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   hu   , hv     !: depth at u- and v-points (meters)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   hu_0 , hv_0   !: refernce depth at u- and v-points (meters)
-
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   hu_0 , hv_0   !: reference depth at u- and v-points (meters)
+   ! bg jchanut tschanges 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ht   , hf     !: depth at t- and f-points (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ht_0 , hf_0   !: reference depth at t- and f-points (meters)
+   ! end jchanut tschanges 
    INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1)
    INTEGER, PUBLIC ::   nlb10              !: shallowest W level Bellow ~10m (nla10 + 1) 
@@ -264 +267 @@
    INTEGER FUNCTION dom_oce_alloc()
       !!----------------------------------------------------------------------
-      INTEGER, DIMENSION(11) :: ierr
+      ! bg jchanut tschanges
+      INTEGER, DIMENSION(13) :: ierr
+      ! end jchanut tschanges
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -314 +319 @@
       ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(11) )
 #endif
+      ! bg jchanut tschanges
+#if defined key_dynspg_ts
+      ALLOCATE( ht(jpi,jpj)  , hf(jpi,jpj)  , STAT=ierr(12) )
+#if defined key_vvl
+      ALLOCATE( ht_0(jpi,jpj), hf_0(jpi,jpj), STAT=ierr(13) )
+#endif
+#endif
+      ! end jchanut tschanges
       !
       dom_oce_alloc = MAXVAL(ierr)

branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

@@ -68 +68 @@
       !!              - 1D configuration, move Coriolis, u and v at T-point
       !!----------------------------------------------------------------------
-      INTEGER ::   jk          ! dummy loop argument
+      INTEGER ::   jj, jk      ! dummy loop arguments
       INTEGER ::   iconf = 0   ! local integers
       !!----------------------------------------------------------------------
@@ -100 +100 @@
       hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1)
 
+      ! bg jchanut tschanges
+#if defined key_dynspg_ts
+      !
+      ht(:,:) = 0._wp                         ! Ocean depth at T-point
+      DO jk = 1, jpk
+         ht(:,:) = ht(:,:) + fse3t(:,:,jk) * tmask(:,:,jk)
+      END DO 
+
+                                              ! Ocean depth at F-point 
+                                              ! Ensure that depth is non-zero over land
+      IF ( .not. ln_sco ) THEN
+         IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
+         ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
+         ENDIF
+         hf(:,:) = gdepw_0(jk+1)
+      ELSE
+         hf(:,:) = hbatf(:,:)
+      END IF
+
+      DO jj = 1, jpjm1
+         hf(:,jj) = hf(:,jj)*(1._wp- umask(:,jj,1) * umask(:,jj+1,1))
+      END DO 
+            
+      DO jk = 1, jpkm1
+         DO jj = 1, jpjm1
+            hf(:,jj) = hf(:,jj) + fse3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
+         END DO
+      END DO
+      CALL lbc_lnk( hf, 'F', 1._wp )  
+#endif 
+      ! end jchanut tschanges
                              CALL dom_stp      ! time step
       IF( nmsh /= 0      )   CALL dom_wri      ! Create a domain file
@@ -218 +249 @@
          WRITE(numout,*) '      ocean time step                      rn_rdt    = ', rn_rdt
          WRITE(numout,*) '      asselin time filter parameter        rn_atfp   = ', rn_atfp
-         WRITE(numout,*) '      time-splitting: nb of sub time-step  nn_baro   = ', nn_baro
+         ! bg jchanut tschanges
+!         WRITE(numout,*) '      time-splitting: nb of sub time-step  nn_baro   = ', nn_baro
+         ! end jchanut tschanges
          WRITE(numout,*) '      acceleration of converge             nn_acc    = ', nn_acc
          WRITE(numout,*) '        nn_acc=1: surface tracer rdt       rn_rdtmin = ', rn_rdtmin

branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

@@ -141 +141 @@
       END DO
       
+#if defined key_dynspg_ts
+      ! bg jchanut tschanges
+      ht_0(:,:) = 0._wp                         ! Reference ocean depth at T-points
+      DO jk = 1, jpk
+         ht_0(:,:) = ht_0(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk)
+      END DO 
+                                                ! Reference ocean depth at F-points 
+                                                ! Ensure that depth is non-zero over land
+      IF ( .not. ln_sco ) THEN
+         IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
+         ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
+         ENDIF
+         hf_0(:,:) = gdepw_0(jk+1)
+      ELSE
+         hf_0(:,:) = hbatf(:,:)
+      END IF
+
+      DO jj = 1, jpjm1
+         hf_0(:,jj) = hf_0(:,jj)*(1._wp- umask(:,jj,1) * umask(:,jj+1,1))
+      END DO  
+           
+      DO jk = 1, jpkm1
+         DO jj = 1, jpjm1
+            hf_0(:,jj) = hf_0(:,jj) + fse3f_0(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
+         END DO
+      END DO
+      CALL lbc_lnk( hf_0, 'F', 1._wp )
+      ! end jchanut tschanges
+#endif  
+
       DO jj = 1, jpjm1                          ! initialise before and now Sea Surface Height at u-, v-, f-points
          DO ji = 1, jpim1   ! NO vector opt.
@@ -192 +222 @@
       INTEGER  ::   iku, ikv     ! local integers    
       INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers
-      REAL(wp) ::   zvt          ! local scalars
+      REAL(wp) ::   zvt, zvtip1, zvtjp1  ! local scalars
       !!----------------------------------------------------------------------
       !
@@ -202 +232 @@
          WRITE(numout,*) '~~~~~~~~~ '
          pe3u_b(:,:,jpk) = fse3u_0(:,:,jpk)
-         pe3v_b(:,:,jpk) = fse3u_0(:,:,jpk)
+         pe3v_b(:,:,jpk) = fse3v_0(:,:,jpk)
       ENDIF
       
@@ -208 +238 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
-               zvt = fse3t_b(ji,jj,jk) * e1e2t(ji,jj)
-               pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1e2t(ji+1,jj) ) / ( e1u(ji,jj) * e2u(ji,jj) )
-               pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e1e2t(ji,jj+1) ) / ( e1v(ji,jj) * e2v(ji,jj) )
+               zvt    = ( fse3t_b(ji  ,jj  ,jk) - fse3t_0(ji  ,jj  ,jk) ) * e1e2t(ji  ,jj  )
+               zvtip1 = ( fse3t_b(ji+1,jj  ,jk) - fse3t_0(ji+1,jj  ,jk) ) * e1e2t(ji+1,jj  )
+               zvtjp1 = ( fse3t_b(ji  ,jj+1,jk) - fse3t_0(ji  ,jj+1,jk) ) * e1e2t(ji  ,jj+1)
+               pe3u_b(ji,jj,jk) = fse3u_0(ji,jj,jk) + 0.5_wp * ( zvt + zvtip1 ) / ( e1u(ji,jj) * e2u(ji,jj) )
+               pe3v_b(ji,jj,jk) = fse3v_0(ji,jj,jk) + 0.5_wp * ( zvt + zvtjp1 ) / ( e1v(ji,jj) * e2v(ji,jj) )
             END DO
          END DO

branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90

@@ -29 +29 @@
 
    REAL(wp), PARAMETER :: gamma1 = 1._wp/3._wp  ! =1/4 quick      ; =1/3  3rd order UBS
-   REAL(wp), PARAMETER :: gamma2 = 1._wp/8._wp  ! =0   2nd order  ; =1/8  4th order centred
+   REAL(wp), PARAMETER :: gamma2 = 1._wp/32._wp ! =0   2nd order  ; =1/32 4th order centred
 
    PUBLIC   dyn_adv_ubs   ! routine called by step.F90
@@ -57 +57 @@
       !!                       = 1/3  3rd order Upstream biased scheme
       !!                gamma2 = 0    2nd order finite differencing 
-      !!                       = 1/8  4th order finite differencing
+      !!                       = 1/32 4th order finite differencing
       !!      For stability reasons, the first term of the fluxes which cor-
       !!      responds to a second order centered scheme is evaluated using  
@@ -64 +64 @@
       !!      before velocity (forward in time). 
       !!      Default value (hard coded in the begining of the module) are 
-      !!      gamma1=1/3 and gamma2=1/8.
+      !!      gamma1=1/3 and gamma2=1/32.
       !!
       !! ** Action : - (ua,va) updated with the 3D advective momentum trends

branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

@@ -17 +17 @@
    !!            3.3  !  2010-09  (D. Storkey, E.O'Dea) Bug fix for BDY module
    !!            3.3  !  2011-03  (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL
+   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
    !!-------------------------------------------------------------------------
   
@@ -42 +43 @@
    USE wrk_nemo        ! Memory Allocation
    USE prtctl          ! Print control
+   USE dynspg_ts       ! Barotropic velocities
+
 #if defined key_agrif
    USE agrif_opa_interp
@@ -103 +106 @@
       REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec   ! local scalars
       REAL(wp) ::   zve3a, zve3n, zve3b, zvf        !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zua, zva
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f 
       !!----------------------------------------------------------------------
@@ -109 +113 @@
       !
       CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
+      IF ( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zua, zva )
       !
       IF( kt == nit000 ) THEN
@@ -127 +132 @@
       !
 #else
+
+#if defined key_dynspg_exp
       ! Next velocity :   Leap-frog time stepping
       ! -------------
@@ -147 +154 @@
          END DO
       ENDIF
-
+#endif
+
+#if defined key_dynspg_ts
+      ! Ensure below that barotropic velocities match time splitting estimate
+      ! Compute actual transport and replace it with ts estimate at "after" time step
+      zua(:,:) = 0._wp
+      zva(:,:) = 0._wp
+      IF (lk_vvl) THEN
+         DO jk = 1, jpkm1
+            zua(:,:) = zua(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
+            zva(:,:) = zva(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)    
+         END DO
+         DO jk = 1, jpkm1
+            ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) &
+            & / ( hu_0(:,:) + sshu_a(:,:) + 1._wp - umask(:,:,1) ) + ua_b(:,:) ) * umask(:,:,jk)
+            va(:,:,jk) = ( va(:,:,jk) - zva(:,:) & 
+            & / ( hv_0(:,:) + sshv_a(:,:) + 1._wp - vmask(:,:,1) ) + va_b(:,:) ) * vmask(:,:,jk)
+         END DO
+      ELSE
+         DO jk = 1, jpkm1
+            zua(:,:) = zua(:,:) + fse3u(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
+            zva(:,:) = zva(:,:) + fse3v(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)    
+         END DO
+         DO jk = 1, jpkm1
+            ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * hur(:,:) + ua_b(:,:) ) *umask(:,:,jk)
+            va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr(:,:) + va_b(:,:) ) *vmask(:,:,jk)
+         END DO
+      ENDIF
+
+      IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN
+         ! Remove advective velocity from "now velocities" 
+         ! prior to asselin filtering     
+         ! In the forward case, this is done below after asselin filtering    
+         DO jk = 1, jpkm1
+            un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
+            vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
+         END DO  
+      ENDIF
+#endif
 
       ! Update after velocity on domain lateral boundaries
@@ -194 +239 @@
             vn(:,:,jk) = va(:,:,jk)
          END DO
+         IF (lk_vvl) THEN
+            DO jk = 1, jpkm1
+               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
+               fse3u_b(:,:,jk) = fse3u_n(:,:,jk)
+               fse3v_b(:,:,jk) = fse3v_n(:,:,jk)
+            ENDDO
+         ENDIF
       ELSE                                             !* Leap-Frog : Asselin filter and swap
          !                                ! =============!
@@ -215 +267 @@
             !                             ! ================!
             !
-            DO jk = 1, jpkm1                 ! Before scale factor at t-points
-               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
-                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
-                  &                         - 2._wp * fse3t_n(:,:,jk)            )
-            END DO
-            zec = atfp * rdt / rau0          ! Add filter correction only at the 1st level of t-point scale factors
-            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
+            IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
+               ! Remove asselin filtering on thicknesses if forward time splitting
+               DO jk = 1, jpkm1   
+                  fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
+               ENDDO
+            ELSE
+
+               DO jk = 1, jpkm1                 ! Before scale factor at t-points
+                  fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
+                     &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
+                     &                         - 2._wp * fse3t_n(:,:,jk)            )
+               END DO
+               zec = atfp * rdt / rau0          ! Add filter correction only at the 1st level of t-point scale factors
+               fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
+            ENDIF
             !
             IF( ln_dynadv_vec ) THEN         ! vector invariant form (no thickness weighted calulation)
@@ -272 +332 @@
          ENDIF
          !
-      ENDIF
+#if defined key_dynspg_ts
+         IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
+         ! Remove asselin filtering of barotropic velocities if forward time splitting
+         ! note that we replace barotropic velocities by advective velocities       
+            zua(:,:) = 0._wp
+            zva(:,:) = 0._wp
+            IF (lk_vvl) THEN
+               DO jk = 1, jpkm1
+                  zua(:,:) = zua(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
+                  zva(:,:) = zva(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)    
+               END DO
+               DO jk = 1, jpkm1
+                  ub(:,:,jk) = ub(:,:,jk) - (zua(:,:) &
+                  & / ( hu_0(:,:) + sshu_n(:,:) + 1._wp - umask(:,:,1) ) - un_b(:,:)) * umask(:,:,jk)
+                  vb(:,:,jk) = vb(:,:,jk) - (zva(:,:) & 
+                  & / ( hv_0(:,:) + sshv_n(:,:) + 1._wp - vmask(:,:,1) ) - vn_b(:,:)) * vmask(:,:,jk)
+               END DO
+            ELSE
+               DO jk = 1, jpkm1
+                  zua(:,:) = zua(:,:) + fse3u(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
+                  zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)    
+               END DO
+               DO jk = 1, jpkm1
+                  ub(:,:,jk) = ub(:,:,jk) - (zua(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)
+                  vb(:,:,jk) = vb(:,:,jk) - (zva(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk)
+               END DO
+            ENDIF
+         ENDIF
+#endif
+         !
+      ENDIF ! neuler =/0
 
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
@@ -278 +368 @@
       ! 
       CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
+      IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva )
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')

branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

@@ -27 +27 @@
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
-   USE solver          ! solver initialization
-   USE wrk_nemo        ! Memory Allocation
-   USE timing          ! Timing
-
+   USE solver         ! solver initialization
+   USE wrk_nemo       ! Memory Allocation
+   USE timing         ! Timing
+   USE dynhpg, ONLY: ln_dynhpg_imp ! semi-implicit hpg flag
 
    IMPLICIT NONE
@@ -99 +99 @@
          ztrdv(:,:,:) = va(:,:,:)
       ENDIF
-
-      IF( ln_apr_dyn ) THEN                   !==  Atmospheric pressure gradient  ==!
+      ! bg jchanut tschanges
+      IF( ( ln_apr_dyn ).AND.(.NOT.lk_dynspg_ts) ) THEN  !==  Atmospheric pressure gradient
+      ! (According to time splitting option, lines below are done latter)
+      ! end jchanut tschanges
          zg_2 = grav * 0.5
          DO jj = 2, jpjm1                          ! gradient of Patm using inverse barometer ssh
@@ -207 +209 @@
          WRITE(numout,*) '     Filtered free surface cst volume       lk_dynspg_flt = ', lk_dynspg_flt
       ENDIF
-
+      !
+      ! bg jchanut tschanges
+      IF( lk_dynspg_ts ) CALL dyn_spg_ts_init( nit000 )
+      ! (do it now, to set nn_baro, used to allocate some arrays later on)
+      ! end jchanut tschanges
+      ! 
       !                        ! allocate dyn_spg arrays
       IF( lk_dynspg_ts ) THEN
@@ -247 +254 @@
       ENDIF
 
-      !                        ! Control of momentum formulation
-      IF( lk_dynspg_ts .AND. lk_vvl ) THEN
-         IF( .NOT.ln_dynadv_vec )   CALL ctl_stop( 'Flux form not implemented for this free surface formulation' )
-      ENDIF
+      ! bg jchanut tschanges
+      !               ! Control of hydrostatic pressure choice
+      IF( lk_dynspg_ts .AND. ln_dynhpg_imp ) THEN
+         CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' )
+      ENDIF
+!      !                        ! Control of momentum formulation
+!      IF( lk_dynspg_ts .AND. lk_vvl ) THEN
+!         IF( .NOT.ln_dynadv_vec )   CALL ctl_stop( 'Flux form not implemented for this free surface formulation' )
+!      ENDIF
+      ! end jchanut tschanges
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_init')

branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

@@ -9 +9 @@
    !!             3.3  ! 2010-09  (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
    !!             3.3  ! 2011-03  (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
+   !!             3.5  ! 2013-07  (J. Chanut) Switch to Forward-backward time stepping
    !!---------------------------------------------------------------------
 #if defined key_dynspg_ts   ||   defined key_esopa
@@ -16 +17 @@
    !!   dyn_spg_ts  : compute surface pressure gradient trend using a time-
    !!                 splitting scheme and add to the general trend 
-   !!   ts_rst      : read/write the time-splitting restart fields in the ocean restart file
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
@@ -24 +24 @@
    USE phycst          ! physical constants
    USE domvvl          ! variable volume
-   USE zdfbfr          ! bottom friction
    USE dynvor          ! vorticity term
-   USE obc_oce         ! Lateral open boundary condition
-   USE obc_par         ! open boundary condition parameters
-   USE obcdta          ! open boundary condition data     
-   USE obcfla          ! Flather open boundary condition  
    USE bdy_par         ! for lk_bdy
    USE bdy_oce         ! Lateral open boundary condition
-   USE bdydta          ! open boundary condition data     
+   USE bdytides        ! open boundary condition data     
    USE bdydyn2d        ! open boundary conditions on barotropic variables
-   USE sbctide
-   USE updtide
+   USE sbctide         ! tides
+   USE updtide         ! tide potential
    USE lib_mpp         ! distributed memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -41 +36 @@
    USE in_out_manager  ! I/O manager
    USE iom             ! IOM library
+   USE restart         ! only for lrst_oce
    USE zdf_oce         ! Vertical diffusion
    USE wrk_nemo        ! Memory Allocation
-   USE timing          ! Timing
+   USE timing          ! Timing    
+   USE sbcapr          ! surface boundary condition: atmospheric pressure
+   USE dynadv, ONLY: ln_dynadv_vec
 
 
@@ -49 +47 @@
    PRIVATE
 
-   PUBLIC dyn_spg_ts        ! routine called by step.F90
-   PUBLIC ts_rst            ! routine called by istate.F90
-   PUBLIC dyn_spg_ts_alloc  ! routine called by dynspg.F90
-
-
+   PUBLIC dyn_spg_ts        ! routine called in dynspg.F90 
+   PUBLIC dyn_spg_ts_alloc  !    "      "     "    "
+   PUBLIC dyn_spg_ts_init   !    "      "     "    "
+
+   ! Potential namelist parameters below to be read in dyn_spg_ts_init
+   LOGICAL,  PUBLIC,  PARAMETER :: ln_bt_fw=.TRUE.        !: Forward integration of barotropic sub-stepping
+   LOGICAL,  PRIVATE, PARAMETER :: ln_bt_av=.TRUE.        !: Time averaging of barotropic variables
+   LOGICAL,  PRIVATE, PARAMETER :: ln_bt_nn_auto=.FALSE.  !: Set number of iterations automatically
+   INTEGER,  PRIVATE, PARAMETER :: nn_bt_flt=1            !: Filter choice
+   REAL(wp), PRIVATE, PARAMETER :: rn_bt_cmax=0.8_wp      !: Max. courant number (used if ln_bt_nn_auto=T)
+   !!INTEGER, PUBLIC :: nn_baro=30                          !: Number of barotropic time steps / baroclinic step (rdt)
+   ! End namelist parameters
+
+   INTEGER, SAVE :: icycle  ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
+   REAL(wp),SAVE :: rdtbt   ! Barotropic time step
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: &
+                    wgtbtp1, &              ! Primary weights used for time filtering of barotropic variables
+                    wgtbtp2                 ! Secondary weights used for time filtering of barotropic variables
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zwz          ! ff/h at F points
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftnw, ftne   ! triad of coriolis parameter
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftsw, ftse   ! (only used with een vorticity scheme)
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_b, vn_b   ! now    averaged velocity
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub_b, vb_b   ! before averaged velocity
+   ! Would be convenient to have arrays below defined whatever the free surface option ?
+   ! These could be computed once for all at the beginning of the each baroclinic time step
+   ! and eventually swapped at the end
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ua_b, va_b     ! after  averaged velocities
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_b, vn_b     ! now    averaged velocities
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub_b, vb_b     ! before averaged velocities
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv, vn_adv ! Advection vel. at "now" barocl. step
+   REAL(wp), PRIVATE,ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub2_b,  vb2_b  ! Advection vel. at "now-0.5" barocl. step
+
+   ! Arrays below are saved to allow testing of the "no time averaging" option
+   ! If this option is not retained, these could be replaced by temporary arrays
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  sshbb_e, sshb_e, & ! Instantaneous barotropic arrays
+                                                   ubb_e, ub_e,     &
+                                                   vbb_e, vb_e
 
    !! * Substitutions
@@ -64 +91 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
-   !! $Id$
+   !! NEMO/OPA 3.5 , NEMO Consortium (2013)
+   !! $Id: dynspg_ts.F90
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -74 +101 @@
       !!                  ***  routine dyn_spg_ts_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( ftnw  (jpi,jpj) , ftne(jpi,jpj) , un_b(jpi,jpj) , vn_b(jpi,jpj) ,     &
-         &      ftsw  (jpi,jpj) , ftse(jpi,jpj) , ub_b(jpi,jpj) , vb_b(jpi,jpj) , STAT= dyn_spg_ts_alloc )
-         !
+      INTEGER :: ierr(3)
+      !!----------------------------------------------------------------------
+      ierr(:) = 0
+
+      ALLOCATE( ub_b(jpi,jpj)  , vb_b(jpi,jpj)   , &
+         &      un_b(jpi,jpj)  , vn_b(jpi,jpj)   , &
+         &      ua_b(jpi,jpj)  , va_b(jpi,jpj)   , &
+         &      un_adv(jpi,jpj), vn_adv(jpi,jpj) , &
+         &      sshb_e(jpi,jpj), sshbb_e(jpi,jpj), &
+         &      ub_e(jpi,jpj)  , vb_e(jpi,jpj)   , &
+         &      ubb_e(jpi,jpj) , vbb_e(jpi,jpj)  , &
+         &      wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), & 
+         &      zwz(jpi,jpj), STAT= ierr(1) )
+
+      IF( ln_bt_fw )      ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj), STAT=ierr(2) )
+
+      IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 
+                             &      ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) )
+
+      dyn_spg_ts_alloc = MAXVAL(ierr(:))
+
       IF( lk_mpp                )   CALL mpp_sum( dyn_spg_ts_alloc )
       IF( dyn_spg_ts_alloc /= 0 )   CALL ctl_warn('dynspg_oce_alloc: failed to allocate arrays')
@@ -82 +127 @@
    END FUNCTION dyn_spg_ts_alloc
 
-
    SUBROUTINE dyn_spg_ts( kt )
       !!----------------------------------------------------------------------
-      !!                  ***  routine dyn_spg_ts  ***
       !!
-      !! ** Purpose :   Compute the now trend due to the surface pressure
-      !!      gradient in case of free surface formulation with time-splitting.
-      !!      Add it to the general trend of momentum equation.
+      !! ** Purpose :   
+      !!      -Compute the now trend due to the explicit time stepping
+      !!      of the quasi-linear barotropic system. Barotropic variables are
+      !!      advanced from internal time steps "n" to "n+1" (if ln_bt_cen=F)
+      !!      or from "n-1" to "n+1" time steps (if ln_bt_cen=T) with a
+      !!      generalized forward-backward (see ref. below) time stepping.
+      !!      -Update the free surface at step "n+1" (ssha, sshu_a, sshv_a).
+      !!      -Compute barotropic advective velocities at step "n" to be used 
+      !!      to advect tracers latter on. These are compliant with discrete
+      !!      continuity equation taken at the baroclinic time steps, thus 
+      !!      ensuring tracers conservation.
       !!
-      !! ** Method  :   Free surface formulation with time-splitting
-      !!      -1- Save the vertically integrated trend. This general trend is
-      !!          held constant over the barotropic integration.
-      !!          The Coriolis force is removed from the general trend as the
-      !!          surface gradient and the Coriolis force are updated within
-      !!          the barotropic integration.
-      !!      -2- Barotropic loop : updates of sea surface height (ssha_e) and 
-      !!          barotropic velocity (ua_e and va_e) through barotropic 
-      !!          momentum and continuity integration. Barotropic former 
-      !!          variables are time averaging over the full barotropic cycle
-      !!          (= 2 * baroclinic time step) and saved in uX_b 
-      !!          and vX_b (X specifying after, now or before).
-      !!      -3- The new general trend becomes :
-      !!          ua = ua - sum_k(ua)/H + ( un_b - ub_b )
+      !! ** Method  :  
       !!
-      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
+      !! ** Action : - Update barotropic velocities: ua_b, va_b
+      !!             - Update trend (ua,va) with barotropic component
+      !!             - Update ssha, sshu_a, sshv_a
+      !!             - Update barotropic advective velocity at kt=now
       !!
-      !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL
+      !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005: 
+      !!              The regional oceanic modeling system (ROMS): 
+      !!              a split-explicit, free-surface,
+      !!              topography-following-coordinate oceanic model. 
+      !!              Ocean Modelling, 9, 347-404. 
       !!---------------------------------------------------------------------
       !
       INTEGER, INTENT(in)  ::   kt   ! ocean time-step index
       !
-      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      INTEGER  ::   icycle           ! local scalar
-      INTEGER  ::   ikbu, ikbv       ! local scalar
-      REAL(wp) ::   zraur, zcoef, z2dt_e, z1_2dt_b, z2dt_bf   ! local scalars
-      REAL(wp) ::   z1_8, zx1, zy1                            !   -      -
-      REAL(wp) ::   z1_4, zx2, zy2                            !   -      -
-      REAL(wp) ::   zu_spg, zu_cor, zu_sld, zu_asp            !   -      -
-      REAL(wp) ::   zv_spg, zv_cor, zv_sld, zv_asp            !   -      -
-      REAL(wp) ::   ua_btm, va_btm                            !   -      -
-      !
-      REAL(wp), POINTER, DIMENSION(:,:) :: zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv 
-      REAL(wp), POINTER, DIMENSION(:,:) :: zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e 
-      REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum
+      LOGICAL  ::   lk_fw_start        ! if true, forward integration 
+      LOGICAL  ::   lk_init         ! if true, special startup of 2d equations
+      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices
+      INTEGER  ::   ikbu, ikbv, noffset      ! local integers
+      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars
+      REAL(wp) ::   zx1, zy1, zx2, zy2          !   -      -
+      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2  !   -      -
+      REAL(wp) ::   zu_spg, zv_spg     !   -      -
+      REAL(wp) ::   zhura, zhvra          !   -      -
+      REAL(wp) ::   za0, za1, za2, za3    !   -      -
+      !
+      REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e
+      REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc
+      REAL(wp), POINTER, DIMENSION(:,:) :: zu_sum, zv_sum, zwx, zwy, zhdiv
+      REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('dyn_spg_ts')
       !
-      CALL wrk_alloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv     )
-      CALL wrk_alloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e   )
-      CALL wrk_alloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum )
-      !
-      IF( kt == nit000 ) THEN             !* initialisation
-         !
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~   free surface with time splitting'
-         IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ',  2*nn_baro
-         !
-         CALL ts_rst( nit000, 'READ' )   ! read or initialize the following fields: un_b, vn_b  
-         !
-         ua_e  (:,:) = un_b (:,:)
-         va_e  (:,:) = vn_b (:,:)
-         hu_e  (:,:) = hu   (:,:)
-         hv_e  (:,:) = hv   (:,:)
-         hur_e (:,:) = hur  (:,:)
-         hvr_e (:,:) = hvr  (:,:)
-         IF( ln_dynvor_een ) THEN
-            ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;   ftsw(1,:) = 0._wp
+      !                                         !* Allocate temporay arrays
+      CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv )
+      CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e  )
+      CALL wrk_alloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc)
+      CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e)
+      !
+      !                                         !* Local constant initialization
+      z1_12 = 1._wp / 12._wp 
+      z1_8  = 0.125_wp                                   
+      z1_4  = 0.25_wp
+      z1_2  = 0.5_wp     
+      zraur = 1._wp / rau0
+      !
+      IF( kt == nit000 .AND. neuler == 0 ) THEN    ! reciprocal of baroclinic time step 
+        z2dt_bf = rdt
+      ELSE
+        z2dt_bf = 2.0_wp * rdt
+      ENDIF
+      z1_2dt_b = 1.0_wp / z2dt_bf 
+      !
+      lk_init = ln_bt_av                           ! if no time averaging, then no specific restart 
+      lk_fw_start = .FALSE.
+      !
+                                                       ! time offset in steps for bdy data update
+      IF (.NOT.ln_bt_fw) THEN ; noffset=-2*nn_baro ; ELSE ;  noffset = 0 ; ENDIF
+      !
+      IF( kt == nit000 ) THEN                !* initialisation
+         !
+         IF (neuler==0) lk_init=.TRUE.
+         !
+         IF (ln_bt_fw.OR.(neuler==0)) THEN
+           lk_fw_start=.TRUE.
+           noffset = 0
+         ELSE
+           lk_fw_start=.FALSE.
+         ENDIF
+         !
+         ! Set averaging weights and cycle length:
+         CALL ts_wgt(ln_bt_av, lk_fw_start, icycle, wgtbtp1, wgtbtp2)
+         !
+         IF ((neuler/=0).AND.(ln_bt_fw)) CALL ts_rst( nit000, 'READ' ) 
+         !
+      ENDIF
+      !
+      ! Set arrays to remove/compute coriolis trend.
+      ! Do it once at kt=nit000 if volume is fixed, else at each long time step.
+      ! Note that these arrays are also used during barotropic loop. These are however frozen
+      ! although they should be updated in variable volume case. Not a big approximation.
+      ! To remove this approximation, copy lines below inside barotropic loop
+      ! and update depths at T-F points (ht and hf resp.) at each barotropic time step
+      !
+      IF ( kt == nit000 .OR. lk_vvl ) THEN
+         IF ( ln_dynvor_een ) THEN
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1
+                  zwz(ji,jj) = ( ht(ji  ,jj+1)*tmask(ji  ,jj+1,1) + & 
+                        &        ht(ji+1,jj+1)*tmask(ji+1,jj+1,1) + &
+                        &        ht(ji  ,jj  )*tmask(ji  ,jj  ,1) + & 
+                        &        ht(ji+1,jj  )*tmask(ji+1,jj  ,1)   &
+                        &      ) * z1_4
+                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zwz(ji,jj)
+               END DO
+            END DO
+            CALL lbc_lnk( zwz, 'F', 1._wp )
+            zwz(:,:) = ff(:,:) * zwz(:,:)
+
+            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
             DO jj = 2, jpj
                DO ji = fs_2, jpi   ! vector opt.
-                  ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3._wp
-                  ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3._wp
-                  ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3._wp
-                  ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3._wp
-               END DO
-            END DO
-         ENDIF
-         !
-      ENDIF
-
-      !                                                     !* Local constant initialization
-      z1_2dt_b = 1._wp / ( 2.0_wp * rdt )                   ! reciprocal of baroclinic time step
-      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt_b = 1.0_wp / rdt    ! reciprocal of baroclinic 
-                                                                        ! time step (euler timestep)
-      z1_8     = 0.125_wp                                   ! coefficient for vorticity estimates
-      z1_4     = 0.25_wp        
-      zraur    = 1._wp / rau0                               ! 1 / volumic mass
-      !
-      zhdiv(:,:) = 0._wp                                    ! barotropic divergence
-      zu_sld = 0._wp   ;   zu_asp = 0._wp                   ! tides trends (lk_tide=F)
-      zv_sld = 0._wp   ;   zv_asp = 0._wp
-
-      IF( kt == nit000 .AND. neuler == 0) THEN              ! for implicit bottom friction
-        z2dt_bf = rdt
-      ELSE
-        z2dt_bf = 2.0_wp * rdt
-      ENDIF
-
+                  ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
+                  ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
+                  ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
+                  ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
+               END DO
+            END DO
+         ELSE
+            zwz(:,:) = 0._wp
+            ! JC: TBC. hf should be greater than 0 
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF( hf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / hf(ji,jj)
+               END DO
+            END DO
+            zwz(:,:) = ff(:,:) * zwz(:,:)
+         ENDIF
+      ENDIF
+      !
+      ! If forward start at previous time step, and centered integration, 
+      ! then update averaging weights:
+      IF ((.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN
+         lk_fw_start=.FALSE.
+         CALL ts_wgt(ln_bt_av, lk_fw_start, icycle, wgtbtp1, wgtbtp2)
+      ENDIF
+                          
       ! -----------------------------------------------------------------------------
       !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
       ! -----------------------------------------------------------------------------
       !      
+      ! Some vertical sums (at now and before time steps) below could be suppressed 
+      ! if one swap barotropic arrays somewhere
+      !
       !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated)
-      !                                   ! --------------------------
-      zua(:,:) = 0._wp   ;   zun(:,:) = 0._wp   ;   ub_b(:,:) = 0._wp
-      zva(:,:) = 0._wp   ;   zvn(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
+      !                                   ! --------------------------------------------------
+      zu_frc(:,:) = 0._wp   ;   ub_b(:,:) = 0._wp  ;  un_b(:,:) = 0._wp
+      zv_frc(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp  ;  vn_b(:,:) = 0._wp
       !
       DO jk = 1, jpkm1
@@ -198 +291 @@
             DO ji = 1, jpi
 #endif
-               !                                                                              ! now trend
-               zua(ji,jj) = zua(ji,jj) + fse3u  (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
-               zva(ji,jj) = zva(ji,jj) + fse3v  (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
-               !                                                                              ! now velocity 
-               zun(ji,jj) = zun(ji,jj) + fse3u  (ji,jj,jk) * un(ji,jj,jk)
-               zvn(ji,jj) = zvn(ji,jj) + fse3v  (ji,jj,jk) * vn(ji,jj,jk)               
-               !
-#if defined key_vvl
-               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk)* ub(ji,jj,jk)   *umask(ji,jj,jk) 
-               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk)* vb(ji,jj,jk)   *vmask(ji,jj,jk)
-#else
-               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_0(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk)
-               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_0(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk)
-#endif
+               !        ! now trend:                                                                   
+               zu_frc(ji,jj) = zu_frc(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
+               zv_frc(ji,jj) = zv_frc(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)         
+               !        ! now bt transp:                   
+               un_b(ji,jj) = un_b(ji,jj) + fse3u(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)        
+               vn_b(ji,jj) = vn_b(ji,jj) + fse3v(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
+               !  ! before bt transp:
+               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk)
+               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk)
             END DO
          END DO
       END DO
-
+      !
+      zu_frc(:,:) = zu_frc(:,:) * hur(:,:)
+      zv_frc(:,:) = zv_frc(:,:) * hvr(:,:)
+      !
+      IF( lk_vvl ) THEN
+          ub_b(:,:) = ub_b(:,:) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) )
+          vb_b(:,:) = vb_b(:,:) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) )
+      ELSE
+          ub_b(:,:) = ub_b(:,:) * hur(:,:)
+          vb_b(:,:) = vb_b(:,:) * hvr(:,:)
+      ENDIF
+      !
       !                                   !* baroclinic momentum trend (remove the vertical mean trend)
-      DO jk = 1, jpkm1                    ! --------------------------
+      DO jk = 1, jpkm1                    ! -----------------------------------------------------------
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
-               va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
+               ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk)
             END DO
          END DO
       END DO
-
-      !                                   !* barotropic Coriolis trends * H (vorticity scheme dependent)
-      !                                   ! ---------------------------====
-      zwx(:,:) = zun(:,:) * e2u(:,:)                   ! now transport 
-      zwy(:,:) = zvn(:,:) * e1v(:,:)
+      !                                   !* barotropic Coriolis trends (vorticity scheme dependent)
+      !                                   ! --------------------------------------------------------
+      zwx(:,:) = un_b(:,:) * e2u(:,:)           ! now transport 
+      zwy(:,:) = vn_b(:,:) * e1v(:,:)
       !
       IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
@@ -239 +337 @@
                zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
                ! energy conserving formulation for planetary vorticity term
-               zcu(ji,jj) = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
-               zcv(ji,jj) =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
-            END DO
-         END DO
-         !
-      ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
+               zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
+               zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
+            END DO
+         END DO
+         !
+      ELSEIF ( ln_dynvor_ens ) THEN             ! enstrophy conserving scheme
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
-               zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
-               zcu(ji,jj)  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
-               zcv(ji,jj)  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
-            END DO
-         END DO
-         !
-      ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme
+               zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
+                 &            + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
+               zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
+                 &            + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
+               zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
+               zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
+            END DO
+         END DO
+         !
+      ELSEIF ( ln_dynvor_een ) THEN             ! enstrophy and energy conserving scheme
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zcu(ji,jj) = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
-                  &                                + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
-               zcv(ji,jj) = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
-                  &                                + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
-            END DO
-         END DO
-         !
-      ENDIF
-
+               zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
+                &                                      + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
+                &                                      + ftse(ji,jj  ) * zwy(ji  ,jj-1) &
+                &                                      + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
+               zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
+                &                                      + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
+                &                                      + ftnw(ji,jj  ) * zwx(ji-1,jj  ) &
+                &                                      + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
+            END DO
+         END DO
+         !
+      ENDIF 
+      !
+      un_b (:,:) = un_b(:,:) * hur(:,:)         ! Revert now transport to barotropic velocities
+      vn_b (:,:) = vn_b(:,:) * hvr(:,:)  
       !                                   !* Right-Hand-Side of the barotropic momentum equation
       !                                   ! ----------------------------------------------------
-      IF( lk_vvl ) THEN                         ! Variable volume : remove both Coriolis and Surface pressure gradient
+      IF( lk_vvl ) THEN                         ! Variable volume : remove surface pressure gradient
          DO jj = 2, jpjm1 
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zcu(ji,jj) = zcu(ji,jj) - grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn(ji+1,jj  )    &
-                  &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hu(ji,jj) / e1u(ji,jj)
-               zcv(ji,jj) = zcv(ji,jj) - grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn(ji  ,jj+1)    &
-                  &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hv(ji,jj) / e2v(ji,jj)
-            END DO
-         END DO
-      ENDIF
-
-      DO jj = 2, jpjm1                             ! Remove coriolis term (and possibly spg) from barotropic trend
+               zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj)
+               zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj)
+            END DO
+         END DO
+      ENDIF
+
+      DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend
          DO ji = fs_2, fs_jpim1
-             zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)
-             zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)
+             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1)
+             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)
           END DO
-      END DO
-
-                    
-      !                                             ! Remove barotropic contribution of bottom friction 
-      !                                             ! from the barotropic transport trend
-      zcoef = -1._wp * z1_2dt_b
-
-      IF(ln_bfrimp) THEN
-      !                                   ! Remove the bottom stress trend from 3-D sea surface level gradient
-      !                                   ! and Coriolis forcing in case of 3D semi-implicit bottom friction 
-        DO jj = 2, jpjm1         
-           DO ji = fs_2, fs_jpim1
-              ikbu = mbku(ji,jj)
-              ikbv = mbkv(ji,jj)
-              ua_btm = zcu(ji,jj) * z2dt_bf * hur(ji,jj) * umask (ji,jj,ikbu)
-              va_btm = zcv(ji,jj) * z2dt_bf * hvr(ji,jj) * vmask (ji,jj,ikbv)
-
-              zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm
-              zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm
-           END DO
-        END DO
-
-      ELSE
-
-# if defined key_vectopt_loop
-        DO jj = 1, 1
-           DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
-# else
-        DO jj = 2, jpjm1
-           DO ji = 2, jpim1
-# endif
-            ! Apply stability criteria for bottom friction
-            !RBbug for vvl and external mode we may need to use varying fse3
-            !!gm  Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx.
-              zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef  )
-              zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef  )
-           END DO
-        END DO
-
-        IF( lk_vvl ) THEN
-           DO jj = 2, jpjm1
-              DO ji = fs_2, fs_jpim1   ! vector opt.
-                 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj)   &
-                    &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1._wp - umask(ji,jj,1) )
-                 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj)   &
-                    &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1._wp - vmask(ji,jj,1) )
-              END DO
-           END DO
-        ELSE
-           DO jj = 2, jpjm1
-              DO ji = fs_2, fs_jpim1   ! vector opt.
-                 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj)
-                 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj)
-              END DO
-           END DO
-        ENDIF
-      END IF    ! end (ln_bfrimp)
-
-                    
-      !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity)
-      !                                   ! -------------------------- 
-      zua(:,:) = zua(:,:) * hur(:,:)
-      zva(:,:) = zva(:,:) * hvr(:,:)
-      !
-      IF( lk_vvl ) THEN
-         ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) )
-         vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) )
-      ELSE
-         ub_b(:,:) = ub_b(:,:) * hur(:,:)
-         vb_b(:,:) = vb_b(:,:) * hvr(:,:)
-      ENDIF
-
+      END DO 
+      !
+      !                 ! Add bottom stress contribution from baroclinic velocities:      
+      IF (ln_bt_fw) THEN
+         DO jj = 2, jpjm1                          
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ikbu = mbku(ji,jj)       
+               ikbv = mbkv(ji,jj)    
+               zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
+               zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
+            END DO
+         END DO
+      ELSE
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ikbu = mbku(ji,jj)       
+               ikbv = mbkv(ji,jj)    
+               zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
+               zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
+            END DO
+         END DO
+      ENDIF
+      !
+      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag
+      zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:)
+      zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:)
+      !       
+      IF (ln_bt_fw) THEN                        ! Add wind forcing
+         zu_frc(:,:) =  zu_frc(:,:) + zraur * utau(:,:) * hur(:,:)
+         zv_frc(:,:) =  zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:)
+      ELSE
+         zu_frc(:,:) =  zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:)
+         zv_frc(:,:) =  zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:)
+      ENDIF  
+      !
+      IF ( ln_apr_dyn ) THEN                    ! Add atm pressure forcing
+         IF (ln_bt_fw) THEN
+            DO jj = 2, jpjm1              
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zu_spg =  grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) /e1u(ji,jj)
+                  zv_spg =  grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj)
+                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
+                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
+               END DO
+            END DO
+         ELSE
+            DO jj = 2, jpjm1              
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zu_spg =  grav * z1_2 * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)    &
+                      &                    + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) /e1u(ji,jj)
+                  zv_spg =  grav * z1_2 * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)    &
+                      &                    + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj)
+                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
+                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
+               END DO
+            END DO
+         ENDIF 
+      ENDIF
+      !                                   !* Right-Hand-Side of the barotropic ssh equation
+      !                                   ! -----------------------------------------------
+      !                                         ! Surface net water flux and rivers
+      IF (ln_bt_fw) THEN
+         zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) )
+      ELSE
+         zssh_frc(:,:) = zraur * z1_2 * (emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:))
+      ENDIF
+#if defined key_asminc
+      !                                         ! Include the IAU weighted SSH increment
+      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
+         zssh_frc(:,:) = zssh_frc(:,:) + ssh_iau(:,:)
+      ENDIF
+#endif
+      !
       ! -----------------------------------------------------------------------
-      !  Phase 2 : Integration of the barotropic equations with time splitting
+      !  Phase 2 : Integration of the barotropic equations 
       ! -----------------------------------------------------------------------
       !
       !                                             ! ==================== !
       !                                             !    Initialisations   !
+      !                                             ! ==================== !  
+      ! Initialize barotropic variables:    
+      IF (ln_bt_fw) THEN                  ! FORWARD integration:  start from NOW fields                             
+         sshn_e (:,:) = sshn (:,:)            
+         zun_e  (:,:) = un_b (:,:)            
+         zvn_e  (:,:) = vn_b (:,:)
+      ELSE                                ! CENTERED integration: start from BEFORE fields
+         sshn_e (:,:) = sshb (:,:)
+         zun_e  (:,:) = ub_b (:,:)         
+         zvn_e  (:,:) = vb_b (:,:)
+      ENDIF
+      !
+      ! Initialize depths:
+      IF ( lk_vvl.AND.(.NOT.ln_bt_fw) ) THEN
+         hu_e   (:,:) = hu_0(:,:) + sshu_b(:,:) 
+         hv_e   (:,:) = hv_0(:,:) + sshv_b(:,:)
+         hur_e  (:,:) = umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) )
+         hvr_e  (:,:) = vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) )
+      ELSE
+         hu_e   (:,:) = hu   (:,:)       
+         hv_e   (:,:) = hv   (:,:) 
+         hur_e  (:,:) = hur  (:,:)    
+         hvr_e  (:,:) = hvr  (:,:)
+      ENDIF
+      !
+      IF (.NOT.lk_vvl) THEN ! Depths at jn+0.5:
+         zhup2_e (:,:) = hu(:,:)
+         zhvp2_e (:,:) = hv(:,:)
+      ENDIF
+      !
+      ! Initialize sums:
+      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)          
+      va_b  (:,:) = 0._wp
+      ssha  (:,:) = 0._wp       ! Sum for after averaged sea level
+      zu_sum(:,:) = 0._wp       ! Sum for now transport issued from ts loop
+      zv_sum(:,:) = 0._wp
       !                                             ! ==================== !
-      icycle = 2  * nn_baro            ! Number of barotropic sub time-step
-      
-      !                                ! Start from NOW field
-      hu_e   (:,:) = hu   (:,:)            ! ocean depth at u- and v-points
-      hv_e   (:,:) = hv   (:,:) 
-      hur_e  (:,:) = hur  (:,:)            ! ocean depth inverted at u- and v-points
-      hvr_e  (:,:) = hvr  (:,:)
-!RBbug     zsshb_e(:,:) = sshn (:,:)  
-      zsshb_e(:,:) = sshn_b(:,:)           ! sea surface height (before and now)
-      sshn_e (:,:) = sshn (:,:)
-      
-      zun_e  (:,:) = un_b (:,:)            ! barotropic velocity (external)
-      zvn_e  (:,:) = vn_b (:,:)
-      zub_e  (:,:) = un_b (:,:)
-      zvb_e  (:,:) = vn_b (:,:)
-
-      zu_sum  (:,:) = un_b (:,:)           ! summation
-      zv_sum  (:,:) = vn_b (:,:)
-      zssh_sum(:,:) = sshn (:,:)
-
-#if defined key_obc
-      ! set ssh corrections to 0
-      ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
-      IF( lp_obc_east  )   sshfoe_b(:,:) = 0._wp
-      IF( lp_obc_west  )   sshfow_b(:,:) = 0._wp
-      IF( lp_obc_south )   sshfos_b(:,:) = 0._wp
-      IF( lp_obc_north )   sshfon_b(:,:) = 0._wp
-#endif
-
-      !                                             ! ==================== !
-      DO jn = 1, icycle                             !  sub-time-step loop  ! (from NOW to AFTER+1)
+      DO jn = 1, icycle                             !  sub-time-step loop  !
          !                                          ! ==================== !
-         z2dt_e = 2. * ( rdt / nn_baro )
-         IF( jn == 1 )   z2dt_e = rdt / nn_baro
-
          !                                                !* Update the forcing (BDY and tides)
          !                                                !  ------------------
-         IF( lk_obc )   CALL obc_dta_bt ( kt, jn   )
-         IF( lk_bdy )   CALL bdy_dta ( kt, jit=jn, time_offset=+1 )
-         IF ( ln_tide_pot .AND. lk_tide) CALL upd_tide( kt, jn )
-
-         !                                                !* after ssh_e
+         ! Update only tidal forcing at open boundaries
+#if defined key_tide
+         IF ( lk_bdy .AND. lk_tide )      CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) )
+         IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, time_offset=noffset )
+#endif
+         !
+         ! Set extrapolation coefficients for predictor step:
+         IF ((jn<3).AND.lk_init) THEN      ! Forward           
+           za1 = 1._wp                                          
+           za2 = 0._wp                        
+           za3 = 0._wp                        
+         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105 
+           za1 =  1.781105_wp              ! za1 =   3/2 +   bet
+           za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet)
+           za3 =  0.281105_wp              ! za3 = bet
+         ENDIF
+
+         ! Extrapolate barotropic velocities at step jit+0.5:
+         ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
+         va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
+
+         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only)
+            !                                             !  ------------------
+            ! Extrapolate Sea Level at step jit+0.5:
+            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
+            !
+            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
+               DO ji = 2, fs_jpim1   ! Vector opt.
+                  zwx(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       &
+                     &              * ( e1t(ji  ,jj) * e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
+                     &              +   e1t(ji+1,jj) * e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
+                  zwy(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       &
+                     &              * ( e1t(ji,jj  ) * e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
+                     &              +   e1t(ji,jj+1) * e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
+               END DO
+            END DO
+            CALL lbc_lnk( zwx, 'U', 1._wp )    ;   CALL lbc_lnk( zwy, 'V', 1._wp )
+            !
+            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)               ! Ocean depth at U- and V-points
+            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
+         ENDIF
+         !                                                !* after ssh
          !                                                !  -----------
-         DO jj = 2, jpjm1                                 ! Horizontal divergence of barotropic transports
+         ! One should enforce volume conservation at open boundaries here
+         ! considering fluxes below:
+         !
+         zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5
+         zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
+         DO jj = 2, jpjm1                                 
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj)     &
-                  &             - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj)     &
-                  &             + e1v(ji,jj  ) * zvn_e(ji,jj  ) * hv_e(ji,jj  )     &
-                  &             - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1)   ) / ( e1t(ji,jj) * e2t(ji,jj) )
-            END DO
-         END DO
-         !
-#if defined key_obc
-         !                                                     ! OBC : zhdiv must be zero behind the open boundary
-!!  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
-         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0._wp      ! east
-         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0._wp      ! west
-         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0._wp      ! north
-         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0._wp      ! south
+               zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   &
+                  &             + zwy(ji,jj) - zwy(ji,jj-1)   &
+                  &           ) / ( e1t(ji,jj) * e2t(ji,jj) )
+            END DO
+         END DO
+         !
+         ! Sum over sub-time-steps to compute advective velocities
+         za2 = wgtbtp2(jn)
+         zu_sum  (:,:) = zu_sum  (:,:) + za2 * ua_e  (:,:) * zhup2_e  (:,:)
+         zv_sum  (:,:) = zv_sum  (:,:) + za2 * va_e  (:,:) * zhvp2_e  (:,:)
+         !
+         ! Set next sea level:
+         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
+         CALL lbc_lnk( ssha_e, 'T',  1._wp )
+
+#if defined key_bdy
+         ! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.)
+         IF (lk_bdy) CALL bdy_ssh( ssha_e )
 #endif
-#if defined key_bdy
-         zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:)               ! BDY mask
-#endif
-         !
-         DO jj = 2, jpjm1                                      ! leap-frog on ssh_e
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 
-            END DO
-         END DO
-
-         !                                                !* after barotropic velocities (vorticity scheme dependent)
-         !                                                !  ---------------------------  
-         zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:)     ! now_e transport
-         zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:)
+         !  
+         ! Sea Surface Height at u-,v-points (vvl case only)
+         IF ( lk_vvl ) THEN                                
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1      ! NO Vector Opt.
+                  sshu_a(ji,jj) = z1_2  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                 &
+                     &                                   * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha_e(ji  ,jj) &
+                     &                                     + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha_e(ji+1,jj) )
+                  sshv_a(ji,jj) = z1_2  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                 &
+                     &                                   * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha_e(ji,jj  ) &
+                     &                                     + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha_e(ji,jj+1) )
+               END DO
+            END DO
+            CALL lbc_lnk( sshu_a, 'U', 1._wp )   ;   CALL lbc_lnk( sshv_a, 'V', 1._wp )
+         ENDIF   
+         !                                 
+         ! Half-step back interpolation of SSH for surface pressure computation:
+         !----------------------------------------------------------------------
+         IF ((jn==1).AND.lk_init) THEN
+           za0=1._wp                        ! Forward-backward
+           za1=0._wp                           
+           za2=0._wp
+           za3=0._wp
+         ELSEIF ((jn==2).AND.lk_init) THEN  ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
+           za0= 1.0833333333333_wp          ! za0 = 1-gam-eps
+           za1=-0.1666666666666_wp          ! za1 = gam
+           za2= 0.0833333333333_wp          ! za2 = eps
+           za3= 0._wp              
+         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 
+           za0=0.614_wp                     ! za0 = 1/2 +   gam + 2*eps    
+           za1=0.285_wp                     ! za1 = 1/2 - 2*gam - 3*eps 
+           za2=0.088_wp                     ! za2 = gam
+           za3=0.013_wp                     ! za3 = eps
+         ENDIF
+
+         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) &
+          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
+
+         !
+         ! Compute associated depths at U and V points:
+         IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN      
+            !                                        
+            DO jj = 2, jpjm1                            
+               DO ji = 2, jpim1
+                  zx1 = z1_2 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       &
+                    &        * ( e1t(ji  ,jj) * e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
+                    &        +   e1t(ji+1,jj) * e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )                 
+                  zy1 = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       &
+                     &       * ( e1t(ji,jj  ) * e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
+                     &       +   e1t(ji,jj+1) * e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
+                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1 
+                  zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
+               END DO
+            END DO
+         ENDIF
+         !
+         ! Add Coriolis trend:
+         ! zwz array below or triads normally depend on sea level with key_vvl and should be updated
+         ! at each time step. We however keep them constant here for optimization.
+         ! Recall that zwx and zwy arrays hold fluxes at this stage:
+         ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5
+         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
          !
          IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      !==  energy conserving or mixed scheme  ==!
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ! surface pressure gradient
-                  IF( lk_vvl) THEN
-                     zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
-                        &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
-                     zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
-                        &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
-                  ELSE
-                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
-                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
-                  ENDIF
-                  ! add tidal astronomical forcing
-                  IF ( ln_tide_pot .AND. lk_tide ) THEN 
-                  zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
-                  zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
-                  ENDIF
-                  ! energy conserving formulation for planetary vorticity term
                   zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
                   zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
                   zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
                   zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
-                  zu_cor = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj)
-                  zv_cor =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj)
-                  ! after velocities with implicit bottom friction
-
-                  IF( ln_bfrimp ) THEN      ! implicit bottom friction
-                     !   A new method to implement the implicit bottom friction. 
-                     !   H. Liu
-                     !   Sept 2011
-                     ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            &
-                      &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            &
-                      &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) )
-                     ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)   
-                     !
-                     va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            &
-                      &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            &
-                      &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) )
-                     va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)   
-                     !
-                  ELSE
-                     ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1)   &
-                      &           / ( 1._wp         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
-                     va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1)   &
-                      &           / ( 1._wp         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
-                  ENDIF
+                  zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
+                  zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
                END DO
             END DO
@@ -491 +655 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                   ! surface pressure gradient
-                  IF( lk_vvl) THEN
-                     zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
-                        &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
-                     zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
-                        &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
-                  ELSE
-                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
-                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
-                  ENDIF
-                  ! add tidal astronomical forcing
-                  IF ( ln_tide_pot .AND. lk_tide ) THEN
-                  zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
-                  zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
-                  ENDIF
-                  ! enstrophy conserving formulation for planetary vorticity term
-                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
-                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
-                  zu_cor  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj)
-                  zv_cor  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj)
-                  ! after velocities with implicit bottom friction
-                  IF( ln_bfrimp ) THEN
-                     !   A new method to implement the implicit bottom friction. 
-                     !   H. Liu
-                     !   Sept 2011
-                     ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            &
-                      &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            &
-                      &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) )
-                     ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)   
-                     !
-                     va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            &
-                      &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            &
-                      &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) )
-                     va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)   
-                     !
-                  ELSE
-                     ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1)   &
-                     &            / ( 1._wp        - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
-                     va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1)   &
-                     &            / ( 1._wp        - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
-                  ENDIF
+                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
+                   &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
+                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
+                   &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
+                  zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
+                  zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
                END DO
             END DO
@@ -538 +667 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ! surface pressure gradient
-                  IF( lk_vvl) THEN
-                     zu_spg = -grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  )    &
-                        &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
-                     zv_spg = -grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
-                        &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
-                  ELSE
-                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
-                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
-                  ENDIF
-                  ! add tidal astronomical forcing
-                  IF ( ln_tide_pot .AND. lk_tide ) THEN
-                  zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
-                  zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
-                  ENDIF
-                  ! energy/enstrophy conserving formulation for planetary vorticity term
-                  zu_cor = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
-                     &                           + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj)
-                  zv_cor = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
-                     &                           + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) * hvr_e(ji,jj)
-                  ! after velocities with implicit bottom friction
-                  IF( ln_bfrimp ) THEN
-                     !   A new method to implement the implicit bottom friction. 
-                     !   H. Liu
-                     !   Sept 2011
-                     ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            &
-                      &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            &
-                      &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) )
-                     ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)   
-                     !
-                     va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            &
-                      &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            &
-                      &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) )
-                     va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)   
-                     !
-                  ELSE
-                     ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1)   &
-                     &            / ( 1._wp        - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
-                     va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1)   &
-                     &            / ( 1._wp        - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
-                  ENDIF
+                  zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
+                     &                                    + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
+                     &                                    + ftse(ji,jj  ) * zwy(ji  ,jj-1) & 
+                     &                                    + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
+                  zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 
+                     &                                    + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
+                     &                                    + ftnw(ji,jj  ) * zwx(ji-1,jj  ) & 
+                     &                                    + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
                END DO
             END DO
             ! 
          ENDIF
-         !                                                     !* domain lateral boundary
-         !                                                     !  -----------------------
-
-                                                               ! OBC open boundaries
-         IF( lk_obc               )   CALL obc_fla_ts ( ua_e, va_e, sshn_e, ssha_e )
-
-                                                               ! BDY open boundaries
-#if defined key_bdy
-         pssh => sshn_e
+         !
+         ! Add tidal astronomical forcing if defined
+         IF ( lk_tide.AND.ln_tide_pot ) THEN
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
+                  zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
+                  zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
+                  zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
+               END DO
+            END DO
+         ENDIF
+         !
+         ! Add bottom stresses:
+         zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)
+         zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)
+         !
+         ! Surface pressure trend:
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ! Add surface pressure gradient
+               zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
+               zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
+               zwx(ji,jj) = zu_spg
+               zwy(ji,jj) = zv_spg
+            END DO
+         END DO
+         !
+         ! Set next velocities:
+         IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN    ! Vector form
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ua_e(ji,jj) = (                                zun_e(ji,jj)   & 
+                            &     + rdtbt * (                      zwx(ji,jj)   &
+                            &                                 + zu_trd(ji,jj)   &
+                            &                                 + zu_frc(ji,jj) ) & 
+                            &   ) * umask(ji,jj,1)
+
+                  va_e(ji,jj) = (                                zvn_e(ji,jj)   &
+                            &     + rdtbt * (                      zwy(ji,jj)   &
+                            &                                 + zv_trd(ji,jj)   &
+                            &                                 + zv_frc(ji,jj) ) &
+                            &   ) * vmask(ji,jj,1)
+               END DO
+            END DO
+
+         ELSE                 ! Flux form
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+
+                  zhura = umask(ji,jj,1)/(hu_0(ji,jj) + sshu_a(ji,jj) + 1._wp - umask(ji,jj,1))
+                  zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + sshv_a(ji,jj) + 1._wp - vmask(ji,jj,1))
+
+                  ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   & 
+                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
+                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   &
+                            &               +      hu(ji,jj)  * zu_frc(ji,jj) ) &
+                            &   ) * zhura
+
+                  va_e(ji,jj) = (                hv_e(ji,jj)  *  zvn_e(ji,jj)   &
+                            &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   &
+                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   &
+                            &               +      hv(ji,jj)  * zv_frc(ji,jj) ) &
+                            &   ) * zhvra
+               END DO
+            END DO
+         ENDIF
+         !
+         IF( lk_vvl ) THEN                             !* Update ocean depth (variable volume case only)
+            !                                          !  ----------------------------------------------        
+            hu_e (:,:) = hu_0(:,:) + sshu_a(:,:)
+            hv_e (:,:) = hv_0(:,:) + sshv_a(:,:)
+            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
+            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
+            !
+         ENDIF
+         !                                                 !* domain lateral boundary
+         !                                                 !  -----------------------
+         !
+         CALL lbc_lnk( ua_e  , 'U', -1._wp )               ! local domain boundaries 
+         CALL lbc_lnk( va_e  , 'V', -1._wp )
+
+#if defined key_bdy  
+
+         pssh => ssha_e
          phur => hur_e
          phvr => hvr_e
          pu2d => ua_e
          pv2d => va_e
-
-         IF( lk_bdy )   CALL bdy_dyn2d( kt ) 
+                                      
+         IF( lk_bdy )   CALL bdy_dyn2d( kt )               ! open boundaries
 #endif
-
-         !
-         CALL lbc_lnk( ua_e  , 'U', -1. )                      ! local domain boundaries 
-         CALL lbc_lnk( va_e  , 'V', -1. )
-         CALL lbc_lnk( ssha_e, 'T',  1. )
-
-         zu_sum  (:,:) = zu_sum  (:,:) + ua_e  (:,:)           ! Sum over sub-time-steps
-         zv_sum  (:,:) = zv_sum  (:,:) + va_e  (:,:) 
-         zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:) 
-
-         !                                                !* Time filter and swap
-         !                                                !  --------------------
-         IF( jn == 1 ) THEN                                     ! Swap only (1st Euler time step)
-            zsshb_e(:,:) = sshn_e(:,:)
-            zub_e  (:,:) = zun_e (:,:)
-            zvb_e  (:,:) = zvn_e (:,:)
-            sshn_e (:,:) = ssha_e(:,:)
-            zun_e  (:,:) = ua_e  (:,:)
-            zvn_e  (:,:) = va_e  (:,:)
-         ELSE                                                   ! Swap + Filter
-            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
-            zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e (:,:)
-            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e (:,:)
-            sshn_e (:,:) = ssha_e(:,:)
-            zun_e  (:,:) = ua_e  (:,:)
-            zvn_e  (:,:) = va_e  (:,:)
-         ENDIF
-
-         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only)
-            !                                             !  ------------------
-            DO jj = 1, jpjm1                                    ! Sea Surface Height at u- & v-points
-               DO ji = 1, fs_jpim1   ! Vector opt.
-                  zsshun_e(ji,jj) = 0.5_wp * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       &
-                     &                     * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    &
-                     &                     +   e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )
-                  zsshvn_e(ji,jj) = 0.5_wp * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       &
-                     &                     * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    &
-                     &                     +   e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )
-               END DO
-            END DO
-            CALL lbc_lnk( zsshun_e, 'U', 1. )                   ! lateral boundaries conditions
-            CALL lbc_lnk( zsshvn_e, 'V', 1. ) 
-            !
-            hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:)              ! Ocean depth at U- and V-points
-            hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:)
-            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
-            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
-            !
-         ENDIF
+         !                                             !* Swap
+         !                                             !  ----
+         ubb_e  (:,:) = ub_e  (:,:)
+         ub_e   (:,:) = zun_e (:,:)
+         zun_e  (:,:) = ua_e  (:,:)
+         !
+         vbb_e  (:,:) = vb_e  (:,:)
+         vb_e   (:,:) = zvn_e (:,:)
+         zvn_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. (.NOT. lk_vvl)) THEN    ! Sum velocities
+            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) 
+            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) 
+         ELSE                                ! Sum transports
+            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:)
+            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:)
+         ENDIF
+         !                                   ! Sum sea level
+         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
          !                                                 ! ==================== !
       END DO                                               !        end loop      !
       !                                                    ! ==================== !
-
-#if defined key_obc
-      IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)     !!gm totally useless ?????
-      IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:)
-      IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:)
-      IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:)
-#endif
-
       ! -----------------------------------------------------------------------------
       ! Phase 3. update the general trend with the barotropic trend
       ! -----------------------------------------------------------------------------
       !
-      !                                   !* Time average ==> after barotropic u, v, ssh
-      zcoef =  1._wp / ( 2 * nn_baro  + 1 ) 
-      zu_sum(:,:) = zcoef * zu_sum  (:,:) 
-      zv_sum(:,:) = zcoef * zv_sum  (:,:) 
-      ! 
-      !                                   !* update the general momentum trend
-      DO jk=1,jpkm1
-         ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b
-         va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b
+      ! At this stage ssha holds a time averaged value
+      !                                                ! Sea Surface Height at u-,v- and f-points
+      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1      ! NO Vector Opt.
+               sshu_a(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
+                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
+                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
+               sshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
+                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
+                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
+            END DO
+         END DO
+         CALL lbc_lnk( sshu_a, 'U', 1._wp )   ;   CALL lbc_lnk( sshv_a, 'V', 1._wp ) ! Boundary conditions
+      ENDIF
+      !
+      ! Set advection velocity correction:
+      IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN     
+         un_adv(:,:) = zu_sum(:,:)*hur(:,:)
+         vn_adv(:,:) = zv_sum(:,:)*hvr(:,:)
+      ELSE
+         un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:)
+         vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:)
+      END IF
+
+      IF (ln_bt_fw) THEN ! Save integrated transport for next computation
+         ub2_b(:,:) = zu_sum(:,:)
+         vb2_b(:,:) = zv_sum(:,:)
+      ENDIF
+      !
+      ! Update barotropic trend:
+      IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN
+         DO jk=1,jpkm1
+            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
+            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
+         END DO
+      ELSE
+         hu_e(:,:) = hu_0(:,:) + sshu_b(:,:)
+         hv_e(:,:) = hv_0(:,:) + sshv_b(:,:)
+         DO jk=1,jpkm1
+            ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_e(:,:) ) * z1_2dt_b
+            va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_e(:,:) ) * z1_2dt_b
+         END DO
+         ! Save barotropic velocities not transport:
+         ua_b  (:,:) =  ua_b(:,:) / ( hu_0(:,:) + sshu_a(:,:) + 1._wp - umask(:,:,1) )
+         va_b  (:,:) =  va_b(:,:) / ( hv_0(:,:) + sshv_a(:,:) + 1._wp - vmask(:,:,1) )
+      ENDIF
+      !
+      !                                   !* write time-spliting arrays in the restart
+      IF(lrst_oce .AND.ln_bt_fw)   CALL ts_rst( kt, 'WRITE' )
+      !
+      CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv )
+      CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e )
+      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc )
+      CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts')
+      !
+   END SUBROUTINE dyn_spg_ts
+
+   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE ts_wgt  ***
+      !!
+      !! ** 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_baro), INTENT(inout) ::   zwgt1, & ! Primary weights
+                                                         zwgt2    ! Secondary weights
+      
+      INTEGER ::  jic, jn, ji                      ! temporary integers
+      REAL(wp) :: za1, za2
+      !!----------------------------------------------------------------------
+
+      zwgt1(:) = 0._wp
+      zwgt2(:) = 0._wp
+
+      ! Set time index when averaged value is requested
+      IF (ll_fw) THEN 
+         jic = nn_baro
+      ELSE
+         jic = 2 * nn_baro
+      ENDIF
+
+      ! Set primary weights:
+      IF (ll_av) THEN
+           ! Define simple boxcar window for primary weights 
+           ! (width = nn_baro, centered around jic)     
+         SELECT CASE ( nn_bt_flt )
+              CASE( 0 )  ! No averaging
+                 zwgt1(jic) = 1._wp
+                 jpit = jic
+
+              CASE( 1 )  ! Boxcar, width = nn_baro
+                 DO jn = 1, 3*nn_baro
+                    za1 = ABS(float(jn-jic))/float(nn_baro) 
+                    IF (za1 < 0.5_wp) THEN
+                      zwgt1(jn) = 1._wp
+                      jpit = jn
+                    ENDIF
+                 ENDDO
+
+              CASE( 2 )  ! Boxcar, width = 2 * nn_baro
+                 DO jn = 1, 3*nn_baro
+                    za1 = ABS(float(jn-jic))/float(nn_baro) 
+                    IF (za1 < 1._wp) THEN
+                      zwgt1(jn) = 1._wp
+                      jpit = jn
+                    ENDIF
+                 ENDDO
+              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
+         END SELECT
+
+      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
       END DO
-      un_b  (:,:) =  zu_sum(:,:) 
-      vn_b  (:,:) =  zv_sum(:,:) 
-      sshn_b(:,:) = zcoef * zssh_sum(:,:) 
-      !
-      !                                   !* write time-spliting arrays in the restart
-      IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' )
-      !
-      CALL wrk_dealloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv     )
-      CALL wrk_dealloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e   )
-      CALL wrk_dealloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts')
-      !
-   END SUBROUTINE dyn_spg_ts
-
+
+      ! 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
+      END DO
+      !
+   END SUBROUTINE ts_wgt
 
    SUBROUTINE ts_rst( kt, cdrw )
@@ -698 +950 @@
       CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
       !
-      INTEGER ::  ji, jk        ! dummy loop indices
       !!----------------------------------------------------------------------
       !
       IF( TRIM(cdrw) == 'READ' ) THEN
-         IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN
-            CALL iom_get( numror, jpdom_autoglo, 'un_b'  , un_b  (:,:) )   ! external velocity issued
-            CALL iom_get( numror, jpdom_autoglo, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
-         ELSE
-            un_b (:,:) = 0._wp
-            vn_b (:,:) = 0._wp
-            ! vertical sum
-            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
-               DO jk = 1, jpkm1
-                  DO ji = 1, jpij
-                     un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
-                     vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
-                  END DO
-               END DO
-            ELSE                             ! No  vector opt.
-               DO jk = 1, jpkm1
-                  un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
-                  vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
-               END DO
-            ENDIF
-            un_b (:,:) = un_b(:,:) * hur(:,:)
-            vn_b (:,:) = vn_b(:,:) * hvr(:,:)
-         ENDIF
-
-         ! Vertically integrated velocity (before)
-         IF (neuler/=0) THEN
-            ub_b (:,:) = 0._wp
-            vb_b (:,:) = 0._wp
-
-            ! vertical sum
-            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
-               DO jk = 1, jpkm1
-                  DO ji = 1, jpij
-                     ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk)
-                     vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk)
-                  END DO
-               END DO
-            ELSE                             ! No  vector opt.
-               DO jk = 1, jpkm1
-                  ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk)
-                  vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk)
-               END DO
-            ENDIF
-
-            IF( lk_vvl ) THEN
-               ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) )
-               vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) )
-            ELSE
-               ub_b(:,:) = ub_b(:,:) * hur(:,:)
-               vb_b(:,:) = vb_b(:,:) * hvr(:,:)
-            ENDIF
-         ELSE                                 ! neuler==0
-            ub_b (:,:) = un_b (:,:)
-            vb_b (:,:) = vn_b (:,:)
-         ENDIF
-
-         IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN
-            CALL iom_get( numror, jpdom_autoglo, 'sshn_b' , sshn_b (:,:) )   ! filtered ssh
-         ELSE
-            sshn_b(:,:) = sshb(:,:)   ! if not in restart set previous time mean to current baroclinic before value   
-         ENDIF 
+         CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:) )   
+         CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:) ) 
+         IF (.NOT.ln_bt_av) THEN
+            CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:) )   
+            CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:) )   
+            CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:) )
+            CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:) ) 
+            CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:) )   
+            CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:) )
+         ENDIF
+      !
       ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
-         CALL iom_rstput( kt, nitrst, numrow, 'un_b'   , un_b  (:,:) )   ! external velocity and ssh
-         CALL iom_rstput( kt, nitrst, numrow, 'vn_b'   , vn_b  (:,:) )   ! issued from barotropic loop
-         CALL iom_rstput( kt, nitrst, numrow, 'sshn_b' , sshn_b(:,:) )   ! 
+         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:) )
+         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:) )
+         !
+         IF (.NOT.ln_bt_av) THEN
+            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) ) 
+            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:) )
+         ENDIF
       ENDIF
       !
    END SUBROUTINE ts_rst
+
+   SUBROUTINE dyn_spg_ts_init( kt )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE dyn_spg_ts_init  ***
+      !!
+      !! ** Purpose : Set time splitting options
+      !!----------------------------------------------------------------------
+      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
+      !
+      INTEGER  :: ji ,jj
+      REAL(wp) :: zxr2, zyr2, zcmax
+      REAL(wp), POINTER, DIMENSION(:,:) :: zcu
+      !!
+!      NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &
+!      &                  nn_baro, rn_bt_cmax, nn_bt_flt
+      !!----------------------------------------------------------------------
+!      REWIND( numnam )              !* Namelist namsplit: split-explicit free surface
+!      READ  ( numnam, namsplit )
+      !         ! Max courant number for ext. grav. waves
+      !
+      CALL wrk_alloc( jpi, jpj, zcu )
+      !
+      IF (lk_vvl) THEN
+         DO jj = 1, jpj
+            DO ji =1, jpi 
+               zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
+               zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
+               zcu(ji,jj) = sqrt(grav*ht_0(ji,jj)*(zxr2 + zyr2) )
+            END DO
+         END DO
+      ELSE
+         DO jj = 1, jpj
+            DO ji =1, jpi 
+               zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
+               zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
+               zcu(ji,jj) = sqrt( grav*ht(ji,jj)*(zxr2 + zyr2) )
+            END DO
+         END DO
+      ENDIF
+
+      zcmax = MAXVAL(zcu(:,:))
+      IF( lk_mpp )   CALL mpp_max( zcmax )
+
+      ! Estimate number of iterations to satisfy a max courant number=0.8 
+      IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
+      
+      rdtbt = rdt / FLOAT(nn_baro)
+      zcmax = zcmax * rdtbt
+                     ! Print results
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
+      IF( ln_bt_nn_auto ) THEN
+         IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.true. Automatically set nn_baro '
+         IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
+      ELSE
+         IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.false.: Use nn_baro in namelist '
+      ENDIF
+      IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
+      IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
+      IF(lwp) WRITE(numout,*) ' Maximum Courant number is   :', zcmax
+
+      IF(ln_bt_av) THEN
+         IF(lwp) WRITE(numout,*) ' ln_bt_av=.true.  => Time averaging over nn_baro time steps is on '
+      ELSE
+         IF(lwp) WRITE(numout,*) ' ln_bt_av=.false. => No time averaging of barotropic variables '
+      ENDIF
+      !
+      !
+      IF(ln_bt_fw) THEN
+         IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true.  => Forward integration of barotropic variables '
+      ELSE
+         IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centered integration of barotropic variables '
+      ENDIF
+      !
+      IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
+      SELECT CASE ( nn_bt_flt )
+           CASE( 0 )  ;   IF(lwp) WRITE(numout,*) '      Dirac'
+           CASE( 1 )  ;   IF(lwp) WRITE(numout,*) '      Boxcar: width = nn_baro'
+           CASE( 2 )  ;   IF(lwp) WRITE(numout,*) '      Boxcar: width = 2*nn_baro' 
+           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
+      END SELECT
+      !
+      IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN
+         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
+      ENDIF
+      IF ( zcmax>0.9_wp ) THEN
+         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )          
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi, jpj, zcu )
+      !
+   END SUBROUTINE dyn_spg_ts_init
 
 #else
@@ -775 +1077 @@
    !!   Default case :   Empty module   No standart free surface cst volume
    !!----------------------------------------------------------------------
+
+   LOGICAL, PUBLIC, PARAMETER :: ln_bt_fw=.FALSE. ! Forward integration of barotropic sub-stepping
+
 CONTAINS
+
    INTEGER FUNCTION dyn_spg_ts_alloc()    ! Dummy function
       dyn_spg_ts_alloc = 0
@@ -782 +1088 @@
       INTEGER, INTENT(in) :: kt
       WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt
-   END SUBROUTINE dyn_spg_ts
+   END SUBROUTINE dyn_spg_ts 
    SUBROUTINE ts_rst( kt, cdrw )          ! Empty routine
       INTEGER         , INTENT(in) ::   kt         ! ocean time-step
       CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
       WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw
-   END SUBROUTINE ts_rst    
+   END SUBROUTINE ts_rst  
+   SUBROUTINE dyn_spg_ts_init( kt )       ! Empty routine
+      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
+      WRITE(*,*) 'dyn_spg_ts_init   : You should not have seen this print! error?', kt
+   END SUBROUTINE dyn_spg_ts_init
 #endif
    

branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

@@ -16 +16 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
+   USE domvvl          ! variable volume
    USE sbc_oce         ! surface boundary condition: ocean
    USE zdf_oce         ! ocean vertical physics
@@ -24 +25 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE dynadv          ! dynamics: vector invariant versus flux form
+   USE dynspg_oce, ONLY: lk_dynspg_ts
+   USE dynspg_ts
 
    IMPLICIT NONE
@@ -29 +33 @@
 
    PUBLIC   dyn_zdf_imp   ! called by step.F90
+
+   REAL(wp) ::  r_vvl     ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise 
 
    !! * Substitutions
@@ -64 +70 @@
       INTEGER  ::   ikbu, ikbv   ! local integers
       REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars
+      REAL(wp) ::   ze3ua, ze3va
       !!----------------------------------------------------------------------
 
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws
-      REAL(wp), POINTER, DIMENSION(:,:)   ::  zavmu, zavmv
       !!----------------------------------------------------------------------
       !
@@ -73 +79 @@
       !
       CALL wrk_alloc( jpi,jpj,jpk, zwi, zwd, zws ) 
-      CALL wrk_alloc( jpi,jpj, zavmu, zavmv ) 
       !
       IF( kt == nit000 ) THEN
@@ -79 +84 @@
          IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
+         !
+         IF( lk_vvl ) THEN   ;    r_vvl = 1._wp       ! Variable volume indicator
+         ELSE                ;    r_vvl = 0._wp       
+         ENDIF
       ENDIF
 
@@ -94 +103 @@
       IF( ln_bfrimp ) THEN
 # if defined key_vectopt_loop
-      DO jj = 1, 1
-         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
+         DO jj = 1, 1
+            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
 # else
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
 # endif
-            ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
-            ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
-            zavmu(ji,jj) = avmu(ji,jj,ikbu+1)
-            zavmv(ji,jj) = avmv(ji,jj,ikbv+1)
-            avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 
-            avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
-         END DO
-      END DO
-      ENDIF
+               ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points 
+               ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
+               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)
+               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
+            END DO
+         END DO
+      ENDIF
+
+#if defined key_dynspg_ts
+      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
+         DO jk = 1, jpkm1
+            ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk)
+            va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk)
+         END DO
+      ELSE                                            ! applied on thickness weighted velocity
+         DO jk = 1, jpkm1
+            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
+               &           + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
+               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
+            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
+               &           + p2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
+               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
+         END DO
+      ENDIF
+
+      IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN
+         ! remove barotropic velocities:
+         DO jk = 1, jpkm1
+            ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk)
+            va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk)
+         ENDDO
+         ! Add bottom stress due to barotropic component only:
+         DO jj = 2, jpjm1        
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
+               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
+               ze3ua =  ( 1. - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu)
+               ze3va =  ( 1. - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv)
+               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua
+               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va
+            END DO
+         END DO
+      ENDIF
+#endif
 
       ! 2. Vertical diffusion on u
@@ -119 +163 @@
          DO jj = 2, jpjm1 
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zcoef = - p2dt / fse3u(ji,jj,jk)
+               ze3ua =  ( 1. - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl   * fse3u_a(ji,jj,jk)   ! after scale factor at T-point
+               zcoef = - p2dt / ze3ua      
                zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
                zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
                zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
                zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
-               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
-            END DO
-         END DO
-      END DO
-      DO jj = 2, jpjm1        ! Surface boudary conditions
+               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
+            END DO
+         END DO
+      END DO
+      DO jj = 2, jpjm1        ! Surface boundary conditions
          DO ji = fs_2, fs_jpim1   ! vector opt.
             zwi(ji,jj,1) = 0._wp
-            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
+            zwd(ji,jj,1) = 1 - zws(ji,jj,1)
          END DO
       END DO
@@ -160 +205 @@
       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * (  ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
-               &                                                       * r1_rau0 / fse3u(ji,jj,1)       )
+            ze3ua =  ( 1. - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1) 
+#if defined key_dynspg_ts
+            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
+               &                                      / ( ze3ua * rau0 ) 
+#else
+            ua(ji,jj,1) = ub(ji,jj,1) + p2dt *(ua(ji,jj,1) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
+               &                                                     / ( fse3u(ji,jj,1) * rau0     ) ) 
+#endif
          END DO
       END DO
@@ -167 +218 @@
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side
+#if defined key_dynspg_ts
+               zrhs = ua(ji,jj,jk)   ! zrhs=right hand side
+#else
+               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)
+#endif
                ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
             END DO
@@ -186 +241 @@
       END DO
 
+#if ! defined key_dynspg_ts
       ! Normalization to obtain the general momentum trend ua
       DO jk = 1, jpkm1
@@ -194 +250 @@
          END DO
       END DO
-
+#endif
 
       ! 3. Vertical diffusion on v
@@ -205 +261 @@
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zcoef = -p2dt / fse3v(ji,jj,jk)
+               ze3va =  ( 1. - r_vvl ) * fse3v_n(ji,jj,jk)  + r_vvl * fse3v_a(ji,jj,jk)   ! after scale factor at T-point
+               zcoef = - p2dt / ze3va
                zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
                zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk)
                zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
                zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
-               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
-            END DO
-         END DO
-      END DO
-      DO jj = 2, jpjm1        ! Surface boudary conditions
+               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
+            END DO
+         END DO
+      END DO
+      DO jj = 2, jpjm1        ! Surface boundary conditions
          DO ji = fs_2, fs_jpim1   ! vector opt.
             zwi(ji,jj,1) = 0._wp
-            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
+            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
          END DO
       END DO
@@ -246 +303 @@
       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            va(ji,jj,1) = vb(ji,jj,1) + p2dt * (  va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
-               &                                                       * r1_rau0 / fse3v(ji,jj,1)       )
+            ze3va =  ( 1. - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1) 
+#if defined key_dynspg_ts            
+            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
+               &                                      / ( ze3va * rau0 ) 
+#else
+            va(ji,jj,1) = vb(ji,jj,1) + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
+               &                                                       / ( fse3v(ji,jj,1) * rau0     )  )
+#endif
          END DO
       END DO
@@ -253 +316 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side
+#if defined key_dynspg_ts
+               zrhs = va(ji,jj,jk)   ! zrhs=right hand side
+#else
+               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)
+#endif
                va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
             END DO
@@ -273 +340 @@
 
       ! Normalization to obtain the general momentum trend va
+#if ! defined key_dynspg_ts
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1   
@@ -280 +348 @@
          END DO
       END DO
-
+#endif
+
+      ! J. Chanut: Lines below are useless ?
       !! restore bottom layer avmu(v) 
       IF( ln_bfrimp ) THEN
@@ -292 +362 @@
             ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
             ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
-            avmu(ji,jj,ikbu+1) = zavmu(ji,jj)
-            avmv(ji,jj,ikbv+1) = zavmv(ji,jj)
+            avmu(ji,jj,ikbu+1) = 0.e0
+            avmv(ji,jj,ikbv+1) = 0.e0
          END DO
       END DO
@@ -299 +369 @@
       !
       CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws) 
-      CALL wrk_dealloc( jpi,jpj, zavmu, zavmv) 
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp')

branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

@@ -8 +8 @@
    !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
    !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
+   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
    !!----------------------------------------------------------------------
 
@@ -28 +29 @@
    USE obc_oce
    USE bdy_oce
+   USE bdy_par         
+   USE bdydyn2d        ! bdy_ssh routine
    USE diaar5, ONLY:   lk_diaar5
    USE iom
-   USE sbcrnf, ONLY: h_rnf, nk_rnf   ! River runoff 
+   USE sbcrnf          ! River runoff 
 #if defined key_agrif
    USE agrif_opa_update
@@ -40 +43 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+
+#if defined key_dynspg_ts
+   ! jchanut: Needed modules for dynamics update:
+   USE eosbn2           ! equation of state                (eos_bn2 routine)
+   USE zpshde           ! partial step: hor. derivative    (zps_hde routine)
+   USE dynadv           ! advection                        (dyn_adv routine)
+   USE dynvor           ! vorticity term                   (dyn_vor routine)
+   USE dynhpg           ! hydrostatic pressure grad.       (dyn_hpg routine)
+   USE dynldf           ! lateral momentum diffusion       (dyn_ldf routine)
+   USE dynspg_oce       ! surface pressure gradient        (dyn_spg routine)
+   USE dynspg           ! surface pressure gradient        (dyn_spg routine)
+   USE dynnept          ! simp. form of Neptune effect(dyn_nept_cor routine)
+   USE asminc           ! assimilation increments      (dyn_asm_inc routine)
+   USE bdydyn3d         ! (bdy_dyn3d_dmp routine)  
+   USE dynspg_ts        ! for advective velocities issued from time splitting
+   USE zdfbfr           ! Bottom friction
+#if defined key_agrif
+   USE agrif_opa_sponge ! Momemtum and tracers sponges
+#endif
+#endif
 
    IMPLICIT NONE
@@ -79 +102 @@
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   indic
       REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars
       REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d, zhdiv
@@ -146 +170 @@
          hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
          hv(:,:) = hv_0(:,:) + sshv_n(:,:)
+         ! bg jchanut tschanges
+#if defined key_dynspg_ts
+         ht(:,:) = ht_0(:,:) + sshn(:,:)              ! now ocean depth (at t- and f-points)
+         hf(:,:) = hf_0(:,:) + sshf_n(:,:)
+#endif
+         ! end jchanut tschanges
          !                                            ! now masked inverse of the ocean depth (at u- and v-points)
          hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) )
@@ -180 +210 @@
 #endif
 #if defined key_bdy
-      ssha(:,:) = ssha(:,:) * bdytmask(:,:)
-      CALL lbc_lnk( ssha, 'T', 1. )                    ! absolutly compulsory !! (jmm)
-#endif
+      ! bg jchanut tschanges
+      ! These lines are not necessary with time splitting since
+      ! boundary condition on sea level is set during ts loop
+      IF (lk_bdy) THEN
+         CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary
+         CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries
+      ENDIF
+#endif
+      ! end jchanut tschanges
 #if defined key_asminc
       !                                                ! Include the IAU weighted SSH increment
@@ -190 +226 @@
       ENDIF
 #endif
-      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
-      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1      ! NO Vector Opt.
-               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
-                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
-                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
-               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
-                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
-                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
-            END DO
-         END DO
-         CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions
-      ENDIF
-
       !                                           !------------------------------!
       !                                           !     Now Vertical Velocity    !
@@ -219 +240 @@
       END DO
 
+      ! bg jchanut tschanges
+#if defined key_dynspg_ts
+      ! In case the time splitting case, update most of all momentum trends here:
+      ! Note that the computation of vertical velocity above, hence "after" sea level is necessary 
+      ! to compute momentum advection for the rhs of barotropic loop:
+                               CALL eos    ( tsn, rhd, rhop )                   ! now in situ density for hpg computation
+
+      IF( ln_zps      )        CALL zps_hde( kt, jpts, tsn, gtsu, gtsv,  &      ! zps: now hor. derivative
+            &                                          rhd, gru , grv  )        ! of t, s, rd at the last ocean level
+                               CALL zdf_bfr( kt )         ! bottom friction (if quadratic)
+
+                               ua(:,:,:) = 0.e0           ! set dynamics trends to zero
+                               va(:,:,:) = 0.e0
+      IF(  ln_asmiau .AND. &
+         & ln_dyninc       )   CALL dyn_asm_inc( kt )     ! apply dynamics assimilation increment
+      IF( ln_neptsimp )        CALL dyn_nept_cor( kt )    ! subtract Neptune velocities (simplified)
+      IF( lk_bdy           )   CALL bdy_dyn3d_dmp( kt )   ! bdy damping trends
+                               CALL dyn_adv( kt )         ! advection (vector or flux form)
+                               CALL dyn_vor( kt )         ! vorticity term including Coriolis
+                               CALL dyn_ldf( kt )         ! lateral mixing
+      IF( ln_neptsimp )        CALL dyn_nept_cor( kt )    ! add Neptune velocities (simplified)
+#if defined key_agrif
+      IF(.NOT. Agrif_Root())   CALL Agrif_Sponge_dyn      ! momentum sponge
+#endif
+                               CALL dyn_hpg( kt )         ! horizontal gradient of Hydrostatic pressure
+                               CALL dyn_spg( kt, indic )  ! surface pressure gradient
+
+                               ua_bak(:,:,:) = ua(:,:,:)  ! save next velocities (not trends !)
+                               va_bak(:,:,:) = va(:,:,:)
+
+      ! At this stage: 
+      ! 1) ssha, sshu_a, sshv_a have been updated. 
+      ! 2) Time averaged barotropic velocities around step kt+1 (ua_b, va_b) as well as 
+      !   "advective" barotropic velocities (un_adv, vn_adv) at step kt in agreement 
+      !   with continuity equation are available.
+      ! 3) 3d velocities (ua, va) hold ua_b, va_b issued from time splitting as barotropic components.
+      ! Below => Update "now" velocities, divergence, then vertical velocity
+      ! NB: hdivn is NOT updated such that hdivb is not updated also. This means that horizontal
+      !     momentum diffusion is still performed on Instantaneous barotropic arrays. I experencied
+      !     some issues with UBS with the current method. Uncomment "divup" below to update the divergence.
+      !
+      CALL wrk_alloc( jpi,jpj,jpk, z3d )
+      !
+      DO jk = 1, jpkm1 
+         ! Correct velocities:                            
+         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk)
+         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk)
+         !
+         ! Compute divergence:
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               z3d(ji,jj,jk) =   &
+                  (  e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk) * un(ji-1,jj,jk)       &
+                   + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk) * vn(ji,jj-1,jk)  )    &
+                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+            END DO  
+         END DO 
+      END DO
+      !
+      IF( ln_rnf      )   CALL sbc_rnf_div( z3d )      ! runoffs (update divergence)
+      !
+      !
+      ! Set new vertical velocities:
+      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
+         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
+         wn(:,:,jk) = wn(:,:,jk+1) -   fse3t_n(:,:,jk) * z3d(:,:,jk)        &    
+            &                      -  (fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )  &
+            &                          * tmask(:,:,jk) * z1_2dt
+#if defined key_bdy
+         ! JC: line below is purely cosmetic I guess:
+         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
+#endif
+      END DO
+      !
+      CALL wrk_dealloc( jpi,jpj,jpk, z3d )
+
+!! bg divup
+!!      !
+!!      DO jk = 1, jpkm1 
+!!         ! Correct velocities:                            
+!!         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk)
+!!         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk)
+!!         !
+!!      END DO
+!!      !
+!!      !
+!!      CALL div_cur( kt )                               ! Horizontal divergence & Relative vorticity
+!!      !
+!!      ! Set new vertical velocities:
+!!      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
+!!         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
+!!         wn(:,:,jk) = wn(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)      &  
+!!            &                      -  (fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )  &
+!!            &                          * tmask(:,:,jk) * z1_2dt
+!!#if defined key_bdy
+!!         ! JC: line below is purely cosmetic I guess:
+!!         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
+!!#endif
+!!      END DO
+!! end divup
+      !
+      !
+      ! End of time splitting case
+#else
+      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
+      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1      ! NO Vector Opt.
+               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
+                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
+                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
+               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
+                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
+                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
+            END DO
+         END DO
+         CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions
+      ENDIF
+#endif
       !                                           !------------------------------!
       !                                           !           outputs            !
@@ -283 +423 @@
          !                    !--------------------------!
          !
-         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
+#if defined key_dynspg_ts
+         IF(( neuler == 0 .AND. kt == nit000 ).OR.(lk_dynspg_ts.AND.ln_bt_fw)) THEN    !** Euler time-stepping: no filter
+#else
+         IF ( neuler == 0 .AND. kt == nit000 ) THEN
+#endif
+            sshb  (:,:) = sshn  (:,:)
             sshn  (:,:) = ssha  (:,:)                       ! now <-- after  (before already = now)
+            sshu_b(:,:) = sshu_n(:,:)
+            sshv_b(:,:) = sshv_n(:,:)
             sshu_n(:,:) = sshu_a(:,:)
             sshv_n(:,:) = sshv_a(:,:)
@@ -336 +483 @@
          !                    !--------------------------!
          !
-         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
+#if defined key_dynspg_ts
+         IF(( neuler == 0 .AND. kt == nit000 ).OR.(lk_dynspg_ts.AND.ln_bt_fw)) THEN    !** Euler time-stepping: no filter
+#else
+         IF( neuler == 0 .AND. kt == nit000 ) THEN
+#endif
+            sshb(:,:) = sshn(:,:)
             sshn(:,:) = ssha(:,:)                           ! now <-- after  (before already = now)
             !

branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90

@@ -4 +4 @@
   !! Initialization of tidal forcing
   !! History :  9.0  !  07  (O. Le Galloudec)  Original code
+  !!            3.5  !  2013-07   (J. Chanut) Compliant with time splitting changes
   !!=================================================================================
 #if defined key_tide
@@ -26 +27 @@
 CONTAINS
 
-  SUBROUTINE upd_tide (kt,kit)
+  SUBROUTINE upd_tide (kt, kit, time_offset)
     !!----------------------------------------------------------------------
     !!                 ***  ROUTINE upd_tide  ***
@@ -32 +33 @@
     !! * Local declarations
 
-    INTEGER, INTENT( in ) ::   kt,kit      ! ocean time-step index
-    INTEGER  :: ji,jj,jk
-    REAL (wp) :: zramp
-    REAL (wp), DIMENSION(nb_harmo) :: zwt 
+    INTEGER, INTENT( in ) ::   kt            ! ocean time-step index
+    INTEGER, INTENT( in ), OPTIONAL  :: kit  ! Barotropic timestep counter (ts option)
+    INTEGER, INTENT( in ), OPTIONAL  :: time_offset ! time offset in units of timesteps
+    INTEGER                          :: time_add    ! time offset in units of timesteps
+    INTEGER  :: ji, jj, jk                   
+    REAL (wp) :: zramp, z_arg, z_t
     !...............................................................................
 
     pot_astro(:,:)=0.e0
     zramp = 1.e0
+    time_add = 0
 
-    IF (lk_dynspg_ts) THEN
-       zwt(:) = omega_tide(:)* ((kt-kt_tide)*rdt + kit*(rdt/REAL(nn_baro,wp)))
-       IF (ln_tide_ramp) THEN
-          zramp = MIN(MAX( ((kt-nit000)*rdt + kit*(rdt/REAL(nn_baro,wp)))/(rdttideramp*rday),0.),1.)
-       ENDIF
-    ELSE
-       zwt(:) = omega_tide(:)*(kt-kt_tide)*rdt
-       IF (ln_tide_ramp) THEN
-          zramp = MIN(MAX( ((kt-nit000)*rdt)/(rdttideramp*rday),0.),1.) 
-       ENDIF  
+    IF( PRESENT(time_offset) ) THEN
+       time_add = time_offset
+    ENDIF
+         
+    IF( PRESENT(kit) ) THEN  
+       z_arg = ((kt-kt_tide) * rdt + (kit+0.5*(time_add-1)) * rdt / REAL(nn_baro,wp) )
+    ELSE                              
+       z_arg = ((kt-kt_tide)+time_add) * rdt
     ENDIF
 
-    do jk=1,nb_harmo
-       do ji=1,jpi
-          do jj=1,jpj
-             pot_astro(ji,jj)=pot_astro(ji,jj) + zramp*(amp_pot(ji,jj,jk)*COS(zwt(jk)+phi_pot(ji,jj,jk)))      
-          enddo
-       enddo
-    enddo
+    IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0.),1.)
+
+    DO jk=1,nb_harmo
+       z_t = z_arg * omega_tide(jk)
+       DO ji=1,jpi
+          DO jj=1,jpj
+             pot_astro(ji,jj)=pot_astro(ji,jj) + zramp * amp_pot(ji,jj,jk)*COS(z_t + phi_pot(ji,jj,jk))      
+          END DO
+       END DO
+    END DO
 
   END SUBROUTINE upd_tide
@@ -68 +73 @@
   !!----------------------------------------------------------------------
 CONTAINS
-  SUBROUTINE upd_tide( kt,kit )          ! Empty routine
-    INTEGER,INTENT (IN) :: kt, kit
+  SUBROUTINE upd_tide( kt)          ! Empty routine
+    INTEGER,INTENT (IN) :: kt
     WRITE(*,*) 'upd_tide: You should not have seen this print! error?', kt
   END SUBROUTINE upd_tide
@@ -78 +83 @@
 
 END MODULE updtide
+

branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

@@ -343 +343 @@
       IF( lk_bdy        )   CALL     bdy_init       ! Open boundaries initialisation
       IF( lk_bdy        )   CALL     bdy_dta_init   ! Open boundaries initialisation of external data arrays
-      IF( lk_bdy        )   CALL     bdytide_init   ! Open boundaries initialisation of tidal harmonic forcing
+      IF( lk_bdy.AND.lk_tide )   CALL     bdytide_init   ! Open boundaries initialisation of tidal harmonic forcing
 
                             CALL dyn_nept_init  ! simplified form of Neptune effect

branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/oce.F90

@@ -22 +22 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   ub   ,  un    , ua     !: i-horizontal velocity        [m/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   vb   ,  vn    , va     !: j-horizontal velocity        [m/s]
+   ! bg jchanut tschanges
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   ua_bak   ,  va_bak     !: Saved trends for mod. ts     [m/s2]
+   ! end jchanut tschanges
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           wn             !: vertical velocity            [m/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   rotb ,  rotn           !: relative vorticity           [s-1]
@@ -70 +73 @@
       !
       ALLOCATE( ub   (jpi,jpj,jpk)      , un   (jpi,jpj,jpk)      , ua(jpi,jpj,jpk)       ,     &
-         &      vb   (jpi,jpj,jpk)      , vn   (jpi,jpj,jpk)      , va(jpi,jpj,jpk)       ,     &      
+         &      vb   (jpi,jpj,jpk)      , vn   (jpi,jpj,jpk)      , va(jpi,jpj,jpk)       ,     & 
+      ! bg jchanut tschanges
+#if defined key_dynspg_ts  
+      ! These temporary arrays are used to save tendencies computed before the time stepping of tracers.
+      ! These could be suppressed if ua and va would not have been used as temporary arrays 
+      ! during tracers' update
+         &      ua_bak(jpi,jpj,jpk)     , va_bak(jpi,jpj,jpk)     ,                             &
+#endif
+      ! end jchanut tschanges
          &      wn   (jpi,jpj,jpk)      ,                                                       &
          &      rotb (jpi,jpj,jpk)      , rotn (jpi,jpj,jpk)      ,                             &   

branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/step.F90

@@ -108 +108 @@
       !
       !  VERTICAL PHYSICS
+      ! bg jchanut tschanges
+      ! One need bottom friction parameter in ssh_wzv routine with time splitting.
+      ! The idea could be to move the call below before ssh_wzv. However, "now" scale factors
+      ! at U-V points (which are set thanks to sshu_n, sshv_n) are actually available in sshwzv.
+      ! These are needed for log bottom friction...
+#if ! defined key_dynspg_ts
                          CALL zdf_bfr( kstp )         ! bottom friction
+#endif
+      ! end jchanut tschanges
 
       !                                               ! Vertical eddy viscosity and diffusivity coefficients
@@ -206 +214 @@
             &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level
 
-      ELSE                                                  ! centered hpg  (eos then time stepping)
+      ELSE   
+                                               ! centered hpg  (eos then time stepping)
+      ! bg jchanut tschanges
+#if ! defined key_dynspg_ts
+      ! eos already called
                              CALL eos    ( tsn, rhd, rhop )      ! now in situ density for hpg computation
          IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &    ! zps: now hor. derivative
             &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level
+#endif
+      ! end jchanut tschanges
          IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
                              CALL tra_nxt( kstp )                ! tracer fields at next time step
@@ -217 +231 @@
       ! Dynamics                                    (tsa used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      ! bg jchanut tschanges
+#if defined key_dynspg_ts      
+! revert to previously computed tendencies:
+! (not using ua, va as temporary arrays during tracers' update could avoid that)
+                               ua(:,:,:) = ua_bak(:,:,:)            
+                               va(:,:,:) = va_bak(:,:,:)
+                               CALL dyn_bfr( kstp )         ! bottom friction
+                               CALL dyn_zdf( kstp )         ! vertical diffusion
+#else
+      ! end jchanut tschanges
                                ua(:,:,:) = 0.e0             ! set dynamics trends to zero
                                va(:,:,:) = 0.e0
@@ -236 +260 @@
                                CALL dyn_zdf( kstp )         ! vertical diffusion
                                CALL dyn_spg( kstp, indic )  ! surface pressure gradient
+      ! bg jchanut tschanges
+#endif
+      ! end jchanut tschanges
                                CALL dyn_nxt( kstp )         ! lateral velocity at next time step