Index: branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 7412)
+++ branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 7421)
@@ -204,8 +204,8 @@
       !!
       !!     ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z).
-      !!         Check value: rho = 1028.21993233072 kg/m^3 for z=3000 dbar, ct=3 Celcius, sa=35.5 g/kg
+      !!         Check value: rho = 1028.21993233072 kg/m^3 for z=3000 dbar, ct=3 Celsius, sa=35.5 g/kg
       !!
       !!     ln_eos80 : polynomial EOS-80 equation of state is used for rho(t,s,z).
-      !!         Check value: rho = 1028.35011066567 kg/m^3 for z=3000 dbar, pt=3 Celcius, sp=35.5 psu
+      !!         Check value: rho = 1028.35011066567 kg/m^3 for z=3000 dbar, pt=3 Celsius, sp=35.5 psu
       !!
       !!     ln_seos : simplified equation of state
@@ -221,5 +221,5 @@
       !!                TEOS-10 Manual, 2010
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius]
       !                                                               ! 2 : salinity               [psu]
       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd   ! in situ density            [-]
@@ -316,5 +316,5 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celcius]
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius]
       !                                                                ! 2 : salinity               [psu]
       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density            [-]
@@ -481,5 +481,5 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
+      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius]
       !                                                           ! 2 : salinity               [psu]
       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                      [m]
@@ -907,6 +907,6 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celcius,psu]
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pab   ! thermal/haline expansion coef.  [Celcius-1,psu-1]
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celsius,psu]
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pab   ! thermal/haline expansion coef.  [Celsius-1,psu-1]
       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::  pn2   ! Brunt-Vaisala frequency squared [1/s^2]
       !
@@ -944,5 +944,5 @@
       !!                 ***  ROUTINE eos_pt_from_ct  ***
       !!
-      !! ** Purpose :   Compute pot.temp. from cons. temp. [Celcius]
+      !! ** Purpose :   Compute pot.temp. from cons. temp. [Celsius]
       !!
       !! ** Method  :   rational approximation (5/3th order) of TEOS-10 algorithm
@@ -952,8 +952,8 @@
       !!                Rational approximation to TEOS10 algorithm (rms error on WOA13 values: 4.0e-5 degC)
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ctmp   ! Cons. Temp [Celcius]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity   [psu]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ctmp   ! Cons. Temp   [Celsius]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity     [psu]
       ! Leave result array automatic rather than making explicitly allocated
-      REAL(wp), DIMENSION(jpi,jpj) ::   ptmp   ! potential temperature [Celcius]
+      REAL(wp), DIMENSION(jpi,jpj) ::   ptmp   ! potential temperature [Celsius]
       !
       INTEGER  ::   ji, jj               ! dummy loop indices
@@ -1003,7 +1003,7 @@
       !!                 ***  ROUTINE eos_fzp  ***
       !!
-      !! ** Purpose :   Compute the freezing point temperature [Celcius]
-      !!
-      !! ** Method  :   UNESCO freezing point (ptf) in Celcius is given by
+      !! ** Purpose :   Compute the freezing point temperature [Celsius]
+      !!
+      !! ** Method  :   UNESCO freezing point (ptf) in Celsius is given by
       !!       ptf(t,z) = (-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s - 7.53e-4*z
       !!       checkvalue: tf=-2.588567 Celsius for s=40psu, z=500m
@@ -1013,5 +1013,5 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   psal   ! salinity   [psu]
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ), OPTIONAL ::   pdep   ! depth      [m]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(out  )           ::   ptf    ! freezing temperature [Celcius]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(out  )           ::   ptf    ! freezing temperature [Celsius]
       !
       INTEGER  ::   ji, jj          ! dummy loop indices
@@ -1056,7 +1056,7 @@
       !!                 ***  ROUTINE eos_fzp  ***
       !!
