Index: trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 1656)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 1662)
@@ -22,4 +22,5 @@
    USE phycst          ! physical constants
    USE domvvl          ! variable volume
+   USE zdfbfr          ! bottom friction
    USE obcdta          ! open boundary condition data     
    USE obcfla          ! Flather open boundary condition  
@@ -91,11 +92,15 @@
       INTEGER, INTENT(in)  ::   kt   ! ocean time-step index
       !!
-      INTEGER  ::   ji, jj, jk, jn             ! dummy loop indices
-      INTEGER  ::   icycle                      ! temporary scalar
-      REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b,   &  ! temporary scalars
-         z1_8, zspgu, zcubt, zx1, zy1,    &  !     "        "
-         z1_4, zspgv, zcvbt, zx2, zy2        !     "        "
+      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
+      INTEGER  ::   icycle           ! temporary scalar
+
+      REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b     ! temporary 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), DIMENSION(jpi,jpj) ::   zhdiv, zsshb_e
-      !! 
+      !!
       REAL(wp), DIMENSION(jpi,jpj) ::   zsshun_e, zsshvn_e   ! 2D workspace
       !!
@@ -144,4 +149,6 @@
       !
       zhdiv(:,:) = 0.e0                                     ! barotropic divergence
+      zu_sld = 0.e0   ;   zu_asp = 0.e0                     ! tides trends (lk_tide=F)
+      zv_sld = 0.e0   ;   zv_asp = 0.e0
 
       ! -----------------------------------------------------------------------------
@@ -245,4 +252,22 @@
       END DO
 
+      !                                            ! Remove barotropic contribution of bottom friction 
+                                                   ! from the barotropic transport trend
+      IF( lk_vvl ) THEN
+       DO jj = 2, jpjm1
+          DO ji = fs_2, fs_jpim1   ! vector opt.
+             zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * zub(ji,jj) / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) )
+             zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * zvb(ji,jj) / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - 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) - bfrua(ji,jj) * zub(ji,jj) * hur(ji,jj)
+             zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * zvb(ji,jj) * hvr(ji,jj)
+          END DO
+       END DO
+      ENDIF
+
       !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity)
       !                                   ! -------------------------- 
@@ -272,5 +297,6 @@
       hur_e  (:,:) = hur  (:,:)            ! ocean depth inverted at u- and v-points
       hvr_e  (:,:) = hvr  (:,:)
-      zsshb_e(:,:) = sshn (:,:)            ! sea surface height (before and now)
+!RBbug     zsshb_e(:,:) = sshn (:,:)  
+      zsshb_e(:,:) = sshn_b(:,:)           ! sea surface height (before and now)
       sshn_e (:,:) = sshn (:,:)
       
@@ -343,11 +369,11 @@
                   ! surface pressure gradient
                   IF( lk_vvl) THEN
-                     zspgu = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
-                        &             - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
-                     zspgv = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
-                        &             - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
+                     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
-                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
-                     zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
+                     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
                   ! energy conserving formulation for planetary vorticity term
@@ -356,9 +382,11 @@
                   zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
                   zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
-                  zcubt = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj)
-                  zcvbt =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj)
-                  ! after barotropic velocity
-                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
-                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
+                  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
+                  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.e0         - 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.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
                END DO
             END DO
@@ -369,20 +397,22 @@
                    ! surface pressure gradient
                   IF( lk_vvl) THEN
-                     zspgu = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
-                        &             - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
-                     zspgv = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
-                        &             - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
+                     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
-                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
-                     zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
+                     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
                   ! 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)
-                  zcubt  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj)
-                  zcvbt  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj)
-                  ! after barotropic velocity
-                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
-                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
+                  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
+                  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.e0         - 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.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
                END DO
             END DO
@@ -393,20 +423,22 @@
                   ! surface pressure gradient
                   IF( lk_vvl) THEN
-                     zspgu = -grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  )    &
-                        &             - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
-                     zspgv = -grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
-                        &             - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
+                     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
-                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
-                     zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
+                     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
                   ! energy/enstrophy conserving formulation for planetary vorticity term
