Index: /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/ARCH/arch-gfortran_ubuntu.fcm
===================================================================
--- /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/ARCH/arch-gfortran_ubuntu.fcm	(revision 2872)
+++ /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/ARCH/arch-gfortran_ubuntu.fcm	(revision 2872)
@@ -0,0 +1,32 @@
+# generic gfortran compiler options for Ubuntu
+# NCDF_INC    netcdf include file
+# NCDF_LIB    netcdf library
+# FC          Fortran compiler command
+# FCFLAGS     Fortran compiler flags
+# FFLAGS      Fortran 77 compiler flags
+# LD          linker
+# LDFLAGS     linker flags, e.g. -L<lib dir> if you have libraries in a
+# FPPFLAGS    pre-processing flags
+# AR          assembler
+# ARFLAGS     assembler flags
+# MK          make
+# USER_INC    additional include files for the compiler,  e.g. -I<include dir>
+# USER_LIB    additional libraries to pass to the linker, e.g. -l<library>
+
+
+%NCDF_INC            -I/usr/include
+%NCDF_LIB            -L/usr/lib -lnetcdff
+
+%FC	             gfortran
+%FCFLAGS             -fdefault-real-8 -O3 -funroll-all-loops -fcray-pointer -ffree-line-length-none
+%FFLAGS              %FCFLAGS
+%LD                  gfortran
+%LDFLAGS
+%FPPFLAGS            -P -C -traditional
+%AR                  ar
+%ARFLAGS             -rs
+#%MK                  gmake
+%MK                  make
+%USER_INC            %NCDF_INC
+%USER_LIB            %NCDF_LIB
+
Index: /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/ARCH/arch-pgf90_mobius.fcm
===================================================================
--- /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/ARCH/arch-pgf90_mobius.fcm	(revision 2872)
+++ /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/ARCH/arch-pgf90_mobius.fcm	(revision 2872)
@@ -0,0 +1,30 @@
+# mpi compiler options for NOCL's cluster mobius
+# NCDF_INC    netcdf include file
+# NCDF_LIB    netcdf library
+# FC          Fortran compiler command
+# FCFLAGS     Fortran compiler flags
+# FFLAGS      Fortran 77 compiler flags
+# LD          linker
+# LDFLAGS     linker flags, e.g. -L<lib dir> if you have libraries in a
+# FPPFLAGS    pre-processing flags
+# AR          assembler
+# ARFLAGS     assembler flags
+# MK          make
+# USER_INC    additional include files for the compiler,  e.g. -I<include dir>
+# USER_LIB    additional libraries to pass to the linker, e.g. -l<library>
+
+
+%NCDF_INC             -I/cm/shared/apps/netcdf/pgi/64/4.1.1/include -I/usr/mpi/qlogic/include
+%NCDF_LIB             -L/cm/shared/apps/pgi/7.1-6/linux86-64/7.1-6/lib -L/usr/mpi/qlogic/lib64  -L/cm/shared/apps/netcdf/pgi/64/4.1.1/lib -lnetcdf
+%FC	             mpif90 -c
+%FCFLAGS             -i4 -r8 -O3 -Mfree  
+%FFLAGS              %FCFLAGS
+%LD                  mpif90
+%LDFLAGS             -Mbounds
+%FPPFLAGS            -P -C -traditional
+%AR                  ar
+%ARFLAGS             -rs 
+%MK                  gmake
+%USER_INC            %NCDF_INC
+%USER_LIB            %NCDF_LIB
+
Index: /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90
===================================================================
--- /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90	(revision 2871)
+++ /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90	(revision 2872)
@@ -13,4 +13,5 @@
    USE dom_oce         ! ocean space and time domain variables 
    USE zdf_oce         ! ocean vertical physics variables
+   USE zdfbfr          ! ocean bottom friction variables
    USE trdmod          ! ocean active dynamics and tracers trends 
    USE trdmod_oce      ! ocean variables trends
@@ -51,37 +52,43 @@
       !!---------------------------------------------------------------------
       !
-      zm1_2dt = - 1._wp / ( 2._wp * rdt )
+      IF( .not. ln_bfrimp) THEN     ! only for explicit bottom friction form
+                                    ! implicit bfr is implemented in dynzdf_imp
+                                    ! H. Liu,  Sept. 2011
 