-      !! ** Purpose :   Compute the freezing point temperature [Celcius]
-      !!
-      !! ** Method  :   UNESCO freezing point (ptf) in Celcius is given by
+      !! ** Purpose :   Compute the freezing point temperature [Celsius]
+      !!
+      !! ** Method  :   UNESCO freezing point (ptf) in Celsius is given by
       !!       ptf(t,z) = (-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s - 7.53e-4*z
       !!       checkvalue: tf=-2.588567 Celsius for s=40psu, z=500m
@@ -1066,5 +1066,5 @@
       REAL(wp), INTENT(in )           ::   psal         ! salinity   [psu]
       REAL(wp), INTENT(in ), OPTIONAL ::   pdep         ! depth      [m]
-      REAL(wp), INTENT(out)           ::   ptf          ! freezing temperature [Celcius]
+      REAL(wp), INTENT(out)           ::   ptf          ! freezing temperature [Celsius]
       !
       REAL(wp) :: zs   ! local scalars
Index: branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 7412)
+++ branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 7421)
@@ -290,13 +290,13 @@
          WRITE(numout,*)
          SELECT CASE ( nadv )
-         CASE( np_NO_adv  )   ;   WRITE(numout,*) '         NO T-S advection'
-         CASE( np_CEN     )   ;   WRITE(numout,*) '         CEN      scheme is used. Horizontal order: ', nn_cen_h,   &
+         CASE( np_NO_adv  )   ;   WRITE(numout,*) '      ===>>   NO T-S advection'
+         CASE( np_CEN     )   ;   WRITE(numout,*) '      ===>>   CEN      scheme is used. Horizontal order: ', nn_cen_h,   &
             &                                                                     ' Vertical   order: ', nn_cen_v
-         CASE( np_FCT     )   ;   WRITE(numout,*) '         FCT      scheme is used. Horizontal order: ', nn_fct_h,   &
+         CASE( np_FCT     )   ;   WRITE(numout,*) '      ===>>   FCT      scheme is used. Horizontal order: ', nn_fct_h,   &
             &                                                                      ' Vertical   order: ', nn_fct_v
-         CASE( np_FCT_zts )   ;   WRITE(numout,*) '         use 2nd order FCT with ', nn_fct_zts,'vertical sub-timestepping'
-         CASE( np_MUS     )   ;   WRITE(numout,*) '         MUSCL    scheme is used'
-         CASE( np_UBS     )   ;   WRITE(numout,*) '         UBS      scheme is used'
-         CASE( np_QCK     )   ;   WRITE(numout,*) '         QUICKEST scheme is used'
+         CASE( np_FCT_zts )   ;   WRITE(numout,*) '      ===>>   use 2nd order FCT with ', nn_fct_zts,'vertical sub-timestepping'
+         CASE( np_MUS     )   ;   WRITE(numout,*) '      ===>>   MUSCL    scheme is used'
+         CASE( np_UBS     )   ;   WRITE(numout,*) '      ===>>   UBS      scheme is used'
+         CASE( np_QCK     )   ;   WRITE(numout,*) '      ===>>   QUICKEST scheme is used'
          END SELECT
       ENDIF
Index: branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen.F90
===================================================================
--- branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen.F90	(revision 7412)
+++ branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen.F90	(revision 7421)
@@ -42,5 +42,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.7 , NEMO Consortium (2014)
-   !! $Id: traadv_cen2.F90 5737 2015-09-13 07:42:41Z gm $
+   !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
Index: branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90
===================================================================
--- branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90	(revision 7412)
+++ branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90	(revision 7421)
@@ -43,4 +43,8 @@
    REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6
 
+   !                                        ! tridiag solver associated indices:
+   INTEGER, PARAMETER ::   np_NH   = 0   ! Neumann homogeneous boundary condition
+   INTEGER, PARAMETER ::   np_CEN2 = 1   ! 2nd order centered  boundary condition
+
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -737,7 +741,7 @@
 
 
-   SUBROUTINE interp_4th_cpt( pt_in, pt_out )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE interp_4th_cpt  ***
+   SUBROUTINE interp_4th_cpt_org( pt_in, pt_out )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE interp_4th_cpt_org  ***
       !! 
       !! **  Purpose :   Compute the interpolation of tracer at w-point