-                  zcubt = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
+                  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)
-                  zcvbt = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
+                  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 barotropic velocity
-                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
-                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
+                  ! after velocities with implicit bottom friction
+                  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.e0         - 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.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
                END DO
             END DO
@@ -499,4 +531,5 @@
       IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' )
       !
+      !
    END SUBROUTINE dyn_spg_ts
 
Index: trunk/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 1656)
+++ trunk/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 1662)
@@ -73,4 +73,5 @@
       INTEGER ::   ji, jj, jk                          ! dummy loop indices
       REAL(wp) ::   zrau0r, z2dtf, zcoef, zzws, zrhs   ! temporary scalars
+      REAL(wp) ::   zzwi                               ! temporary scalars
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi        ! temporary workspace arrays
       !!----------------------------------------------------------------------
@@ -89,12 +90,15 @@
       ! ---------------------------
       ! Matrix and second member construction
-      ! bottom boundary condition: only zws must be masked as avmu can take
+      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
       ! non zero value at the ocean bottom depending on the bottom friction
-      ! used (see zdfmix.F)
+      ! used but the bottom velocities have already been updated with the bottom
+      ! friction velocity in dyn_bfr using values from the previous timestep. There
+      ! is no need to include these in the implicit calculation.
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1 
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zcoef = - p2dt / fse3u(ji,jj,jk)
-               zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  )
+               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)
@@ -183,12 +187,15 @@
       ! ---------------------------
       ! Matrix and second member construction
-      ! bottom boundary condition: only zws must be masked as avmv can take
+      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
       ! non zero value at the ocean bottom depending on the bottom friction
-      ! used (see zdfmix.F)
+      ! used but the bottom velocities have already been updated with the bottom
+      ! friction velocity in dyn_bfr using values from the previous timestep. There
+      ! is no need to include these in the implicit calculation.
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zcoef = -p2dt / fse3v(ji,jj,jk)
-               zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  )
+               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)
Index: trunk/NEMO/OPA_SRC/TRD/trdmod.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdmod.F90	(revision 1656)
+++ trunk/NEMO/OPA_SRC/TRD/trdmod.F90	(revision 1662)
@@ -127,13 +127,5 @@
                      ztswu(ji,jj) = utau(ji,jj) / ( fse3u(ji,jj,1)*rau0 )
                      ztswv(ji,jj) = vtau(ji,jj) / ( fse3v(ji,jj,1)*rau0 )
-                     ! save bottom friction momentum fluxes
-                     ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) )
-                     ikbv   = MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) )
-                     ikbum1 = MAX( ikbu-1, 1 )
-                     ikbvm1 = MAX( ikbv-1, 1 )
-                     zua = ua(ji,jj,ikbum1) * r2dt + ub(ji,jj,ikbum1)
-                     zva = va(ji,jj,ikbvm1) * r2dt + vb(ji,jj,ikbvm1)
-                     ztbfu(ji,jj) = - avmu(ji,jj,ikbu) * zua / ( fse3u(ji,jj,ikbum1)*fse3uw(ji,jj,ikbu) )
-                     ztbfv(ji,jj) = - avmv(ji,jj,ikbv) * zva / ( fse3v(ji,jj,ikbvm1)*fse3vw(ji,jj,ikbv) )
+                     ! bottom friction contribution now handled explicitly
                      !
                      ptrdx(ji,jj,1     ) = ptrdx(ji,jj,1     ) - ztswu(ji,jj)
@@ -146,5 +138,6 @@
                CALL trd_icp( ptrdx, ptrdy, jpicpd_zdf, ctype )   
                CALL trd_icp( ztswu, ztswv, jpicpd_swf, ctype )                               ! wind stress forcing term
-               CALL trd_icp( ztbfu, ztbfv, jpicpd_bfr, ctype )                               ! bottom friction term
+               ! bottom friction contribution now handled explicitly
+            CASE ( jpdyn_trd_bfr )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_bfr, ctype )     ! bottom friction term
             END SELECT
             !
Index: trunk/NEMO/OPA_SRC/ZDF/zdfbfr.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 1656)
+++ trunk/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 1662)
@@ -4,6 +4,7 @@
    !! Ocean physics: Bottom friction
    !!======================================================================
