Index: /branches/2011/dev_r2802_NOCL_bfrimp/DOC/TexFiles/Chapters/Chap_ZDF.tex
===================================================================
--- /branches/2011/dev_r2802_NOCL_bfrimp/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 3056)
+++ /branches/2011/dev_r2802_NOCL_bfrimp/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 3057)
@@ -1,4 +1,4 @@
 % ================================================================
-% Chapter Ñ Vertical Ocean Physics (ZDF)
+% Chapter  Vertical Ocean Physics (ZDF)
 % ================================================================
 \chapter{Vertical Ocean Physics (ZDF)}
@@ -539,5 +539,5 @@
 the clipping factor is of crucial importance for the entrainment depth predicted in 
 stably stratified situations, and that its value has to be chosen in accordance 
-with the algebraic model for the turbulent ßuxes. The clipping is only activated 
+with the algebraic model for the turbulent uxes. The clipping is only activated 
 if \np{ln\_length\_lim}=true, and the $c_{lim}$ is set to the \np{rn\_clim\_galp} value.
 
@@ -982,4 +982,57 @@
 
 % -------------------------------------------------------------------------------------------------------------
+%       Implicit Bottom Friction
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Implicit Bottom Friction}
+\label{ZDF_bfr_imp}
+
+An Implicit form of bottom friction has been devised in the NEMO to improve
+model stability. We recommend this option for shelf sea and coastal ocean applications, especially 
+for split-explicit time splitting. This option can be invoked by setting \np{ln\_bfrimp} 
+to \textit{true} in the \textit{nambfr} namelist and \np{ln\_zdfexp} to \textit{false} 
+in the \textit{namzdf} namelist. 
+
+This implementation is realised in \mdl{dynzdf\_imp} and \mdl{dynspg\_ts}. In \mdl{dynzdf\_imp}, the 
+bottom boundary condition is implemented implicitly.
+
+\begin{equation} \label{Eq_dynzdf_bfr}
+\left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{mbk}
+	 = \binom{c_{b}^{u}u^{n+1}_{mbk}}{c_{b}^{v}v^{n+1}_{mbk}}
+\end{equation}
+
+where $mbk$ is the layer number of the bottom layer. superscript $n+1$ means the velocity used in the
+friction formulea is to be calculated, so, it is implicit.
+
+If split-explicit time splitting is used, care must be taken to avoid the double counting of
+the bottom friction in the 2-D barotropic momentum equations. As NEMO only updates the barotropic 
+pressure gradient and Coriolis' forcing terms in the 2-D barotropic calculation, we need to remove
+the bottom friction induced by these two terms which has been included in the 3-D momentum trend 
+and update it with the latest value. On the other hand, the bottom friction contributed by the
+other terms (e.g. the advection term, viscosity term) has been included in the 3-D momentum equations
+and should not be added in the 2-D barotropic mode.
+
+The implementation of the implicit bottom friction in \mdl{dynspg\_ts} is done in two steps as the
+following:
+
+\begin{equation} \label{Eq_dynspg_ts_bfr1}
+\frac{\textbf{U}_{med}-\textbf{U}^{m-1}}{2\Delta t}=-g\nabla\eta-f\textbf{k}\times\textbf{U}^{m}+c_{b}
+\left(\textbf{U}_{med}-\textbf{U}^{m-1}\right)
+\end{equation}
+\begin{equation} \label{Eq_dynspg_ts_bfr2}
+\frac{\textbf{U}^{m+1}-\textbf{U}_{med}}{2\Delta t}=\textbf{T}+
+\left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{U}^{'}\right)-
+2\Delta t_{bc}c_{b}\left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{u}_{b}\right)
+\end{equation}
+
+where $\textbf{T}$ is the vertical integrated 3-D momentum trend. We assume the leap-frog time-stepping
+is used here. $\Delta t$ is the barotropic mode time step and $\Delta t_{bc}$ is the baroclinic mode time step.
+ $c_{b}$ is the friction coefficient. $\eta$ is the sea surface level calculated in the barotropic loops
+while $\eta^{'}$ is the sea surface level used in the 3-D baroclinic mode. $\textbf{u}_{b}$ is the bottom
+layer horizontal velocity.
+
+
+
+
+% -------------------------------------------------------------------------------------------------------------
 %       Bottom Friction with split-explicit time splitting
 % -------------------------------------------------------------------------------------------------------------
@@ -1091,5 +1144,5 @@
 The essential goal of the parameterization is to represent the momentum 
 exchange between the barotropic tides and the unrepresented internal waves 
-induced by the tidal ßow over rough topography in a stratified ocean. 
+induced by the tidal ow over rough topography in a stratified ocean. 
 In the current version of \NEMO, the map is built from the output of 
 the barotropic global ocean tide model MOG2D-G \citep{Carrere_Lyard_GRL03}.
Index: /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 3056)
+++ /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 3057)
@@ -120,10 +120,10 @@
       INTEGER  ::   icycle           ! local scalar
       INTEGER  ::   ikbu, ikbv       ! local scalar