@@ -770,7 +774,7 @@
       END DO
       !
-      jk=2                                            ! Switch to second order centered at top
-      DO jj=1,jpj
-         DO ji=1,jpi
+      jk = 2                                          ! Switch to second order centered at top
+      DO jj = 1, jpj
+         DO ji = 1, jpi
             zwd (ji,jj,jk) = 1._wp
             zwi (ji,jj,jk) = 0._wp
@@ -820,6 +824,178 @@
       END DO
       !    
+   END SUBROUTINE interp_4th_cpt_org
+   
+
+   SUBROUTINE interp_4th_cpt( pt_in, pt_out )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE interp_4th_cpt  ***
+      !! 
+      !! **  Purpose :   Compute the interpolation of tracer at w-point
+      !!
+      !! **  Method  :   4th order compact interpolation
+      !!----------------------------------------------------------------------
+      REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pt_in    ! field at t-point
+      REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pt_out   ! field interpolated at w-point
+      !
+      INTEGER ::   ji, jj, jk   ! dummy loop integers
+      INTEGER ::   ikt, ikb     ! local integers
+      REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt
+      !!----------------------------------------------------------------------
+      !
+      !                      !==  build the three diagonal matrix & the RHS  ==!
+      !
+      DO jk = 3, jpkm1                 ! interior (from jk=3 to jpk-1)
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp                 !       diagonal
+               zwi (ji,jj,jk) =         wmask(ji,jj,jk)                         ! lower diagonal
+               zws (ji,jj,jk) =         wmask(ji,jj,jk)                         ! upper diagonal
+               zwrm(ji,jj,jk) = 3._wp * wmask(ji,jj,jk)                     &   ! RHS
+                  &           *       ( pt_in(ji,jj,jk) + pt_in(ji,jj,jk-1) )
+            END DO
+         END DO
+      END DO
+      !
+!!gm
+!      SELECT CASE( kbc )               !* boundary condition
+!      CASE( np_NH   )   ! Neumann homogeneous at top & bottom
+!      CASE( np_CEN2 )   ! 2nd order centered  at top & bottom
+!      END SELECT
+!!gm  
+      !
+      DO jj = 2, jpjm1                 ! 2nd order centered at top & bottom
+         DO ji = fs_2, fs_jpim1
+            ikt = mikt(ji,jj) + 1            ! w-point below the 1st  wet point
+            ikb = mbkt(ji,jj)                !     -   above the last wet point
+            !
+            zwd (ji,jj,ikt) = 1._wp          ! top
+            zwi (ji,jj,ikt) = 0._wp
+            zws (ji,jj,ikt) = 0._wp
+            zwrm(ji,jj,ikt) = 0.5_wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
+            !
+            zwd (ji,jj,ikb) = 1._wp          ! bottom
+            zwi (ji,jj,ikb) = 0._wp
+            zws (ji,jj,ikb) = 0._wp
+            zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )            
+         END DO
+      END DO   
+      !
+      !                       !==  tridiagonal solver  ==!
+      !
+      DO jj = 2, jpjm1              !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1
+         DO ji = fs_2, fs_jpim1
+            zwt(ji,jj,2) = zwd(ji,jj,2)
+         END DO
+      END DO
+      DO jk = 3, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      !
+      DO jj = 2, jpjm1              !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
+         DO ji = fs_2, fs_jpim1
+            pt_out(ji,jj,2) = zwrm(ji,jj,2)
+         END DO
+      END DO
+      DO jk = 3, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)             
+            END DO
+         END DO
+      END DO
+
+      DO jj = 2, jpjm1              !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk
+         DO ji = fs_2, fs_jpim1
+            pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
+         END DO
+      END DO
+      DO jk = jpk-2, 2, -1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !    
    END SUBROUTINE interp_4th_cpt