-   !! History :  8.0  ! 1997-06  (G. Madec, A.-M. Treguier)  Original code
+   !! History :  OPA  ! 1997-06  (G. Madec, A.-M. Treguier)  Original code
    !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
+   !!            3.2  ! 2001-09  (A.C.Coward)  Correction to include barotropic contribution
    !!----------------------------------------------------------------------
 
@@ -11,4 +12,6 @@
    !!   zdf_bfr      : update momentum Kz at the ocean bottom due to the type of bottom friction chosen
    !!   zdf_bfr_init : read in namelist and control the bottom friction parameters.
+   !!   zdf_bfr_2d   : read in namelist and control the bottom friction
+   !!                  parameters.
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers variables
@@ -29,4 +32,10 @@
    REAL(wp) ::   rn_bfri2 = 1.0e-3_wp   ! bottom drag coefficient (non linear case)
    REAL(wp) ::   rn_bfeb2 = 2.5e-3_wp   ! background bottom turbulent kinetic energy  [m2/s2]
+   REAL(wp) ::   rn_bfrien = 30_wp      ! local factor to enhance coefficient bfri
+   REAL(wp), DIMENSION(jpi,jpj) :: &
+                 bfrcoef2d = 1.e-3_wp   ! 2D bottom drag coefficient
+   LOGICAL  ::   ln_bfr2d = .false.     ! logical switch for 2D enhancement
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & 
+      &          bfrua , bfrva          ! Bottom friction coefficients set in zdfbfr
 
    !! * Substitutions
@@ -47,27 +56,37 @@
       !!              Kz at the ocean bottom.
       !!
-      !! ** Method  :   Update the value of avmu and avmv at the ocean bottom 
-      !!       level following the chosen friction type (no-slip, free-slip, 
-      !!       linear, or quadratic)
+      !! ** Method  :   Calculate and store part of the momentum trend due    
+      !!              to bottom friction following the chosen friction type 
+      !!              (free-slip, linear, or quadratic). The component
+      !!              calculated here is multiplied by the bottom velocity in
+      !!              dyn_bfr to provide the trend term.
+      !!
+      !! ** Action  :   bfrua , bfrva   bottom friction coefficients
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
       !!
-      INTEGER  ::   ji, jj   ! dummy loop indexes
-      INTEGER  ::   ikbu, ikbv, ikbum1, ikbvm1   ! temporary integers
-      REAL(wp) ::   zvu, zuv, zecu, zecv         ! temporary scalars
-      !!----------------------------------------------------------------------
-
-      IF( kt == nit000 )   CALL zdf_bfr_init
-
-      !                               ! --------------------------------------
-      SELECT CASE (nn_bfr)            ! Compute avmu, avmv at the ocean bottom
-      !                               ! --------------------------------------
-      !
-      CASE( 0 )                            !==  no-slip boundary condition  ==!
+      INTEGER  ::   ji, jj         ! dummy loop indexes
+      INTEGER  ::   ikbu, ikbum1   ! temporary integers
+      INTEGER  ::   ikbv, ikbvm1   ! temporary integers
+      REAL(wp) ::   zvu, zuv, zecu, zecv   ! temporary scalars
+      !!----------------------------------------------------------------------
+
+
+      IF( kt == nit000 )   CALL zdf_bfr_init   ! initialisation
+
+      IF( nn_bfr == 2 ) THEN                 ! quadratic botton friction
+         ! Calculate and store the quadratic bottom friction terms bfrua and bfrva
+         ! where bfrUa = C_d*SQRT(u_bot^2 + v_bot^2 + e_b) {U=[u,v]}
+         ! from these the trend due to bottom friction:  -F_h/e3U  can be calculated
+         ! where -F_h/e3U_bot = bfrUa*Ub/e3U_bot {U=[u,v]}
+         !
 # if defined key_vectopt_loop
          DO jj = 1, 1
+!CDIR NOVERRCHK
             DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
 # else
+!CDIR NOVERRCHK
          DO jj = 2, jpjm1
+!CDIR NOVERRCHK
             DO ji = 2, jpim1
 # endif