-      IF( l_trddyn )   THEN                      ! temporary save of ua and va trends
-         ztrduv(:,:,:,1) = ua(:,:,:)
-         ztrduv(:,:,:,2) = va(:,:,:)
-      ENDIF
+        zm1_2dt = - 1._wp / ( 2._wp * rdt )
+
+        IF( l_trddyn )   THEN                      ! temporary save of ua and va trends
+           ztrduv(:,:,:,1) = ua(:,:,:)
+           ztrduv(:,:,:,2) = va(:,:,:)
+        ENDIF
+
 
 # if defined key_vectopt_loop
-      DO jj = 1, 1
-         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
+        DO jj = 1, 1
+           DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
 # else
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1
 # endif
-            ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels
-            ikbv = mbkv(ji,jj)
-            !
-            ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt)
-            ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  bfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu)
-            va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  bfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv)
-         END DO
-      END DO
+              ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels
+              ikbv = mbkv(ji,jj)
+              !
+              ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt)
+              ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  bfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu)
+              va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  bfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv)
+           END DO
+        END DO
 
-      !
-      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
-         ztrduv(:,:,:,1) = ua(:,:,:) - ztrduv(:,:,:,1)
-         ztrduv(:,:,:,2) = va(:,:,:) - ztrduv(:,:,:,2)
-         CALL trd_mod( ztrduv(:,:,:,1), ztrduv(:,:,:,2), jpdyn_trd_bfr, 'DYN', kt )
-      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( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
+           ztrduv(:,:,:,1) = ua(:,:,:) - ztrduv(:,:,:,1)
+           ztrduv(:,:,:,2) = va(:,:,:) - ztrduv(:,:,:,2)
+           CALL trd_mod( ztrduv(:,:,:,1), ztrduv(:,:,:,2), jpdyn_trd_bfr, 'DYN', kt )
+        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' )
+        !
+      ENDIF     ! end explicit bottom friction
    END SUBROUTINE dyn_bfr
 
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 2871)
+++ /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
Index: /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 2871)
+++ /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 2872)
@@ -20,4 +20,5 @@
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! MPP library
+   USE zdfbfr          ! bottom friction setup
 
    IMPLICIT NONE
@@ -61,6 +62,12 @@
       REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step
       !!
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars
+      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) ::   bfr_imp = 1._wp
+      !!----------------------------------------------------------------------
       !!----------------------------------------------------------------------
 
@@ -73,4 +80,5 @@
          IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
+         IF(.NOT.ln_bfrimp) bfr_imp = 0._wp
       ENDIF
 
@@ -80,30 +88,43 @@
 
       ! 1. Vertical diffusion on u
+
+      ! Vertical diffusion on u&v
       ! ---------------------------
       ! 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.
-      !
-      DO jk = 1, jpkm1        ! Matrix
-         DO jj = 2, jpjm1 
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+      !! 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)
+            ! Apply stability criteria on absolute value  : Min abs(bfr) => Max (bfr)
+!               zbfru = MAX( bfrua(ji,jj), fse3u(ji,jj,ikbum1)*zinv )
+               zbfru = bfrua(ji,jj)
+
+            DO jk = 1, ikbum1
                zcoef = - p2dt / fse3u(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)
-               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,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 boudary conditions
             zwi(ji,jj,1) = 0._wp
             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
-         END DO
-      END DO
+
+      ! Bottom boudary conditions  ! H. Liu, May, 2010
+!            !commented out to be consisent with v3.2, h.liu
+!            z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 2.0 * bfr_imp
+            z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 1.0_wp * bfr_imp
+            zws(ji,jj,ikbum1) = 0._wp
+            zwd(ji,jj,ikbum1) = 1._wp - zwi(ji,jj,ikbum1) - z2dtf 
 
       ! Matrix inversion starting from the first level
@@ -121,37 +142,23 @@
       !   The solution (the after velocity) is in ua
       !-----------------------------------------------------------------------
-      !
-      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.
+
+      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
+            DO jk = 2, ikbum1
                zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
             END DO
-         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.
+
+      ! 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   
                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
-         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)
+
+
+      ! thrid 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
@@ -170,29 +177,30 @@
       ! 2. Vertical diffusion on v
       ! ---------------------------
-      ! 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 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        ! Matrix
+
+      DO ji = fs_2, fs_jpim1   ! vector opt.
          DO jj = 2, jpjm1   
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+            ikbvm1 = mbkv(ji,jj)
+            ! Apply stability criteria on absolute value  : Min abs(bfr) => Max (bfr)
+!               zbfrv = MAX( bfrva(ji,jj), fse3v(ji,jj,ikbvm1)*zinv )
+               zbfrv = bfrva(ji,jj)
+
+            DO jk = 1, ikbvm1
                zcoef = -p2dt / fse3v(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)
