MODULE dynbfr
   !!==============================================================================
   !!                 ***  MODULE  dynbfr  ***
   !! Ocean dynamics :  bottom friction component of the momentum mixing trend
   !!==============================================================================
   !! History :  3.2  ! 2008-11  (A. C. Coward)  Original code
   !!            3.4  ! 2011-09  (H. Liu) Make it consistent with semi-implicit Bottom friction (ln_drgimp =T) 
   !!            4.0  ! 2017-05  (G. Madec)  drag coef. defined at t-point (zdfdrg.F90)
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   dyn_bfr       : Update the momentum trend with the bottom friction contribution
   !!----------------------------------------------------------------------
   USE oce            ! ocean dynamics and tracers variables
   USE dom_oce        ! ocean space and time domain variables 
   USE zdf_oce        ! vertical physics: variables
   USE zdfdrg         ! vertical physics: top/bottom drag coef.
   USE trd_oce        ! trends: ocean variables
   USE trddyn         ! trend manager: dynamics
   !
   USE in_out_manager ! I/O manager
   USE prtctl         ! Print control
   USE timing         ! Timing

   IMPLICIT NONE
   PRIVATE

   PUBLIC   dyn_bfr   !  routine called by step.F90

   !! * Substitutions
#  include "vectopt_loop_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
   !! $Id$
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS
   
   SUBROUTINE dyn_bfr( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE dyn_bfr  ***
      !!
      !! ** Purpose :   compute the bottom friction ocean dynamics physics.
      !!
      !!              only for explicit bottom friction form
      !!              implicit bfr is implemented in dynzdf_imp
      !!
      !! ** Action  :   (ua,va)   momentum trend increased by bottom friction trend
      !!---------------------------------------------------------------------
      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
      !! 
      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
      !!---------------------------------------------------------------------
      !
      IF( ln_timing )   CALL timing_start('dyn_bfr')
      !
!!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(:,:,:) = ua(:,:,:)
            ztrdv(:,:,:) = va(:,:,:)
      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)
            !
            ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * ub(ji,jj,ikbu)
            va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * vb(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)
               !
               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * ub(ji,jj,ikbu)
               va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * vb(ji,jj,ikbv)
           END DO
         END DO
      ENDIF
      !
      IF( l_trddyn ) THEN      ! trends: send trends to trddyn for further diagnostics
         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
         ztrdv(:,:,:) = va(:,:,:) - 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=ua, clinfo1=' bfr  - Ua: ', mask1=umask,               &
         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
      !
      IF( ln_timing )   CALL timing_stop('dyn_bfr')
      !
   END SUBROUTINE dyn_bfr

   !!==============================================================================
END MODULE dynbfr