@@ -76,74 +95,24 @@
                ikbum1 = MAX( ikbu-1, 1 )
                ikbvm1 = MAX( ikbv-1, 1 )
-               avmu(ji,jj,ikbu) = 2. * avmu(ji,jj,ikbum1)
-               avmv(ji,jj,ikbv) = 2. * avmv(ji,jj,ikbvm1)
+               !
+               zvu  = 0.25 * (  vn(ji,jj  ,ikbum1) + vn(ji+1,jj  ,ikbum1)     &
+                  &           + vn(ji,jj-1,ikbum1) + vn(ji+1,jj-1,ikbum1)  )
+               zuv  = 0.25 * (  un(ji,jj  ,ikbvm1) + un(ji-1,jj  ,ikbvm1)     &
+                  &           + un(ji,jj+1,ikbvm1) + un(ji-1,jj+1,ikbvm1)  )
+               !
+               zecu = SQRT(  un(ji,jj,ikbum1) * un(ji,jj,ikbum1) + zvu*zvu + rn_bfeb2  )
+               zecv = SQRT(  vn(ji,jj,ikbvm1) * vn(ji,jj,ikbvm1) + zuv*zuv + rn_bfeb2  )
+               !
+               bfrua(ji,jj) = - bfrcoef2d(ji,jj) * zecu 
+               bfrva(ji,jj) = - bfrcoef2d(ji,jj) * zecv
             END DO
          END DO
-
-      CASE( 1 )                            !==  linear botton friction  ==!
-# if defined key_vectopt_loop
-         DO jj = 1, 1
-            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
-# else
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-# endif
-               ikbu = MIN( mbathy(ji+1,jj), mbathy(ji,jj) )
-               ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )
-               avmu(ji,jj,ikbu) = rn_bfri1 * fse3uw(ji,jj,ikbu)
-               avmv(ji,jj,ikbv) = rn_bfri1 * fse3vw(ji,jj,ikbv)
-            END DO
-         END DO
-
-      CASE( 2 )                            !==  quadratic botton friction  ==!
-# if defined key_vectopt_loop
-         DO jj = 1, 1
-!CDIR NOVERRCHK
-            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
-# else
-!CDIR NOVERRCHK
-         DO jj = 2, jpjm1
-!CDIR NOVERRCHK
-            DO ji = 2, jpim1
-# endif
-               ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) )
-               ikbv   = MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) )
-               ikbum1 = MAX( ikbu-1, 1 )
-               ikbvm1 = MAX( ikbv-1, 1 )
-               
-               zvu  = 0.25 * (  vn(ji,jj  ,ikbum1) + vn(ji+1,jj  ,ikbum1)     &
-                              + vn(ji,jj-1,ikbum1) + vn(ji+1,jj-1,ikbum1)  )
-               
-               zuv  = 0.25 * (  un(ji,jj  ,ikbvm1) + un(ji-1,jj  ,ikbvm1)     &
-                              + un(ji,jj+1,ikbvm1) + un(ji-1,jj+1,ikbvm1)  )
-               
-               zecu = SQRT(  un(ji,jj,ikbum1) * un(ji,jj,ikbum1) + zvu*zvu + rn_bfeb2  )
-               zecv = SQRT(  vn(ji,jj,ikbvm1) * vn(ji,jj,ikbvm1) + zuv*zuv + rn_bfeb2  )
-               
-               avmu(ji,jj,ikbu) = rn_bfri2 * zecu * fse3uw(ji,jj,ikbu)
-               avmv(ji,jj,ikbv) = rn_bfri2 * zecv * fse3vw(ji,jj,ikbv)
-            END DO
-         END DO
-
-      CASE( 3 )                            !==  free-slip boundary condition  ==!
-# if defined key_vectopt_loop
-         DO jj = 1, 1
-            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
-# else
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-# endif
-               ikbu = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) )
-               ikbv = MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) )
-               avmu(ji,jj,ikbu) = 0.e0
-               avmv(ji,jj,ikbv) = 0.e0
-            END DO
-         END DO
-         !
-      END SELECT
-      CALL lbc_lnk( avmu, 'U', 1. )   ;    CALL lbc_lnk( avmv, 'V', 1. )   ! Lateral boundary condition   (unchanged sign)
-
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=avmu, clinfo1=' bfr  - u: ', mask1=umask,             &
-         &                       tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask,ovlap=1, kdim=jpk )
+         !
+      ENDIF
+      !
+      CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition
+
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr  - u: ', mask1=umask,        &
+         &                       tab2d_2=bfrva, clinfo2=       ' v: ', mask2=vmask,ovlap=1 )
       !
    END SUBROUTINE zdf_bfr