-      REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b     ! local 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) ::   ua_btm, va_btm                   !   -      -
+      REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b, z2dt_bf     ! local 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) ::   ua_btm, va_btm                            !   -      -
       !!----------------------------------------------------------------------
 
@@ -171,4 +171,10 @@
       zu_sld = 0._wp   ;   zu_asp = 0._wp                   ! tides trends (lk_tide=F)
       zv_sld = 0._wp   ;   zv_asp = 0._wp
+
+      IF( kt == nit000 .AND. neuler == 0) THEN              ! for implicit bottom friction
+        z2dt_bf = rdt
+      ELSE
+        z2dt_bf = 2.0_wp * rdt
+      ENDIF
 
       ! -----------------------------------------------------------------------------
@@ -207,4 +213,13 @@
       END DO
 
+      !                                   !* baroclinic momentum trend (remove the vertical mean trend)
+      DO jk = 1, jpkm1                    ! --------------------------
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
+               va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
+            END DO
+         END DO
+      END DO
 
       !                                   !* barotropic Coriolis trends * H (vorticity scheme dependent)
@@ -253,5 +268,4 @@
          DO jj = 2, jpjm1 
             DO ji = fs_2, fs_jpim1   ! vector opt.
-      !     remove the rhd term according to J. Chanut's suggestion            
                zcu(ji,jj) = zcu(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * hu(ji,jj) / e1u(ji,jj)
                zcv(ji,jj) = zcv(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * hv(ji,jj) / e2v(ji,jj)
@@ -260,37 +274,9 @@
       ENDIF
 
-     IF(ln_bfrimp.AND.lk_vvl) THEN
-!    Remove the bottom stress trend from 3-D sea surface level gradient and Coriolis force           
-!    H. Liu, Sept. 2011
-      DO jj = 2, jpjm1         
-         DO ji = fs_2, fs_jpim1
-            ikbu = mbku(ji,jj)
-            ikbv = mbkv(ji,jj)
-            ua_btm = zcu(ji,jj) * z2dt_b * hur(ji,jj) * umask (ji,jj,ikbu)
-            va_btm = zcv(ji,jj) * z2dt_b * hvr(ji,jj) * vmask (ji,jj,ikbv)
-
-            ua(ji,jj,ikbu) = ua(ji,jj,ikbu) - bfrua(ji,jj) * ua_btm / fse3u(ji,jj,ikbu)
-            va(ji,jj,ikbv) = va(ji,jj,ikbv) - bfrva(ji,jj) * va_btm / fse3v(ji,jj,ikbv)
-            
-            zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm
-            zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm
-         END DO
-      END DO
-     END IF
-      !                                   !* baroclinic momentum trend (remove the vertical mean trend)
-      DO jk = 1, jpkm1                    ! --------------------------
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
-               va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
-            END DO
-         END DO
-      END DO
-
       DO jj = 2, jpjm1                             ! Remove coriolis term (and possibly spg) from barotropic trend
          DO ji = fs_2, fs_jpim1
-            zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)
-            zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)
-         END DO
+             zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)
+             zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)
+          END DO
       END DO
 
@@ -299,39 +285,57 @@
       !                                             ! from the barotropic transport trend
       zcoef = -1. / z2dt_b
