Index: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90	(revision 3258)
+++ branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90	(revision 3259)
@@ -5,4 +5,6 @@
    !!==============================================================================
    !! 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_bfrimp = .true.) 
    !!----------------------------------------------------------------------
 
@@ -57,5 +59,4 @@
       IF( .NOT.ln_bfrimp) THEN     ! only for explicit bottom friction form
                                     ! implicit bfr is implemented in dynzdf_imp
-                                    ! H. Liu,  Sept. 2011
 
         zm1_2dt = - 1._wp / ( 2._wp * rdt )
Index: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 3258)
+++ branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 3259)
@@ -8,4 +8,5 @@
    !!   NEMO     0.5  !  2002-08  (G. Madec)  F90: Free form and module
    !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
+   !!            3.4  !  2012-01  (H. Liu) Semi-implicit bottom friction
    !!----------------------------------------------------------------------
 
@@ -20,5 +21,5 @@
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! MPP library
-   USE zdfbfr          ! bottom friction setup
+   USE zdfbfr          ! Bottom friction setup
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
@@ -60,11 +61,9 @@
       REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step
       !!
-      INTEGER  ::   ji, jj, jk     ! dummy loop indices
-      INTEGER  ::   ikbum1, ikbvm1 ! local variable
-      REAL(wp) ::   z1_p2dt, z2dtf, zcoef, zzwi, zzws, zrhs ! local scalars
-
-      !! * Local variables for implicit bottom friction.    H. Liu
-      REAL(wp) ::   zbfru, zbfrv 
-      REAL(wp) ::   zbfr_imp = 0._wp                        ! toggle (SAVE'd by assignment) 
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   ikbu, ikbv   ! local integers
+      REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars
+      !!----------------------------------------------------------------------
+
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws
       !!----------------------------------------------------------------------
@@ -78,5 +77,4 @@
          IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
-         IF(ln_bfrimp) zbfr_imp = 1._wp
       ENDIF
 
@@ -85,42 +83,50 @@
       z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
 
-      ! 1. Vertical diffusion on u
-
-      ! Vertical diffusion on u&v
+      ! 1. Apply semi-implicit bottom friction
+      ! --------------------------------------
+      ! Only needed for semi-implicit bottom friction setup. The explicit
+      ! bottom friction has been included in "u(v)a" which act as the R.H.S
+      ! column vector of the tri-diagonal matrix equation
+      !
+      IF( ln_bfrimp ) THEN
+# 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 = 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
+
+      ! 2. Vertical diffusion on u
       ! ---------------------------
       ! Matrix and second member construction
-      !! 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 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.
-
-      ! The code has been modified here to implicitly implement bottom
-      ! friction: u(v)mask is not necessary here anymore. 
-      ! H. Liu, April 2010.
-
-      ! 1. Vertical diffusion on u
-      DO jj = 2, jpjm1 
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            ikbum1 = mbku(ji,jj)
-               zbfru = bfrua(ji,jj)
-
-            DO jk = 1, ikbum1
+      ! 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.
+      !
+      DO jk = 1, jpkm1        ! Matrix
+         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  )
-               zws(ji,jj,jk) = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
-               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk)
-            END DO
-
-      ! Surface boundary conditions
+               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
+         DO ji = fs_2, fs_jpim1   ! vector opt.
             zwi(ji,jj,1) = 0._wp
             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
-
-      ! Bottom boundary conditions  ! H. Liu, May, 2010
-!           !commented out to be consistent with v3.2, h.liu
-!           z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 2.0_wp * zbfr_imp
-            z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 1.0_wp * zbfr_imp
-            zws(ji,jj,ikbum1) = 0._wp
-            zwd(ji,jj,ikbum1) = 1._wp - zwi(ji,jj,ikbum1) - z2dtf 
+         END DO
+      END DO
 
       ! Matrix inversion starting from the first level
@@ -138,23 +144,37 @@
       !   The solution (the after velocity) is in ua
       !-----------------------------------------------------------------------
-
-      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
-            DO jk = 2, ikbum1
+      !
+      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
+         DO jj = 2, jpjm1   
+            DO ji = fs_2, fs_jpim1   ! vector opt.
                zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
             END DO
-
-      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
-            z2dtf = 0.5_wp * p2dt / ( fse3u(ji,jj,1) * rau0 )
-            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ua(ji,jj,1) + z2dtf * (utau_b(ji,jj) + utau(ji,jj))
-           DO jk = 2, ikbum1   
+         END DO
+      END DO
+      !
+      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) )   &
+               &                                                       / ( fse3u(ji,jj,1) * rau0       )  )
+         END DO
+      END DO
+      DO jk = 2, jpkm1
+         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
                ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
             END DO