@@ -157,38 +126,135 @@
       !!
       !! ** Method  :   Read the nammbf namelist and check their consistency
-      !!----------------------------------------------------------------------
-      NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfeb2
-      !!----------------------------------------------------------------------
-
-      REWIND( numnam )              ! Read Namelist nambfr : bottom momentum boundary condition
-      READ  ( numnam, nambfr )
-
-      IF(lwp) WRITE(numout,*)       ! Parameter print
+      !!              called at the first timestep (nit000)
+      !!----------------------------------------------------------------------
+      USE iom   ! I/O module for ehanced bottom friction file
+      !!
+      NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfeb2, ln_bfr2d, rn_bfrien
+      INTEGER ::   inum       ! logical unit for enhanced bottom friction file
+      INTEGER ::   ji, jj     ! dummy loop indexes
+      INTEGER ::   ikbu, ikbv ! temporary integers
+      INTEGER ::   ictu, ictv ! temporary integers
+      REAL(wp) ::  rminbfr, rmaxbfr
+                              ! temporary reals
+      !!----------------------------------------------------------------------
+
+      REWIND ( numnam )               !* Read Namelist nam_bfr : bottom momentum boundary condition
+      READ   ( numnam, nambfr )
+
+
+      !                               !* Parameter control and print
+      IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'zdf_bfr : momentum bottom friction'
       IF(lwp) WRITE(numout,*) '~~~~~~~'
-      IF(lwp) WRITE(numout,*) '   Namelist nambfr : set bottom friction parameters'
-
-      SELECT CASE (nn_bfr)          ! Parameter control
-      !
+      IF(lwp) WRITE(numout,*) '          Namelist nam_bfr : set bottom friction parameters'
+
+      SELECT CASE (nn_bfr)
+
       CASE( 0 )
-         IF(lwp) WRITE(numout,*) '      no-slip '
+         IF(lwp) WRITE(numout,*) '            free-slip '
+         !
+         bfrua(:,:) = 0.e0
+         bfrva(:,:) = 0.e0
          !
       CASE( 1 )
-         IF(lwp) WRITE(numout,*) '      linear botton friction      nn_bfr    = ', nn_bfr
-         IF(lwp) WRITE(numout,*) '      friction coef.              rn_bfri1  = ', rn_bfri1
+         IF(lwp) WRITE(numout,*) '            linear botton friction'
+         IF(lwp) WRITE(numout,*) '            friction coef.   rn_bfri1  = ', rn_bfri1
+         IF(lwp)  THEN
+            IF( ln_bfr2d ) THEN
+                 WRITE(numout,*) '            read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d
+                 WRITE(numout,*) '            coef rn_bfri2 enhancement factor   rn_bfrien  = ',rn_bfrien
+            ENDIF
+         ENDIF
+         !
+         bfrcoef2d(:,:) = rn_bfri1  ! initialize bfrcoef2d to the namelist variable
+
+         IF (ln_bfr2d) THEN 
+            ! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
+            CALL iom_open('bfr_coef.nc',inum)
+            CALL iom_get (inum, jpdom_data, 'bfr_coef',bfrcoef2d,1) ! bfrcoef2d is used as tmp array
+            CALL iom_close(inum)
+            bfrcoef2d(:,:)= rn_bfri1*(1+ rn_bfrien*bfrcoef2d(:,:))
+         ENDIF
+         bfrua(:,:) = - bfrcoef2d(:,:)
+         bfrva(:,:) = - bfrcoef2d(:,:)
          !
       CASE( 2 )