-   
+
+
+   SUBROUTINE tridia_solver( pD, pU, pL, pRHS, pt_out , klev )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tridia_solver  ***
+      !! 
+      !! **  Purpose :   solve a symmetric 3diagonal system
+      !!
+      !! **  Method  :   solve M.t_out = RHS(t)  where M is a tri diagonal matrix ( jpk*jpk )
+      !!     
+      !!             ( D_1 U_1  0   0   0  )( t_1 )   ( RHS_1 )
+      !!             ( L_2 D_2 U_2  0   0  )( t_2 )   ( RHS_2 )
+      !!             (  0  L_3 D_3 U_3  0  )( t_3 ) = ( RHS_3 )
+      !!             (        ...          )( ... )   ( ...  )
+      !!             (  0   0   0  L_k D_k )( t_k )   ( RHS_k )
+      !!     
+      !!        M is decomposed in the product of an upper and lower triangular matrix.
+      !!        The tri-diagonals matrix is given as input 3D arrays:   pD, pU, pL 
+      !!        (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal).
+      !!        The solution is pta.
+      !!        The 3d array zwt is used as a work space array.
+      !!----------------------------------------------------------------------
+      REAL(wp),DIMENSION(:,:,:), INTENT(in   ) ::   pD, pU, PL    ! 3-diagonal matrix
+      REAL(wp),DIMENSION(:,:,:), INTENT(in   ) ::   pRHS          ! Right-Hand-Side
+      REAL(wp),DIMENSION(:,:,:), INTENT(  out) ::   pt_out        !!gm field at level=F(klev)
+      INTEGER                  , INTENT(in   ) ::   klev          ! =1 pt_out at w-level 
+      !                                                           ! =0 pt at t-level
+      INTEGER ::   ji, jj, jk   ! dummy loop integers
+      INTEGER ::   kstart       ! local indices
+      REAL(wp),DIMENSION(jpi,jpj,jpk) ::   zwt   ! 3D work array
+      !!----------------------------------------------------------------------
+      !
+      kstart =  1  + klev
+      !
+      DO jj = 2, jpjm1              !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1
+         DO ji = fs_2, fs_jpim1
+            zwt(ji,jj,kstart) = pD(ji,jj,kstart)
+         END DO
+      END DO
+      DO jk = kstart+1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      !
+      DO jj = 2, jpjm1              !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
+         DO ji = fs_2, fs_jpim1
+            pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart)
+         END DO
+      END DO
+      DO jk = kstart+1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)             
+            END DO
+         END DO
+      END DO
+
+      DO jj = 2, jpjm1              !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk
+         DO ji = fs_2, fs_jpim1
+            pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
+         END DO
+      END DO
+      DO jk = jpk-2, kstart, -1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+   END SUBROUTINE tridia_solver
+
    !!======================================================================
 END MODULE traadv_fct
Index: branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mle.F90
===================================================================
--- branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mle.F90	(revision 7412)
+++ branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mle.F90	(revision 7421)
@@ -308,9 +308,9 @@
          WRITE(numout,*)
          IF( ln_mle ) THEN
-            WRITE(numout,*) '   Mixed Layer Eddy induced transport added to tracer advection'
-            IF( nn_mle == 0 )   WRITE(numout,*) '   Fox-Kemper et al 2010 formulation'
-            IF( nn_mle == 1 )   WRITE(numout,*) '   New formulation'
+            WRITE(numout,*) '      ===>>   Mixed Layer Eddy induced transport added to tracer advection'
+            IF( nn_mle == 0 )   WRITE(numout,*) '              Fox-Kemper et al 2010 formulation'
+            IF( nn_mle == 1 )   WRITE(numout,*) '              New formulation'
          ELSE
-            WRITE(numout,*) '   Mixed Layer Eddy parametrisation NOT used'
+            WRITE(numout,*) '      ===>>   Mixed Layer Eddy parametrisation NOT used'
          ENDIF
       ENDIF
@@ -329,6 +329,6 @@
             DO jj = 2, jpj                           ! "coriolis+ time^-1" at u- & v-points
                DO ji = fs_2, jpi   ! vector opt.
