Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90	(revision 9249)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90	(revision 9250)
@@ -62,5 +62,5 @@
          ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
          ztrdu(:,:,:) = ua(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:)
+         ztrdv(:,:,:) = va(:,:,:)
       ENDIF
 
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 9249)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 9250)
@@ -91,5 +91,9 @@
       ELSEIF(               kt <= nit000 + 1 ) THEN   ;   r2dt = 2. * rdt   ! = 2 rdt (leapfrog)
       ENDIF
-
+      !
+      !                             !* explicit top/bottom drag case
+      IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, ub, vb, ua, va )  ! add top/bottom friction trend to (ua,va)
+      !
+      !
       IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
          ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfdrg.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfdrg.F90	(revision 9249)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfdrg.F90	(revision 9250)
@@ -16,4 +16,5 @@
    !!----------------------------------------------------------------------
    !!   zdf_drg       : update bottom friction coefficient (non-linear bottom friction only)
+   !!   zdf_drg_exp   : compute the top & bottom friction in explicit case
    !!   zdf_drg_init  : read in namdrg namelist and control the bottom friction parameters.
    !!       drg_init  :
@@ -23,4 +24,6 @@
    USE dom_oce        ! ocean space and time domain variables
    USE zdf_oce        ! ocean vertical physics variables
+   USE trd_oce        ! trends: ocean variables
+   USE trddyn         ! trend manager: dynamics
    !
    USE in_out_manager ! I/O manager
@@ -34,4 +37,5 @@
 
    PUBLIC   zdf_drg         ! called by zdf_phy
+   PUBLIC   zdf_drg_exp     ! called by dyn_zdf
    PUBLIC   zdf_drg_init    ! called by zdf_phy_init
 
@@ -140,4 +144,79 @@
 
 
+   SUBROUTINE zdf_drg_exp( kt, pub, pvb, pua, pva )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_drg_exp  ***
+      !!
+      !! ** Purpose :   compute and add the explicit top and bottom frictions.
+      !!
+      !! ** Method  :   in explicit case, 
+      !!
+      !!              NB: in implicit case the calculation is performed in dynzdf.F90
+      !!
+      !! ** Action  :   (pua,pva)   momentum trend increased by top & bottom friction trend
+      !!---------------------------------------------------------------------
+      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pub, pvb   ! the two components of the before velocity
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! the two components of the velocity tendency
+      !! 
+      INTEGER  ::   ji, jj       ! dummy loop indexes
+      INTEGER  ::   ikbu, ikbv   ! local integers
+      REAL(wp) ::   zm1_2dt      ! local scalar
+      REAL(wp) ::   zCdu, zCdv   !   -      -
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv
+      !!---------------------------------------------------------------------
+      !
+!!gm bug : time step is only rdt (not 2 rdt if euler start !)
+      zm1_2dt = - 1._wp / ( 2._wp * rdt )
+
+      IF( l_trddyn ) THEN      ! trends: store the input trends
+         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
+         ztrdu(:,:,:) = pua(:,:,:)
+         ztrdv(:,:,:) = pva(:,:,:)
+      ENDIF
+
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            ikbu = mbku(ji,jj)          ! deepest wet ocean u- & v-levels
+            ikbv = mbkv(ji,jj)
+            !
+            ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt)
+            zCdu = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / e3u_n(ji,jj,ikbu)
+            zCdv = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / e3v_n(ji,jj,ikbv)
+            !
+            pua(ji,jj,ikbu) = pua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * pub(ji,jj,ikbu)
+            pva(ji,jj,ikbv) = pva(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * pvb(ji,jj,ikbv)
+         END DO
+      END DO
+      !
+      IF( ln_isfcav ) THEN        ! ocean cavities
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               ikbu = miku(ji,jj)          ! first wet ocean u- & v-levels
+               ikbv = mikv(ji,jj)
+               !
+               ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt)
+               zCdu = 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / e3u_n(ji,jj,ikbu)    ! NB: Cdtop masked
+               zCdv = 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / e3v_n(ji,jj,ikbv)
+               !
+               pua(ji,jj,ikbu) = pua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * pub(ji,jj,ikbu)
+               pva(ji,jj,ikbv) = pva(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * pvb(ji,jj,ikbv)
+           END DO
+         END DO
+      ENDIF
+      !
+      IF( l_trddyn ) THEN      ! trends: send trends to trddyn for further diagnostics
+         ztrdu(:,:,:) = pua(:,:,:) - ztrdu(:,:,:)
+         ztrdv(:,:,:) = pva(:,:,:) - ztrdv(:,:,:)
+         CALL trd_dyn( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_bfr, kt )
+         DEALLOCATE( ztrdu, ztrdv )
+      ENDIF
+      !                                          ! print mean trends (used for debugging)
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pua, clinfo1=' bfr  - Ua: ', mask1=umask,               &
+         &                       tab3d_2=pva, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      !
+   END SUBROUTINE zdf_drg_exp
+
+
    SUBROUTINE zdf_drg_init
       !!----------------------------------------------------------------------
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 9249)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 9250)
@@ -199,6 +199,4 @@
                             CALL wzv        ( kstp )              ! now cross-level velocity 
       ENDIF
-
-      IF( .NOT.ln_drgimp)   CALL dyn_bfr    ( kstp )  ! bottom friction
       
                          CALL dyn_zdf       ( kstp )  ! vertical diffusion
@@ -276,5 +274,5 @@
 !!jc2: dynnxt must be the latest call. e3t_b are indeed updated in that routine
                          CALL tra_nxt       ( kstp )  ! finalize (bcs) tracer fields at next time step and swap
-                         CALL dyn_nxt       ( kstp )  ! finalize (bcs) velocities at next time step and swap
+                         CALL dyn_nxt       ( kstp )  ! finalize (bcs) velocities at next time step and swap (always called after tra_nxt)
                          CALL ssh_swp       ( kstp )  ! swap of sea surface height
       IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 9249)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 9250)
@@ -37,5 +37,4 @@
    USE divhor          ! horizontal divergence            (div_hor routine)
    USE dynadv          ! advection                        (dyn_adv routine)
-   USE dynbfr          ! Bottom friction terms            (dyn_bfr routine)
    USE dynvor          ! vorticity term                   (dyn_vor routine)
    USE dynhpg          ! hydrostatic pressure grad.       (dyn_hpg routine)