-         IF(lwp) WRITE(numout,*) '      quadratic botton friction   nn_bfr    = ', nn_bfr
-         IF(lwp) WRITE(numout,*) '      friction coef.              rn_bfri2  = ', rn_bfri2
-         IF(lwp) WRITE(numout,*) '      background KE               rn_bfeb2  = ', rn_bfeb2
-         !
-      CASE( 3 )
-         IF(lwp) WRITE(numout,*) '      free-slip '
+         IF(lwp) WRITE(numout,*) '            quadratic botton friction'
+         IF(lwp) WRITE(numout,*) '            friction coef.   rn_bfri2  = ', rn_bfri2
+         IF(lwp) WRITE(numout,*) '            background tke   rn_bfeb2  = ', rn_bfeb2
+         IF(lwp)  THEN
+            IF( ln_bfr2d ) THEN
+                 WRITE(numout,*) '            read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d
+                 WRITE(numout,*) '            coef rn_bfri2 enhancement factor   rn_bfrien  = ',rn_bfrien
+            ENDIF
+         ENDIF
+         bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable
+
+         IF (ln_bfr2d) THEN 
+            ! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
+            CALL iom_open('bfr_coef.nc',inum)
+            CALL iom_get (inum, jpdom_data, 'bfr_coef',bfrcoef2d,1) ! bfrcoef2d is used as tmp array
+            CALL iom_close(inum)
+            bfrcoef2d(:,:)= rn_bfri2*(1+ rn_bfrien*bfrcoef2d(:,:))
+         ENDIF
+
          !
       CASE DEFAULT
-         WRITE(ctmp1,*) 'bad flag value for nn_bfr = ', nn_bfr
+         WRITE(ctmp1,*) '         bad flag value for nn_bfr = ', nn_bfr
          CALL ctl_stop( ctmp1 )
          !
       END SELECT
+      !
+      ! Basic stability check on bottom friction coefficient
+      !
+      ictu = 0               ! counter for stability criterion breaches at U-pts
+      ictv = 0               ! counter for stability criterion breaches at V-pts
+      rminbfr =  1.e10_wp    ! initialise tracker for minimum of bottom friction coefficient
+      rmaxbfr = -1.e10_wp    ! initialise tracker for maximum of bottom friction coefficient
+      !
+#  if defined key_vectopt_loop
+      DO jj = 1, 1
+!CDIR NOVERRCHK
+         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
+#  else
+!CDIR NOVERRCHK
+      DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+         DO ji = 2, jpim1
+#  endif
+             ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) )
+             ikbv   = MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) )
+             IF ( bfrcoef2d(ji,jj).gt.2.0*fse3u(ji,jj,ikbu)/rdt ) then
+                IF ( ln_ctl ) THEN
+                   write(numout,*) 'BFR ',narea,nimpp+ji,njmpp+jj,ikbu
+                   write(numout,*) 'BFR ',bfrcoef2d(ji,jj),2.0*fse3u(ji,jj,ikbu)/rdt
+                ENDIF
+                ictu = ictu + 1
+             ENDIF
+             IF ( bfrcoef2d(ji,jj).gt.2.0*fse3v(ji,jj,ikbv)/rdt ) then
+                 IF ( ln_ctl ) THEN
+                     write(numout,*) 'BFR ',narea,nimpp+ji,njmpp+jj,ikbv
+                     write(numout,*) 'BFR ',bfrcoef2d(ji,jj),2.0*fse3v(ji,jj,ikbv)/rdt
+                 ENDIF
+                 ictv = ictv + 1
+             ENDIF
+             bfrcoef2d(ji,jj) = MIN(2.0*fse3u(ji,jj,ikbu)/rdt, bfrcoef2d(ji,jj))
+             bfrcoef2d(ji,jj) = MIN(2.0*fse3v(ji,jj,ikbv)/rdt, bfrcoef2d(ji,jj))
+             rminbfr = MIN(rminbfr, bfrcoef2d(ji,jj))
+             rmaxbfr = MAX(rmaxbfr, bfrcoef2d(ji,jj))
+         END DO
+      END DO
+      IF ( lk_mpp ) THEN
+         CALL mpp_sum( ictu )
+         CALL mpp_sum( ictv )
+         CALL mpp_min( rminbfr )
+         CALL mpp_max( rmaxbfr )
+      ENDIF
+      IF ( lwp .AND. ictu + ictv .GT. 0 ) THEN
+         WRITE(numout,*) ' Bottom friction stability check failed at ', ictu, ' U-points '
+         WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points '
+         WRITE(numout,*) ' Bottom friction coefficient reset where necessary'
+         WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', rminbfr, ' to ', rmaxbfr
+      ENDIF
       !
    END SUBROUTINE zdf_bfr_init