-                  zfu = ( ff(ji,jj) + ff(ji,jj-1) ) * 0.5_wp
-                  zfv = ( ff(ji,jj) + ff(ji-1,jj) ) * 0.5_wp
+                  zfu = ( ff_f(ji,jj) + ff_f(ji,jj-1) ) * 0.5_wp
+                  zfv = ( ff_f(ji,jj) + ff_f(ji-1,jj) ) * 0.5_wp
                   rfu(ji,jj) = SQRT(  zfu * zfu + z1_t2 )
                   rfv(ji,jj) = SQRT(  zfv * zfv + z1_t2 )
@@ -347,6 +347,5 @@
          !
          z1_t2 = 1._wp / ( rn_time * rn_time )
-         r1_ft(:,:) = 2._wp * omega * SIN( rad * gphit(:,:) )
-         r1_ft(:,:) = 1._wp / SQRT(  r1_ft(:,:) * r1_ft(:,:) + z1_t2 )
+         r1_ft(:,:) = 1._wp / SQRT(  ff_t(:,:) * ff_t(:,:) + z1_t2  )
          !
       ENDIF
Index: branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90
===================================================================
--- branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90	(revision 7412)
+++ branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90	(revision 7421)
@@ -40,5 +40,5 @@
    
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   upsmsk   !: mixed upstream/centered scheme near some straits
-   !                                                           !  and in closed seas (orca 2 and 4 configurations)
+   !                                                           !  and in closed seas (orca 2 and 1 configurations)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xind     !: mixed upstream/centered index
    
Index: branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90
===================================================================
--- branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 7412)
+++ branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 7421)
@@ -176,5 +176,5 @@
             ! fill sf_chl with sn_chl and control print
             CALL fld_fill( sf_qgh, (/ sn_qgh /), cn_dir, 'tra_bbc_init',   &
-               &          'bottom temperature boundary condition', 'nambbc' )
+               &          'bottom temperature boundary condition', 'nambbc', no_print )
 
             CALL fld_read( nit000, 1, sf_qgh )                         ! Read qgh data
Index: branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
===================================================================
--- branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 7412)
+++ branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 7421)
@@ -519,9 +519,9 @@
          WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
-         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
-         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
-         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
-         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
+         WRITE(numout,*) '   Namelist nambbl : set bbl parameters'
+         WRITE(numout,*) '      diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
+         WRITE(numout,*) '      advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
+         WRITE(numout,*) '      diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
+         WRITE(numout,*) '      advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
       ENDIF
 
@@ -545,5 +545,5 @@
       CALL wrk_dealloc( jpi, jpj, zmbk )
 
-                                        !* sign of grad(H) at u- and v-points
+      !                                 !* sign of grad(H) at u- and v-points
       mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0
       DO jj = 1, jpjm1
@@ -553,5 +553,5 @@
          END DO
       END DO
-
+      !
       DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
          DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
@@ -561,32 +561,9 @@
       END DO
       CALL lbc_lnk( e3u_bbl_0, 'U', 1. )   ;   CALL lbc_lnk( e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
-
+      !
       !                             !* masked diffusive flux coefficients
       ahu_bbl_0(:,:) = rn_ahtbbl * e2_e1u(:,:) * e3u_bbl_0(:,:) * umask(:,:,1)
       ahv_bbl_0(:,:) = rn_ahtbbl * e1_e2v(:,:) * e3v_bbl_0(:,:) * vmask(:,:,1)
 
-
-      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl
-         !
-         SELECT CASE ( jp_cfg )
-         CASE ( 2 )                          ! ORCA_R2
-            ij0 = 102   ;   ij1 = 102              ! Gibraltar enhancement of BBL
-            ii0 = 139   ;   ii1 = 140
-            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
-            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
-            !
-            ij0 =  88   ;   ij1 =  88              ! Red Sea enhancement of BBL
-            ii0 = 161   ;   ii1 = 162
-            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
-            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
-            !
-         CASE ( 4 )                          ! ORCA_R4
-            ij0 =  52   ;   ij1 =  52              ! Gibraltar enhancement of BBL
-            ii0 =  70   ;   ii1 =  71
-            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
-            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
-         END SELECT
-         !
-      ENDIF
       !
       IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init')
Index: branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
===================================================================
--- branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 7412)
+++ branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 7421)
@@ -192,5 +192,5 @@
          WRITE(numout,*)
          WRITE(numout,*) 'tra_dmp_init : T and S newtonian relaxation'