-               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,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 boudary conditions
             zwi(ji,jj,1) = 0._wp
             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
-         END DO
-      END DO
+
+      ! Bottom boudary conditions  
+      ! Bottom boudary conditions  ! H. Liu, May, 2010
+!            !commented out to be consisent with v3.2, h.liu
+!            z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 2.0 * bfr_imp
+            z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 1.0_wp * bfr_imp
+            zws(ji,jj,ikbvm1) = 0._wp
+            zwd(ji,jj,ikbvm1) = 1._wp - zwi(ji,jj,ikbvm1) - z2dtf 
 
       ! Matrix inversion
@@ -206,42 +214,30 @@
       !        (  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
       !-----------------------------------------------------------------------
-      !
-      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.
+
+      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
+            DO jk = 2, ikbvm1
                zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
             END DO
-         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.
+
+      ! 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
                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
-         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
+
+      ! thrid 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
@@ -258,4 +254,5 @@
       IF( wrk_not_released(3, 3) )   CALL ctl_stop('dyn_zdf_imp: failed to release workspace array')
       !
+
    END SUBROUTINE dyn_zdf_imp
 
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 2871)
+++ /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 2872)
@@ -36,4 +36,5 @@
    REAL(wp) ::   rn_bfrien = 30._wp      ! local factor to enhance coefficient bfri
    LOGICAL  ::   ln_bfr2d  = .false.     ! logical switch for 2D enhancement
+   LOGICAL, PUBLIC  ::   ln_bfrimp = .false.     ! logical switch for implicit bottom friction
    
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  bfrcoef2d   ! 2D bottom drag coefficient
@@ -142,5 +143,5 @@
       REAL(wp) ::  zfru, zfrv         !    -         -
       !!
-      NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfeb2, ln_bfr2d, rn_bfrien
+      NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfeb2, ln_bfr2d, rn_bfrien, ln_bfrimp
       !!----------------------------------------------------------------------
 