Index: trunk/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 1656)
+++ trunk/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 1662)
@@ -48,4 +48,5 @@
    USE in_out_manager ! I/O manager
    USE iom            ! I/O manager library
+   USE zdfbfr          ! bottom friction
 
    IMPLICIT NONE
@@ -182,4 +183,7 @@
       REAL(wp) ::   zus   , zwlc  , zind      !    -         -
       REAL(wp) ::   zzd_up, zzd_lw            !    -         -
+      INTEGER  ::   ikbu, ikbv, ikbum1, ikbvm1      ! temporary scalar
+      INTEGER  ::   ikbt, ikbumm1, ikbvmm1          ! temporary scalar
+      REAL(wp) ::   zebot                           ! temporary scalars
       INTEGER , DIMENSION(jpi,jpj)     ::   imlc    ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc    !  -      -
@@ -209,5 +213,26 @@
       !                     !  Bottom boundary condition on tke
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-!!gm to be added soon
+      !                     en(bot)   = 0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
+!CDIR NOVERRCHK
+      DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            ikbu = MIN( mbathy(ji+1,jj), mbathy(ji,jj) )
+            ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )
+            ikbum1 = MAX( ikbu-1, 1 )
+            ikbvm1 = MAX( ikbv-1, 1 )
+            ikbu = MIN( mbathy(ji,jj), mbathy(ji-1,jj) )
+            ikbv = MIN( mbathy(ji,jj), mbathy(ji,jj-1) )
+            ikbumm1 = MAX( ikbu-1, 1 )
+            ikbvmm1 = MAX( ikbv-1, 1 )
+            ikbt = MAX( mbathy(ji,jj), 1 )
+            ztx2 = bfrua(ji-1,jj) * fse3u(ji-1,jj  ,ikbumm1) + &
+                   bfrua(ji  ,jj) * fse3u(ji  ,jj  ,ikbum1 )
+            zty2 = bfrva(ji,jj  ) * fse3v(ji  ,jj  ,ikbvm1) + &
+                   bfrva(ji,jj-1) * fse3v(ji  ,jj-1,ikbvmm1 )
+            zebot = 0.25_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 )
+            en (ji,jj,ikbt) = MAX( zebot, rn_emin ) * tmask(ji,jj,1)
+         END DO
+      END DO
       !
       !
Index: trunk/NEMO/OPA_SRC/step.F90
===================================================================
--- trunk/NEMO/OPA_SRC/step.F90	(revision 1656)
+++ trunk/NEMO/OPA_SRC/step.F90	(revision 1662)
@@ -64,4 +64,5 @@
 
    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)
@@ -197,4 +198,5 @@
       !
       !  VERTICAL PHYSICS   
+                         CALL zdf_bfr( kstp )         ! bottom friction
       !                                               ! Vertical eddy viscosity and diffusivity coefficients
       IF( lk_zdfric     )   CALL zdf_ric    ( kstp )  ! Richardson number dependent Kz
@@ -216,7 +218,4 @@
       IF( lk_zdfddm .AND. .NOT. lk_zdfkpp )   &
          &               CALL zdf_ddm( kstp )         ! double diffusive mixing
-
-                         CALL zdf_bfr( kstp )         ! bottom friction
-
                          CALL zdf_mxl( kstp )         ! mixed layer depth
 
@@ -306,4 +305,5 @@
 #endif
                                CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure
+                               CALL dyn_bfr( kstp )           ! bottom friction   
                                CALL dyn_zdf( kstp )         ! vertical diffusion
                                indic=0
