diff --git a/src/OCE/DYN/dynspg_ts.F90 b/src/OCE/DYN/dynspg_ts.F90
index 69d04f57218a85c7a88020cfa1fd5135c39723ac..db3a442ad6872da49202845798d18b2d461d566c 100644
--- a/src/OCE/DYN/dynspg_ts.F90
+++ b/src/OCE/DYN/dynspg_ts.F90
@@ -911,7 +911,7 @@ CONTAINS
ENDIF
CALL iom_put( "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) ) ! barotropic i-current
- CALL iom_put( "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) ) ! barotropic i-current
+ CALL iom_put( "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) ) ! barotropic j-current
! !== END Phase 3 for MLF time integration ==!
#endif
diff --git a/src/OCE/stp2d.F90 b/src/OCE/stp2d.F90
index d865a9685fde16887fb647eec721c89bfeb432ba..b7a99bd0bf575369ecdfd1b0cd92ad894515b903 100644
--- a/src/OCE/stp2d.F90
+++ b/src/OCE/stp2d.F90
@@ -208,13 +208,13 @@ CONTAINS
! !== Net water flux forcing ==! (applied to a water column)
!
IF (ln_bt_fw) THEN ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2)
- sshe_rhs(:,:) = r1_rho0 * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) )
+ sshe_rhs(:,:) = r1_rho0 * ( emp(:,:) - rnf(:,:) - fwfisf_cav(:,:) - fwfisf_par(:,:) )
ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW)
zztmp = r1_rho0 * r1_2
sshe_rhs(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) &
& - rnf(:,:) - rnf_b(:,:) &
- & + fwfisf_cav(:,:) + fwfisf_cav_b(:,:) &
- & + fwfisf_par(:,:) + fwfisf_par_b(:,:) )
+ & - fwfisf_cav(:,:) - fwfisf_cav_b(:,:) &
+ & - fwfisf_par(:,:) - fwfisf_par_b(:,:) )
ENDIF
!
! !== Stokes drift divergence ==! (if exist)
diff --git a/src/OCE/stprk3.F90 b/src/OCE/stprk3.F90
index c027b91a4d840644c423f7ed6c80334119809c80..139cd57c62b73d552ad63ce9435d728f3ab1cb82 100644
--- a/src/OCE/stprk3.F90
+++ b/src/OCE/stprk3.F90
@@ -184,7 +184,7 @@ CONTAINS
! RK3 : single first external mode computation
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
- CALL stp_2D( kstp, Nbb, Nbb, Naa, Nrhs ) ! out: ssh, (uu_b,vv_b) and (un_adv,vn_adv) at Naa
+ CALL stp_2D( kstp, Nbb, Nbb, Naa, Nrhs ) ! out: ssh, (uu_b,vv_b) at Naa and (un_adv,vn_adv) between Nbb and Naa
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
diff --git a/src/OCE/stprk3_stg.F90 b/src/OCE/stprk3_stg.F90
index 9f9357cefa63bc18ef62b5208248efd88dd45b2a..38e90efd0f779c10d716a730531a7cb88b14963f 100644
--- a/src/OCE/stprk3_stg.F90
+++ b/src/OCE/stprk3_stg.F90
@@ -123,7 +123,7 @@ CONTAINS
!
! !== ssh/h0 ratio at Kaa ==!
!
- IF( .NOT.lk_linssh ) THEN ! "after" ssh/h_0 ratio at t,u,v-column
+ IF( .NOT.lk_linssh ) THEN ! "after" ssh/h_0 ratio at t,u,v-column computed at N+1 stored in r3.a
!
ALLOCATE( r3ta(jpi,jpj) , r3ua(jpi,jpj) , r3va(jpi,jpj) , r3fa(jpi,jpj) , r3fb(jpi,jpj) )
!
@@ -131,7 +131,7 @@ CONTAINS
CALL dom_qco_r3c_RK3( ssha, r3ta, r3ua, r3va, r3fa )
!
CALL lbc_lnk( 'stprk3_stg', r3ua, 'U', 1._wp, r3va, 'V', 1._wp, r3fa, 'F', 1._wp )
- !
+ ! !
r3t(:,:,Kaa) = r2_3 * r3t(:,:,Kbb) + r1_3 * r3ta(:,:) ! at N+1/3 (Kaa)
r3u(:,:,Kaa) = r2_3 * r3u(:,:,Kbb) + r1_3 * r3ua(:,:)
r3v(:,:,Kaa) = r2_3 * r3v(:,:,Kbb) + r1_3 * r3va(:,:)
@@ -194,7 +194,12 @@ CONTAINS
zaV(ji,jj,jk) = vv(ji,jj,jk,Kmm) + zvb(ji,jj)*vmask(ji,jj,jk)
END_3D
! !- vertical components -! ww
- CALL wzv ( kstp, Kbb, Kmm, Kaa, zaU, zaV, ww ) ! ww cross-level velocity
+ !
+ IF( ln_dynadv_vec ) THEN ! ww cross-level velocity consistent with uu/vv at Kmm
+ CALL wzv ( kstp, Kbb, Kmm, Kaa, uu(:,:,:,Kmm), vv(:,:,:,Kmm), ww )
+ ELSE ! ww cross-level velocity consistent with zaU/zaV
+ CALL wzv ( kstp, Kbb, Kmm, Kaa, zaU, zaV, ww )
+ ENDIF
!!st IF( ln_zad_Aimp ) CALL wAimp( kstp, ww, wi ) ! Adaptive-implicit vertical advection partitioning
! !==>>> implicite for stages 1 & 2 ????
@@ -209,7 +214,11 @@ CONTAINS
!
!===>>>>>> Modify dyn_adv_... dyn_keg routines so that Krhs to zero useless
! ! advection (VF or FF) ==> RHS
- CALL dyn_adv( kstp, Kbb, Kmm , uu, vv, Krhs, zaU, zaV, ww )
+ IF( ln_dynadv_vec ) THEN ! uu and vv used for momentum advection
+ CALL dyn_adv( kstp, Kbb, Kmm , uu, vv, Krhs)
+ ELSE ! advective velocity used for momentum advection
+ CALL dyn_adv( kstp, Kbb, Kmm , uu, vv, Krhs, zaU, zaV, ww )
+ ENDIF
! ! Coriolis / vorticity ==> RHS
CALL dyn_vor( kstp, Kmm , uu, vv, Krhs )
!
@@ -219,11 +228,7 @@ CONTAINS
!===>>>>>> Modify dyn_hpg & dyn_hpg_... routines : rhd computed in dyn_hpg and pass in argument to dyn_hpg_...
-!!st IF( kstg == 3 ) THEN
-! CALL eos ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0 ) ! now in situ density for hpg computation
-! ELSE
- CALL eos ( ts, Kmm, rhd ) ! Kmm in situ density anomaly for hpg computation
-! ENDIF
+ CALL eos ( ts, Kmm, rhd ) ! Kmm in situ density anomaly for hpg computation
!!gm end
CALL dyn_hpg( kstp, Kmm , uu, vv, Krhs )
@@ -237,13 +242,15 @@ CONTAINS
! vv(ji,jj,jk,Krhs) = vv(ji,jj,jk,Krhs) - grav * ( ssh(ji ,jj+1,Kmm) - ssh(ji,jj,Kmm) )
! END_3D
!!gm
-
+
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! RHS of tracers : ADV only using (zaU,zaV,ww)
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
-# if defined key_top
!
+ ! ! Advective velocity needed for tracers advection - already computed if ln_dynadv_vec=F
+ IF( ln_dynadv_vec ) CALL wzv ( kstp, Kbb, Kmm, Kaa, zaU, zaV, ww )
+ !
+# if defined key_top
! !== Passive Tracer ==!
!
SELECT CASE( kstg )
@@ -308,11 +315,7 @@ CONTAINS
!===>>>>>> Modify tra_adv_... routines so that Krhs to zero useless
DO jn = 1, jpts
-!!st does not work due to lbc_lnk north fold : need to be merged with the trunk LBC pb changes for the namelist
-!! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-!! ts(ji,jj,jk,jn,Krhs) = 0._wp ! set tracer trends to zero
-!! END_3D
- ts(:,:,:,jn,Krhs) = 0._wp
+ ts(:,:,:,jn,Krhs) = 0._wp ! set tracer trends to zero (:,:,:) needed otherwise it does not work (?)
END DO
CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0 ) ! now in potential density for tra_mle computation