-         WRITE(numout,*) '~~~~~~~~~~~'
+         WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namtra_dmp : set relaxation parameters'
          WRITE(numout,*) '      Apply relaxation   or not       ln_tradmp = ', ln_tradmp
Index: branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90
===================================================================
--- branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 7412)
+++ branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 7421)
@@ -110,8 +110,7 @@
          WRITE(numout,*)
          WRITE(numout,*) 'tra_ldf_init : lateral tracer diffusive operator'
-         WRITE(numout,*) '~~~~~~~~~~~'
+         WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namtra_ldf: already read in ldftra module'
          WRITE(numout,*) '      see ldf_tra_init report for lateral mixing parameters'
-         WRITE(numout,*)
       ENDIF
       !                                   ! use of lateral operator or not
@@ -187,11 +186,11 @@
          WRITE(numout,*)
          SELECT CASE( nldf )
-         CASE( np_no_ldf )   ;   WRITE(numout,*) '   NO lateral diffusion'
-         CASE( np_lap    )   ;   WRITE(numout,*) '   laplacian iso-level operator'
-         CASE( np_lap_i  )   ;   WRITE(numout,*) '   Rotated laplacian operator (standard)'
-         CASE( np_lap_it )   ;   WRITE(numout,*) '   Rotated laplacian operator (triad)'
-         CASE( np_blp    )   ;   WRITE(numout,*) '   bilaplacian iso-level operator'
-         CASE( np_blp_i  )   ;   WRITE(numout,*) '   Rotated bilaplacian operator (standard)'
-         CASE( np_blp_it )   ;   WRITE(numout,*) '   Rotated bilaplacian operator (triad)'
+         CASE( np_no_ldf )   ;   WRITE(numout,*) '      ===>>   NO lateral diffusion'
+         CASE( np_lap    )   ;   WRITE(numout,*) '      ===>>   laplacian iso-level operator'
+         CASE( np_lap_i  )   ;   WRITE(numout,*) '      ===>>   Rotated laplacian operator (standard)'
+         CASE( np_lap_it )   ;   WRITE(numout,*) '      ===>>   Rotated laplacian operator (triad)'
+         CASE( np_blp    )   ;   WRITE(numout,*) '      ===>>   bilaplacian iso-level operator'
+         CASE( np_blp_i  )   ;   WRITE(numout,*) '      ===>>   Rotated bilaplacian operator (standard)'
+         CASE( np_blp_it )   ;   WRITE(numout,*) '      ===>>   Rotated bilaplacian operator (triad)'
          END SELECT
       ENDIF
Index: branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 7412)
+++ branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 7421)
@@ -406,5 +406,5 @@
             !                                        ! fill sf_chl with sn_chl and control print
             CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   &
-               &           'Solar penetration function of read chlorophyll', 'namtra_qsr' )
+               &           'Solar penetration function of read chlorophyll', 'namtra_qsr' , no_print )
          ENDIF
          IF( nqsr == np_RGB ) THEN                 ! constant Chl
Index: branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90
===================================================================
--- branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 7412)
+++ branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 7421)
@@ -141,6 +141,6 @@
          WRITE(numout,*) 'tra_zdf_init : vertical tracer physics scheme'
          WRITE(numout,*) '~~~~~~~~~~~'
-         IF( nzdf ==  0 )   WRITE(numout,*) '              Explicit time-splitting scheme'
-         IF( nzdf ==  1 )   WRITE(numout,*) '              Implicit (euler backward) scheme'
+         IF( nzdf ==  0 )   WRITE(numout,*) '      ===>>   Explicit time-splitting scheme'
+         IF( nzdf ==  1 )   WRITE(numout,*) '      ===>>   Implicit (euler backward) scheme'
       ENDIF
       !
