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 2837)
+++ branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 2872)
@@ -119,4 +119,5 @@
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
       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                   !   -      -
@@ -124,4 +125,5 @@
       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                   !   -      -
       !!----------------------------------------------------------------------
 
@@ -147,11 +149,11 @@
          hvr_e (:,:) = hvr  (:,:)
          IF( ln_dynvor_een ) THEN
-            ftne(1,:) = 0.e0   ;   ftnw(1,:) = 0.e0   ;   ftse(1,:) = 0.e0   ;   ftsw(1,:) = 0.e0
+            ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;   ftsw(1,:) = 0._wp
             DO jj = 2, jpj
                DO ji = fs_2, jpi   ! vector opt.
-                  ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3.
-                  ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3.
-                  ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3.
-                  ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3.
+                  ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3._wp
+                  ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3._wp
+                  ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3._wp
+                  ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3._wp
                END DO
             END DO
@@ -161,12 +163,12 @@
 
       !                                   !* Local constant initialization
-      z2dt_b = 2.0 * rdt                                    ! baroclinic time step
-      z1_8 = 0.5 * 0.25                                     ! coefficient for vorticity estimates
-      z1_4 = 0.5 * 0.5
-      zraur  = 1. / rau0                                    ! 1 / volumic mass
-      !
-      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
+      z2dt_b = 2.0_wp * rdt                                 ! baroclinic time step
+      z1_8   = 0.5_wp * 0.25_wp                             ! coefficient for vorticity estimates
+      z1_4   = 0.5_wp * 0.5_wp
+      zraur  = 1._wp / rau0                                 ! 1 / volumic mass
+      !
+      zhdiv(:,:) = 0._wp                                    ! barotropic divergence
+      zu_sld = 0._wp   ;   zu_asp = 0._wp                   ! tides trends (lk_tide=F)
+      zv_sld = 0._wp   ;   zv_asp = 0._wp
 
       ! -----------------------------------------------------------------------------
@@ -176,6 +178,6 @@
       !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated)
       !                                   ! --------------------------
-      zua(:,:) = 0.e0   ;   zun(:,:) = 0.e0   ;   ub_b(:,:) = 0.e0
-      zva(:,:) = 0.e0   ;   zvn(:,:) = 0.e0   ;   vb_b(:,:) = 0.e0
+      zua(:,:) = 0._wp   ;   zun(:,:) = 0._wp   ;   ub_b(:,:) = 0._wp
+      zva(:,:) = 0._wp   ;   zvn(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
       !
       DO jk = 1, jpkm1
@@ -205,13 +207,4 @@
       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)
@@ -260,11 +253,38 @@
          DO jj = 2, jpjm1 
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zcu(ji,jj) = zcu(ji,jj) - grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn(ji+1,jj  )    &
-                  &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hu(ji,jj) / e1u(ji,jj)
-               zcv(ji,jj) = zcv(ji,jj) - grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn(ji  ,jj+1)    &
-                  &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hv(ji,jj) / e2v(ji,jj)
+      !     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)
             END DO
          END DO
       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
@@ -279,4 +299,5 @@
       !                                             ! from the barotropic transport trend
       zcoef = -1. / z2dt_b
+      IF(.NOT.ln_bfrimp) THEN
 # if defined key_vectopt_loop
       DO jj = 1, 1
@@ -311,4 +332,5 @@
          END DO
       ENDIF
+      END IF
 
       !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity)
@@ -410,13 +432,6 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! surface pressure gradient
-                  IF( lk_vvl) THEN
-                     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
-                     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
+                  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)
                   ! energy conserving formulation for planetary vorticity term
                   zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
@@ -427,8 +442,21 @@
                   zv_cor =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj)
                   ! 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)   &
+                     &         / ( 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)   &
+                     &         / ( 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)   
+                  ELSE
                   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) )
+                  ENDIF
                END DO
             END DO
@@ -438,13 +466,6 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                    ! surface pressure gradient
-                  IF( lk_vvl) THEN
-                     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
-                     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
+                  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)
                   ! 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)
@@ -453,8 +474,20 @@
                   zv_cor  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj)
                   ! 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)   &
+                     &         / ( 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)   &
+                     &         / ( 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)   
+                  ELSE
                   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) )
+                  ENDIF
                END DO
             END DO
@@ -464,13 +497,6 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! surface pressure gradient
-                  IF( lk_vvl) THEN
-                     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
-                     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
+                  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)
                   ! energy/enstrophy conserving formulation for planetary vorticity term
                   zu_cor = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
@@ -479,8 +505,20 @@
                      &                           + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) * hvr_e(ji,jj)
                   ! 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)   &
+                     &         / ( 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)   &
+                     &         / ( 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)   
+                  ELSE
                   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) )
+                  ENDIF
                END DO
             END DO