Index: /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/par_AMM12.h90
===================================================================
--- /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/par_AMM12.h90	(revision 2872)
+++ /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/par_AMM12.h90	(revision 2872)
@@ -0,0 +1,83 @@
+   !!---------------------------------------------------------------------
+   !!                     ***  par_ATLARC_R033.h90  ***  
+   !!   Ocean Domain : 0.33 degrees resolution atlantic and arctic oceans
+   !!                  (ATLARC_R033 configuration)
+   !!---------------------------------------------------------------------
+   CHARACTER (len=16), PARAMETER ::   &
+      cp_cfg = "amm"         !: name of the configuration
+   INTEGER, PARAMETER ::     &
+      jp_cfg = 011  ,        &  !: resolution of the configuration (degrees)
+      ! Original data size
+      jpidta  = 198,        &  !: first horizontal dimension > or = to jpi
+      jpjdta  = 224,        &  !: second                     > or = to jpj
+      jpkdta  = 33,         &  !: number of levels           > or = to jpk
+      ! total domain matrix size
+      jpiglo  = jpidta,      &  !: first  dimension of global domain --> i
+      jpjglo  = jpjdta,      &  !: second dimension of global domain --> j
+      ! starting position of the zoom
+      jpizoom =   1   ,      &  !: left bottom (i,j) indices of the zoom
+      jpjzoom =   1   ,      &  !: in data indices
+      ! Domain characteristics
+      jperio  =    0  ,      &  !: lateral cond. type (between 0 and 6)
+      jpisl   =  150  ,      &  !: number of islands (irrelevent)
+      jpnisl  = 3000            !: maximum number of points per island
+
+   !!  Values set to pp_not_used indicates that this parameter is not used in THIS config.
+   !!  Values set to pp_to_be_computed  indicates that variables will be computed in domzgr
+   REAL(wp),PARAMETER      ::  pp_not_used = 999999_wp , &
+      &                        pp_to_be_computed = 0._wp
+   !!
+   !! Coefficients associated with the horizontal coordinate system (jphgr_msh /= 0 )
+   !!
+   INTEGER, PARAMETER ::     & !
+      jphgr_msh = 0            !: type of horizontal mesh
+      !                        !  = 0 curvilinear coordinate on the sphere
+      !                        !      read in coordinate.nc file
+      !                        !  = 1 geographical mesh on the sphere
+      !                        !      with regular grid-spacing
+      !                        !  = 2 f-plane with regular grid-spacing
+      !                        !  = 3 beta-plane with regular grid-spacing
+      !                        !  = 4 Mercator grid with T/U point at the equator  with
+      !                        !      isotropic resolution (e1_deg)
+
+      !   ppglam0 , ppgphi0: coordinates of the lower leftmost T point of the grid.
+      !   The mercator grid starts only approximately at gphi0 because
+      !   of the constraint that the equator be a T point.
+   REAL(wp), PARAMETER ::       &  !
+      ppglam0  = pp_not_used,   &  !: longitude of first raw and column T-point (jphgr_msh = 1)
+      ppgphi0  = pp_not_used,   &  !: latitude  of first raw and column T-point (jphgr_msh = 1)
+      !                            !  latitude for the Coriolis or Beta parameter (jphgr_msh = 2 or 3)
+      ppe1_deg = pp_not_used,   &  !: zonal      grid-spacing (degrees)
+      ppe2_deg = pp_not_used,   &  !: meridional grid-spacing (degrees)
+      !
+      ppe1_m   = pp_not_used,   &  !: zonal      grid-spacing (meters )
+      ppe2_m   = pp_not_used       !: meridional grid-spacing (meters )
+
+   !!  Coefficients associated with the vertical coordinate system
+
+   REAL(wp), PARAMETER  ::       &
+      &     ppsur = pp_to_be_computed ,  &  !: Computed in domzgr, set ppdzmin and pphmax below
+      &     ppa0  = pp_to_be_computed ,  &  !:    "           "
+      &     ppa1  = pp_to_be_computed ,  &  !:    "           "
+      !
+      &     ppkth =  23.563_wp        ,  &  !: (non dimensional): gives the approximate
+      !                                     !: layer number above which  stretching will
+      !                                     !: be maximum. Usually of order jpk/2.
+      &     ppacr =    9.00000000000_wp     !: (non dimensional): stretching factor
+      !                                     !: for the grid. The highest zacr, the smallest
+      !                                     !: the stretching.
+
+   !!
+   !!  If both ppa0 ppa1 and ppsur are specified to 0, then
+   !!  they are computed from ppdzmin, pphmax , ppkth, ppacr in dom_zgr
+   !!
+   REAL(wp), PARAMETER ::        &
+      &     ppdzmin = 6._wp           ,  &  !: (meters) vertical thickness of the top layer
+      &     pphmax  = 5720._wp              !: (meters) Maximum depth of the ocean gdepw(jpk)
+   LOGICAL, PARAMETER ::         &
+      &     ldbletanh = .FALSE.
+   REAL(wp), PARAMETER ::        &
+      &     ppa2    = pp_not_used           ,  &  !: Double tanh function parameters
+      &     ppkth2  = pp_not_used           ,  &  !:
+      &     ppacr2  = pp_not_used                 !:
+   !!---------------------------------------------------------------------
Index: /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/par_AMM7.h90
===================================================================
--- /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/par_AMM7.h90	(revision 2872)
+++ /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/par_AMM7.h90	(revision 2872)
@@ -0,0 +1,83 @@
+   !!---------------------------------------------------------------------
+   !!                     ***  par_ATLARC_R033.h90  ***  
+   !!   Ocean Domain : 0.33 degrees resolution atlantic and arctic oceans
+   !!                  (ATLARC_R033 configuration)
+   !!---------------------------------------------------------------------
+   CHARACTER (len=16), PARAMETER ::   &
+      cp_cfg = "amm"         !: name of the configuration
+   INTEGER, PARAMETER ::     &
+      jp_cfg = 011  ,        &  !: resolution of the configuration (degrees)
+      ! Original data size
+      jpidta  = 297,        &  !: first horizontal dimension > or = to jpi
+      jpjdta  = 375,        &  !: second                     > or = to jpj
+      jpkdta  = 33,         &  !: number of levels           > or = to jpk
+      ! total domain matrix size
+      jpiglo  = jpidta,      &  !: first  dimension of global domain --> i
+      jpjglo  = jpjdta,      &  !: second dimension of global domain --> j
+      ! starting position of the zoom
+      jpizoom =   1   ,      &  !: left bottom (i,j) indices of the zoom
+      jpjzoom =   1   ,      &  !: in data indices
+      ! Domain characteristics
+      jperio  =    0  ,      &  !: lateral cond. type (between 0 and 6)
+      jpisl   =  150  ,      &  !: number of islands (irrelevent)
+      jpnisl  = 3000            !: maximum number of points per island
+
+   !!  Values set to pp_not_used indicates that this parameter is not used in THIS config.
+   !!  Values set to pp_to_be_computed  indicates that variables will be computed in domzgr
+   REAL(wp),PARAMETER      ::  pp_not_used = 999999_wp , &
+      &                        pp_to_be_computed = 0._wp
+   !!
+   !! Coefficients associated with the horizontal coordinate system (jphgr_msh /= 0 )
+   !!
+   INTEGER, PARAMETER ::     & !
+      jphgr_msh = 0            !: type of horizontal mesh
+      !                        !  = 0 curvilinear coordinate on the sphere
+      !                        !      read in coordinate.nc file
+      !                        !  = 1 geographical mesh on the sphere
+      !                        !      with regular grid-spacing
+      !                        !  = 2 f-plane with regular grid-spacing
+      !                        !  = 3 beta-plane with regular grid-spacing
+      !                        !  = 4 Mercator grid with T/U point at the equator  with
+      !                        !      isotropic resolution (e1_deg)
+
+      !   ppglam0 , ppgphi0: coordinates of the lower leftmost T point of the grid.
+      !   The mercator grid starts only approximately at gphi0 because
+      !   of the constraint that the equator be a T point.
+   REAL(wp), PARAMETER ::       &  !
+      ppglam0  = pp_not_used,   &  !: longitude of first raw and column T-point (jphgr_msh = 1)
+      ppgphi0  = pp_not_used,   &  !: latitude  of first raw and column T-point (jphgr_msh = 1)
+      !                            !  latitude for the Coriolis or Beta parameter (jphgr_msh = 2 or 3)
+      ppe1_deg = pp_not_used,   &  !: zonal      grid-spacing (degrees)
+      ppe2_deg = pp_not_used,   &  !: meridional grid-spacing (degrees)
+      !
+      ppe1_m   = pp_not_used,   &  !: zonal      grid-spacing (meters )
+      ppe2_m   = pp_not_used       !: meridional grid-spacing (meters )
+
+   !!  Coefficients associated with the vertical coordinate system
+
+   REAL(wp), PARAMETER  ::       &
+      &     ppsur = pp_to_be_computed ,  &  !: Computed in domzgr, set ppdzmin and pphmax below
+      &     ppa0  = pp_to_be_computed ,  &  !:    "           "
+      &     ppa1  = pp_to_be_computed ,  &  !:    "           "
+      !
+      &     ppkth =  23.563_wp        ,  &  !: (non dimensional): gives the approximate
+      !                                     !: layer number above which  stretching will
+      !                                     !: be maximum. Usually of order jpk/2.
+      &     ppacr =    9.00000000000_wp     !: (non dimensional): stretching factor
+      !                                     !: for the grid. The highest zacr, the smallest
+      !                                     !: the stretching.
+
+   !!
+   !!  If both ppa0 ppa1 and ppsur are specified to 0, then
+   !!  they are computed from ppdzmin, pphmax , ppkth, ppacr in dom_zgr
+   !!
+   REAL(wp), PARAMETER ::        &
+      &     ppdzmin = 6._wp           ,  &  !: (meters) vertical thickness of the top layer
+      &     pphmax  = 5720._wp              !: (meters) Maximum depth of the ocean gdepw(jpk)
+   LOGICAL, PARAMETER ::         &
+      &     ldbletanh = .FALSE.
+   REAL(wp), PARAMETER ::        &
+      &     ppa2    = pp_not_used           ,  &  !: Double tanh function parameters
+      &     ppkth2  = pp_not_used           ,  &  !:
+      &     ppacr2  = pp_not_used                 !:
+   !!---------------------------------------------------------------------
Index: /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/par_oce.F90
===================================================================
--- /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/par_oce.F90	(revision 2871)
+++ /branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/par_oce.F90	(revision 2872)
@@ -81,4 +81,14 @@
    !!---------------------------------------------------------------------
 #             include "par_POMME_R025.h90"
+#elif defined key_amm12 
+   !!--------------------------------------------------------------------- 
+   !!   'key_amm12'         :   Atlantic Margin Model (~12km)     :      AMM 
+   !!--------------------------------------------------------------------- 
+#             include "par_AMM12.h90" 
+#elif defined key_amm7 
+   !!--------------------------------------------------------------------- 
+   !!   'key_amm7'         :   Atlantic Margin Model (~7km)     :      AMM 
+   !!--------------------------------------------------------------------- 
+#             include "par_AMM7.h90" 
 #else
    !!---------------------------------------------------------------------