-
-
-      ! third recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk
-            ua(ji,jj,ikbum1) = ua(ji,jj,ikbum1) / zwd(ji,jj,ikbum1)
-            DO jk = ikbum1-1, 1, -1
-               ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
+         END DO
+      END DO
+      !
+      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
+         END DO
+      END DO
+      DO jk = jpk-2, 1, -1
+         DO jj = 2, jpjm1   
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
             END DO
          END DO
@@ -171,29 +191,28 @@
 
 
-      ! 2. Vertical diffusion on v
+      ! 3. Vertical diffusion on v
       ! ---------------------------
-
-      DO ji = fs_2, fs_jpim1   ! vector opt.
-         DO jj = 2, jpjm1   
-            ikbvm1 = mbkv(ji,jj)
-               zbfrv = bfrva(ji,jj)
-
-            DO jk = 1, ikbvm1
+      ! Matrix and second member construction
+      ! 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
+      !
+      DO jk = 1, jpkm1        ! Matrix
+         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  )
-               zws(ji,jj,jk) = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
-               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk)
-            END DO
-
-      ! Surface boundary conditions
+               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
+         DO ji = fs_2, fs_jpim1   ! vector opt.
             zwi(ji,jj,1) = 0._wp
             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
-
-      ! Bottom boundary conditions  ! H. Liu, May, 2010
-!           !commented out to be consistent with v3.2, h.liu
-!           z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 2.0_wp * zbfr_imp
-            z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 1.0_wp * zbfr_imp
-            zws(ji,jj,ikbvm1) = 0._wp
-            zwd(ji,jj,ikbvm1) = 1._wp - zwi(ji,jj,ikbvm1) - z2dtf 
+         END DO
+      END DO
 
       ! Matrix inversion
@@ -207,30 +226,42 @@
       !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
       !
-      !   m is decomposed in the product of an upper and lower triangular
-      !   matrix
+      !   m is decomposed in the product of an upper and lower triangular matrix
       !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
       !   The solution (after velocity) is in 2d array va
       !-----------------------------------------------------------------------
-
-      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
-            DO jk = 2, ikbvm1
+      !
+      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
+         DO jj = 2, jpjm1   
+            DO ji = fs_2, fs_jpim1   ! vector opt.
                zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
             END DO
-
-      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
-            z2dtf = 0.5_wp * p2dt / ( fse3v(ji,jj,1)*rau0 )
-            va(ji,jj,1) = vb(ji,jj,1) + p2dt * va(ji,jj,1) + z2dtf * (vtau_b(ji,jj) + vtau(ji,jj))
-            DO jk = 2, ikbvm1
+         END DO
+      END DO
+      !
+      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) )   &
+               &                                                       / ( fse3v(ji,jj,1) * rau0       )  )
+         END DO
+      END DO
+      DO jk = 2, jpkm1
+         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
                va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
             END DO
-
-      ! third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk
-            va(ji,jj,ikbvm1) = va(ji,jj,ikbvm1) / zwd(ji,jj,ikbvm1)
-
-            DO jk = ikbvm1-1, 1, -1
-               va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
-            END DO
-
+         END DO
+      END DO
+      !
+      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
+         END DO
+      END DO
+      DO jk = jpk-2, 1, -1
+         DO jj = 2, jpjm1   
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
+            END DO
          END DO
       END DO
Index: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 3258)
+++ branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 3259)
@@ -123,4 +123,5 @@
             &                       tab2d_2=bfrva, clinfo2=       ' v: ', mask2=vmask,ovlap=1 )
       ENDIF
+
       !
       IF( nn_timing == 1 )  CALL timing_stop('zdf_bfr')
@@ -170,5 +171,5 @@
             WRITE(numout,*) 'Implicit bottom friction can only be used when ln_zdfexp=.false.'
             WRITE(numout,*) '         but you set: ln_bfrimp=.true. and ln_zdfexp=.true.!!!!'
-            WRITE(ctmp1,*)  '         bad ln_bfrimp value = .true.' 
+            WRITE(ctmp1,*)  '         set either "ln_zdfexp = .false" or "ln_bfrimp = .false."
             CALL ctl_stop( ctmp1 )
          END IF
@@ -272,4 +273,5 @@
          CALL mpp_max( zmaxbfr )
       ENDIF
+      IF( .NOT.ln_bfrimp) THEN
       IF( lwp .AND. ictu + ictv > 0 ) THEN
          WRITE(numout,*) ' Bottom friction stability check failed at ', ictu, ' U-points '
@@ -278,4 +280,5 @@
          WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary'
       ENDIF
+      ENDIF
       !
       IF( nn_timing == 1 )  CALL timing_stop('zdf_bfr_init')