-      IF(.NOT.ln_bfrimp) THEN
+
+      IF(ln_bfrimp) THEN
+      !                                   ! Remove the bottom stress trend from 3-D sea surface level gradient
+      !                                   ! and Coriolis forcing in case of 3D semi-implicit bottom friction 
+        DO jj = 2, jpjm1         
+           DO ji = fs_2, fs_jpim1
+              ikbu = mbku(ji,jj)
+              ikbv = mbkv(ji,jj)
+              ua_btm = zcu(ji,jj) * z2dt_bf * hur(ji,jj) * umask (ji,jj,ikbu)
+              va_btm = zcv(ji,jj) * z2dt_bf * hvr(ji,jj) * vmask (ji,jj,ikbv)
+
+              zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm
+              zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm
+           END DO
+        END DO
+
+      ELSE
+
 # if defined key_vectopt_loop
-      DO jj = 1, 1
-         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+        DO jj = 1, 1
+           DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
 # else
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1
 # endif
             ! Apply stability criteria for bottom friction
             !RBbug for vvl and external mode we may need to use varying fse3
             !!gm  Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx.
-            zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef  )
-            zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef  )
-         END DO
-      END DO
-
-      IF( lk_vvl ) THEN
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj)   &
-                  &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) )
-               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(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) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj)
-               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj)
-            END DO
-         END DO
-      ENDIF
-      END IF
-
+              zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef  )
+              zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef  )
+           END DO
+        END DO
+
+        IF( lk_vvl ) THEN
+           DO jj = 2, jpjm1
+              DO ji = fs_2, fs_jpim1   ! vector opt.
+                 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj)   &
+                    &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) )
+                 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(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) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj)
+                 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj)
+              END DO
+           END DO
+        ENDIF
+      END IF    ! end (ln_bfrimp)
+
+                    
       !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity)
       !                                   ! -------------------------- 
@@ -443,12 +447,11 @@
                   ! after velocities with implicit bottom friction
 
-                  IF(ln_bfrimp) THEN
-                  !   A new method to implement the implicit bottom friction. 
-                  !   H. Liu
-                  !   May 2010
-                  ua_e(ji,jj) =  zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )  * umask(ji,jj,1)   &
+                  IF(ln_bfrimp) THEN      !implicit bottom friction
+                  ua_e(ji,jj) = ((1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj)) * zub_e(ji,jj) +   &
+                     &          z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )) * umask(ji,jj,1)   &
                      &         / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
                   ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)   
-                  va_e(ji,jj) =  zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )  * vmask(ji,jj,1)   &
+                  va_e(ji,jj) = ((1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj)) * zvb_e(ji,jj) +   &
+                     &          z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )) * vmask(ji,jj,1)   &
                      &         / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
                   va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)   
@@ -478,8 +481,10 @@
                   !   H. Liu
                   !   May 2010
-                  ua_e(ji,jj) =  zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )  * umask(ji,jj,1)   &
+                  ua_e(ji,jj) = ((1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj)) * zub_e(ji,jj) +   & 
+                     &          z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )) * umask(ji,jj,1)   &
                      &         / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
                   ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)   
-                  va_e(ji,jj) =  zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )  * vmask(ji,jj,1)   &
+                  va_e(ji,jj) = ((1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj)) * zvb_e(ji,jj) +   &
+                     &          z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )) * vmask(ji,jj,1)   &
                      &         / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
                   va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)   
@@ -509,8 +514,10 @@
                   !   H. Liu
                   !   May 2010
-                  ua_e(ji,jj) =  zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )  * umask(ji,jj,1)   &
+                  ua_e(ji,jj) = ((1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj)) * zub_e(ji,jj) +   &
+                     &          z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )) * umask(ji,jj,1)   &
                      &         / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
                   ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)   
-                  va_e(ji,jj) =  zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )  * vmask(ji,jj,1)   &
+                  va_e(ji,jj) = ((1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj)) * zvb_e(ji,jj) +   &
+                     &          z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )) * vmask(ji,jj,1)   &
                      &         / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
                   va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)   
Index: /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
===================================================================
--- /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 3056)
+++ /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 3057)
@@ -157,4 +157,15 @@
       !                              ! allocate zdfbfr arrays
       IF( zdf_bfr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_bfr_init : unable to allocate arrays' )
+
+      !                              ! Make sure ln_zdfexp=.false. when use implicit bfr
+      IF( ln_bfrimp .AND. ln_zdfexp ) THEN
+         IF(lwp) THEN
+            WRITE(numout,*)
+            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.' 
+            CALL ctl_stop( ctmp1 )
+         END IF
+      END IF
 
       SELECT CASE (nn_bfr)
