Index: branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_2/limhdf_2.F90
===================================================================
--- branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_2/limhdf_2.F90	(revision 4619)
+++ branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_2/limhdf_2.F90	(revision 4680)
@@ -86,10 +86,8 @@
       zdiv0(:, 1 ) = 0._wp
       zdiv0(:,jpj) = 0._wp
-      IF( .NOT.lk_vopt_loop ) THEN
-         zflu (jpi,:) = 0._wp   
-         zflv (jpi,:) = 0._wp
-         zdiv0(1,  :) = 0._wp
-         zdiv0(jpi,:) = 0._wp
-      ENDIF
+      zflu (jpi,:) = 0._wp   
+      zflv (jpi,:) = 0._wp
+      zdiv0(1,  :) = 0._wp
+      zdiv0(jpi,:) = 0._wp
 
       zconv = 1._wp           !==  horizontal diffusion using a Crant-Nicholson scheme  ==!
Index: branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_2/limistate_2.F90
===================================================================
--- branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_2/limistate_2.F90	(revision 4619)
+++ branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_2/limistate_2.F90	(revision 4680)
@@ -14,5 +14,4 @@
    !!   'key_lim2' :                                  LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
    !!   lim_istate_2      :  Initialisation of diagnostics ice variables
    !!   lim_istate_init_2 :  initialization of ice state and namelist read
@@ -34,6 +33,6 @@
    PUBLIC lim_istate_2      ! routine called by lim_init_2.F90
 
-   !!! **  namelist (namiceini) **
-   LOGICAL  ::   ln_limini   !: Ice initialization state
+   !                        !! **  namelist (namiceini) **
+   LOGICAL  ::   ln_limini   ! Ice initialization state
    REAL(wp) ::   ttest       ! threshold water temperature for initial sea ice
    REAL(wp) ::   hninn       ! initial snow thickness in the north
@@ -51,5 +50,4 @@
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -71,5 +69,5 @@
       IF( .NOT. ln_limini ) THEN  
          
-         tfu(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius]
+         tfu(:,:) = eos_fzp( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius]
 
          DO jj = 1, jpj
Index: branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
===================================================================
--- branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90	(revision 4619)
+++ branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90	(revision 4680)
@@ -83,10 +83,8 @@
       zdiv0(:, 1 ) = 0._wp
       zdiv0(:,jpj) = 0._wp
-      IF( .NOT.lk_vopt_loop ) THEN
-         zflu (jpi,:) = 0._wp   
-         zflv (jpi,:) = 0._wp
-         zdiv0(1,  :) = 0._wp
-         zdiv0(jpi,:) = 0._wp
-      ENDIF
+      zflu (jpi,:) = 0._wp   
+      zflv (jpi,:) = 0._wp
+      zdiv0(1,  :) = 0._wp
+      zdiv0(jpi,:) = 0._wp
 
       zconv = 1._wp           !==  horizontal diffusion using a Crant-Nicholson scheme  ==!
Index: branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 4619)
+++ branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 4680)
@@ -36,5 +36,4 @@
    PUBLIC   lim_istate      ! routine called by lim_init.F90
 
-   !! * Module variables
    !                          !!** init namelist (namiceini) **
    REAL(wp) ::   ttest   ! threshold water temperature for initial sea ice
@@ -53,5 +52,4 @@
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -121,5 +119,5 @@
 
       ! Basal temperature is set to the freezing point of seawater in Celsius
-      t_bo(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius]
+      t_bo(:,:) = eos_fzp( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius]
 
       DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest
Index: branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
===================================================================
--- branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 4619)
+++ branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 4680)
@@ -551,6 +551,6 @@
       DO jj = 1, jpjm1
          DO ji = 1, jpim1
-            mgrhu(ji,jj) = INT(  SIGN( 1.e0, fsdept_0(ji+1,jj,mbkt(ji+1,jj)) - fsdept_0(ji,jj,mbkt(ji,jj)) )  )
-            mgrhv(ji,jj) = INT(  SIGN( 1.e0, fsdept_0(ji,jj+1,mbkt(ji,jj+1)) - fsdept_0(ji,jj,mbkt(ji,jj)) )  )
+            mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
+            mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
          END DO
       END DO
@@ -558,6 +558,6 @@
       DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
          DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
-            e3u_bbl_0(ji,jj) = MIN( fse3u_0(ji,jj,mbkt(ji+1,jj  )), fse3u_0(ji,jj,mbkt(ji,jj)) )
-            e3v_bbl_0(ji,jj) = MIN( fse3v_0(ji,jj,mbkt(ji  ,jj+1)), fse3v_0(ji,jj,mbkt(ji,jj)) )
+            e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
+            e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
          END DO
       END DO
Index: branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90
===================================================================
--- branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90	(revision 4619)
+++ branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90	(revision 4680)
@@ -2,5 +2,5 @@
    !!==============================================================================
    !!                       ***  MODULE  tranpc  ***
-   !! Ocean active tracers:  non penetrative convection scheme
+   !! Ocean active tracers:  non penetrative convective adjustment scheme
    !!==============================================================================
    !! History :  1.0  ! 1990-09  (G. Madec)  Original code
@@ -9,4 +9,5 @@
    !!            3.0  ! 2008-06  (G. Madec)  applied on ta, sa and called before tranxt in step.F90
    !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
+   !!            3.7  ! 2014-06  (L. Brodeau) new algorithm based on local Brunt-Vaisala freq.
    !!----------------------------------------------------------------------
 
@@ -14,10 +15,11 @@
    !!   tra_npc : apply the non penetrative convection scheme
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and active tracers 
+   USE oce             ! ocean dynamics and active tracers
    USE dom_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
    USE zdf_oce         ! ocean vertical physics
    USE trd_oce         ! ocean active tracer trends
    USE trdtra          ! ocean active tracer trends
-   USE eosbn2          ! equation of state (eos routine) 
+   USE eosbn2          ! equation of state (eos routine)
    USE lbclnk          ! lateral boundary conditions (or mpp link)
    USE in_out_manager  ! I/O manager
@@ -29,11 +31,12 @@
    PRIVATE
 
-   PUBLIC   tra_npc       ! routine called by step.F90
+   PUBLIC   tra_npc    ! routine called by step.F90
 
    !! * Substitutions
 #  include "domzgr_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$ 
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
+   !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -44,44 +47,48 @@
       !!                  ***  ROUTINE tranpc  ***
       !!
-      !! ** Purpose :   Non penetrative convective adjustment scheme. solve 
+      !! ** Purpose : Non-penetrative convective adjustment scheme. solve
       !!      the static instability of the water column on after fields
       !!      while conserving heat and salt contents.
       !!
-      !! ** Method  :   The algorithm used converges in a maximium of jpk 
-      !!      iterations. instabilities are treated when the vertical density
-      !!      gradient is less than 1.e-5.
-      !!      l_trdtra=T: the trend associated with this algorithm is saved.
+      !! ** Method  : updated algorithm able to deal with non-linear equation of state
+      !!              (i.e. static stability computed locally)
       !!
       !! ** Action  : - (ta,sa) after the application od the npc scheme
-      !!              - save the associated trends (ttrd,strd) ('key_trdtra')
+      !!              - send the associated trends for on-line diagnostics (l_trdtra=T)
       !!
-      !! References : Madec, et al., 1991, JPO, 21, 9, 1349-1371.
+      !! References :     Madec, et al., 1991, JPO, 21, 9, 1349-1371.
       !!----------------------------------------------------------------------
-      !
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   inpcc        ! number of statically instable water column
-      INTEGER  ::   inpci        ! number of iteration for npc scheme
-      INTEGER  ::   jiter, jkdown, jkp        ! ???
-      INTEGER  ::   ikbot, ik, ikup, ikdown   ! ???
-      REAL(wp) ::   ze3tot, zta, zsa, zraua, ze3dwn
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zwx, zwy, zwz
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds, zrhop
+      INTEGER  ::   jiter, ikbot, ik, ikup, ikdown, ilayer, ikm   ! local integers
+      LOGICAL  ::   l_bottom_reached, l_column_treated
+      REAL(wp) ::   zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z
+      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_r2dt
+      REAL(wp), POINTER, DIMENSION(:)       ::   zvn2   ! vertical profile of N2 at 1 given point...
+      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvts   ! vertical profile of T and S at 1 given point...
+      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvab   ! vertical profile of alpha and beta
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zn2    ! N^2 
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zab    ! alpha and beta
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   ztrdt, ztrds   ! 3D workspace
+      !
+      !!LB debug:
+      LOGICAL, PARAMETER :: l_LB_debug = .FALSE.
+      INTEGER :: ilc1, jlc1, klc1, nncpu
+      LOGICAL :: lp_monitor_point
+      !!LB debug.
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('tra_npc')
       !
-      CALL wrk_alloc(jpi, jpj, jpk, zrhop )
-      CALL wrk_alloc(jpi, jpk, zwx, zwy, zwz )
-      !
       IF( MOD( kt, nn_npc ) == 0 ) THEN
-
-         inpcc = 0
-         inpci = 0
-
-         CALL eos( tsa, rhd, zrhop, fsdept_n(:,:,:) )         ! Potential density
-
-         IF( l_trdtra )   THEN                    !* Save ta and sa trends
+         !
+         CALL wrk_alloc( jpi, jpj, jpk, zn2 )    ! N2
+         CALL wrk_alloc( jpi, jpj, jpk, 2, zab ) ! Alpha and Beta
+         CALL wrk_alloc( jpk, 2, zvts, zvab )    ! 1D column vector at point ji,jj
+         CALL wrk_alloc( jpk, zvn2 )             ! 1D column vector at point ji,jj
+
+         IF( l_trdtra )   THEN                    !* Save initial after fields
             CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
             ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
@@ -89,135 +96,240 @@
          ENDIF
 
-         !                                                ! ===============
-         DO jj = 1, jpj                                   !  Vertical slab
-            !                                             ! ===============
-            !  Static instability pointer 
-            ! ----------------------------
-            DO jk = 1, jpkm1
-               DO ji = 1, jpi
-                  zwx(ji,jk) = ( zrhop(ji,jj,jk) - zrhop(ji,jj,jk+1) ) * tmask(ji,jj,jk+1)
-               END DO
-            END DO
-
-            ! 1.1 do not consider the boundary points
-
-            ! even if east-west cyclic b. c. do not considere ji=1 or jpi
-            DO jk = 1, jpkm1
-               zwx( 1 ,jk) = 0.e0
-               zwx(jpi,jk) = 0.e0
-            END DO
-            ! even if south-symmetric b. c. used, do not considere jj=1
-            IF( jj == 1 )   zwx(:,:) = 0.e0
-
-            DO jk = 1, jpkm1
-               DO ji = 1, jpi
-                  zwx(ji,jk) = 1.
-                  IF( zwx(ji,jk) < 1.e-5 ) zwx(ji,jk) = 0.e0
-               END DO
-            END DO
-
-            zwy(:,1) = 0.e0
-            DO ji = 1, jpi
-               DO jk = 1, jpkm1
-                  zwy(ji,1) = zwy(ji,1) + zwx(ji,jk)
-               END DO
-            END DO
-
-            zwz(1,1) = 0.e0
-            DO ji = 1, jpi
-               zwz(1,1) = zwz(1,1) + zwy(ji,1)
-            END DO
-
-            inpcc = inpcc + NINT( zwz(1,1) )
-
-
-            ! 2. Vertical mixing for each instable portion of the density profil
-            ! ------------------------------------------------------------------
-
-            IF( zwz(1,1) /= 0.e0 ) THEN         ! -->> the density profil is statically instable :
-               DO ji = 1, jpi
-                  IF( zwy(ji,1) /= 0.e0 ) THEN
+         !LB debug:
+         IF( lwp .AND. l_LB_debug ) THEN
+            WRITE(numout,*) '' ; WRITE(numout,*) ''
+            WRITE(numout,*) 'LOLO: entering tra_npc, kt, narea =', kt, narea
+         ENDIF
+         !LBdebug: Monitoring of 1 column subject to convection...
+         IF( l_LB_debug ) THEN
+            ! Location of 1 known convection spot to follow what's happening in the water column
+            ilc1 = 54 ;  jlc1 = 15 ; ! Labrador ORCA1 4x4 cpus:
+            nncpu = 15  ; ! the CPU domain contains the convection spot
+            !ilc1 = 14 ;  jlc1 = 13 ; ! Labrador ORCA1 8x8 cpus:
+            !nncpu = 54  ; ! the CPU domain contains the convection spot
+            klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...          
+         ENDIF
+         !LBdebug.
+
+         CALL eos_rab( tsa, zab )         ! after alpha and beta
+         CALL bn2    ( tsa, zab, zn2 )    ! after Brunt-Vaisala
+        
+         inpcc = 0
+
+         DO jj = 2, jpjm1                 ! interior column only
+            DO ji = fs_2, fs_jpim1
+               !
+               IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points
+                  !                                     ! consider one ocean column 
+                  zvts(:,jp_tem) = tsa(ji,jj,:,jp_tem)      ! temperature
+                  zvts(:,jp_sal) = tsa(ji,jj,:,jp_sal)      ! salinity
+
+                  zvab(:,jp_tem)  = zab(ji,jj,:,jp_tem)     ! Alpha 
+                  zvab(:,jp_sal)  = zab(ji,jj,:,jp_sal)     ! Beta  
+                  zvn2(:)         = zn2(ji,jj,:)            ! N^2 
+                 
+                  IF( l_LB_debug ) THEN                  !LB debug:
+                     lp_monitor_point = .FALSE.
+                     IF( ( ji == ilc1 ).AND.( jj == jlc1 ) ) lp_monitor_point = .TRUE.
+                     ! writing only if on CPU domain where conv region is:
+                     lp_monitor_point = (narea == nncpu).AND.lp_monitor_point 
+                     
+                     IF(lp_monitor_point) THEN
+                        WRITE(numout,*) '' ;WRITE(numout,*) '' ;
+                       WRITE(numout,'("Time step = ",i6.6," !!!")')  kt
+                        WRITE(numout,'(" *** BEFORE anything, N^2 for point ",i3,",",i3,":" )') , ji,jj
+                        DO jk = 1, klc1
+                           WRITE(numout,*) jk, zvn2(jk)
+                        END DO
+                        WRITE(numout,*) ' '
+                     ENDIF
+                  ENDIF                                  !LB debug  end
+
+                  ikbot = mbkt(ji,jj)   ! ikbot: ocean bottom T-level
+                  ik = 1                ! because N2 is irrelevant at the surface level (will start at ik=2)
+                  ilayer = 0
+                  jiter  = 0
+                  l_column_treated = .FALSE.
+                 
+                  DO WHILE ( .NOT. l_column_treated )
                      !
-                     ikbot = mbkt(ji,jj)        ! ikbot: ocean bottom T-level
+                     jiter = jiter + 1
+                    
+                     IF( jiter >= 400 ) EXIT
+                    
+                     l_bottom_reached = .FALSE.
+
+                     DO WHILE ( .NOT. l_bottom_reached )
+
+                        ik = ik + 1
+                       
+                        !! Checking level ik for instability
+                        !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+                        IF( zvn2(ik) < 0. ) THEN ! Instability found!
+
+                           ikm  = ik              ! first level whith negative N2
+                           ilayer = ilayer + 1    ! yet another layer found....
+                           IF(jiter == 1) inpcc = inpcc + 1
+
+                           IF(l_LB_debug .AND. lp_monitor_point) &
+                              & WRITE(numout,*) 'Negative N2 at ik =', ikm, ' layer nb.', ilayer, &
+                              & ' inpcc =', inpcc
+
+                           !! Case we mix with upper regions where N2==0:
+                           !! All the points above ikup where N2 == 0 must also be mixed => we go
+                           !! upward to find a new ikup, where the layer doesn't have N2==0
+                           ikup = ikm
+                           DO jk = ikm, 2, -1
+                              ikup = ikup - 1
+                              IF( (zvn2(jk-1) > 0.).OR.(ikup == 1) ) EXIT
+                           END DO
+                          
+                           ! adjusting ikup if the upper part of the unstable column was neutral (N2=0)
+                           IF((zvn2(ikup+1) == 0.).AND.(ikup /= 1)) ikup = ikup+1 ;
+
+                          
+                           IF(lp_monitor_point) WRITE(numout,*) ' => ikup is =', ikup, ' layer nb.', ilayer
+                          
+                           zsum_temp = 0._wp
+                           zsum_sali = 0._wp
+                           zsum_alfa = 0._wp
+                           zsum_beta = 0._wp
+                           zsum_z    = 0._wp
+                                                    
+                           DO jk = ikup, ikbot+1      ! Inside the instable (and overlying neutral) portion of the column
+                              !
+                              IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '     -> summing for jk =', jk
+                              !
+                              zdz       = fse3t(ji,jj,jk)
+                              zsum_temp = zsum_temp + zvts(jk,jp_tem)*zdz
+                              zsum_sali = zsum_sali + zvts(jk,jp_sal)*zdz
+                              zsum_alfa = zsum_alfa + zvab(jk,jp_tem)*zdz
+                              zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz
+                              zsum_z    = zsum_z    + zdz
+                              !
+                              !! EXIT if we found the bottom of the unstable portion of the water column    
+                              IF( (zvn2(jk+1) > 0.).OR.(jk == ikbot ).OR.((jk==ikm).AND.(zvn2(jk+1) == 0.)) )   EXIT
+                           END DO
+                          
+                           !ik     = jk !LB remove?
+                           ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative N2
+                          
+                           IF(l_LB_debug .AND. lp_monitor_point) &
+                              &    WRITE(numout,*) '  => ikdown =', ikdown, '  layer nb.', ilayer
+                          
+                           ! Mixing Temperature and salinity between ikup and ikdown:
+                           zta   = zsum_temp/zsum_z
+                           zsa   = zsum_sali/zsum_z
+                           zalfa = zsum_alfa/zsum_z
+                           zbeta = zsum_beta/zsum_z
+
+                           IF(l_LB_debug .AND. lp_monitor_point) THEN
+                              WRITE(numout,*) '  => Mean temp. in that portion =', zta
+                              WRITE(numout,*) '  => Mean sali. in that portion =', zsa
+                              WRITE(numout,*) '  => Mean Alpha in that portion =', zalfa
+                              WRITE(numout,*) '  => Mean Beta  in that portion =', zbeta
+                           ENDIF
+
+                           !! Homogenaizing the temperature, salinity, alpha and beta in this portion of the column
+                           DO jk = ikup, ikdown
+                              zvts(jk,jp_tem) = zta
+                              zvts(jk,jp_sal) = zsa
+                              zvab(jk,jp_tem) = zalfa
+                              zvab(jk,jp_sal) = zbeta
+                           END DO
+                           !
+                           !! Before updating N2, it is possible that another unstable
+                           !! layer exists underneath the one we just homogeneized!
+                           ik = ikdown
+                           ! 
+                        ENDIF  ! IF( zvn2(ik+1) < 0. ) THEN
+                        !
+                        IF( ik == ikbot ) l_bottom_reached = .TRUE.
+                        !
+                     END DO ! DO WHILE ( .NOT. l_bottom_reached )
+
+                     IF( ik /= ikbot )   STOP 'ERROR: tranpc.F90 => PROBLEM #1'
+                    
+                     ! ******* At this stage ik == ikbot ! *******
+                    
+                     IF( ilayer > 0 ) THEN
+                        !! least an unstable layer has been found
+                        !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion
+                        !! => Need to re-compute N2! will use Alpha and Beta!
+                        !
+                        DO jk = ikup+1, ikdown+1   ! we must go 1 point deeper than ikdown!     
+                           !! Doing exactly as in eosbn2.F90:
+                           !! * Except that we only are interested in the sign of N2 !!!
+                           !!   => just considering the vertical gradient of density
+                           zrw =   (fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk)) &
+                               & / (fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk))
+                           zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw
+                           zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw
+                          
+                           !zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     &
+                           !     &           - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )  &
+                           !     &       / fse3w(ji,jj,jk) * tmask(ji,jj,jk)
+                           zvn2(jk) = (  zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     &
+                                &      - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )  
+                        END DO
+
+                        IF(l_LB_debug .AND. lp_monitor_point) THEN
+                           WRITE(numout, '(" *** After iteration #",i3.3,", N^2 for point ",i3,",",i3,":" )') &
+                              & jiter, ji,jj
+                           DO jk = 1, klc1
+                              WRITE(numout,*) jk, zvn2(jk)
+                           END DO
+                           WRITE(numout,*) ' '
+                        ENDIF
+
+                        ik     = 1  ! starting again at the surface for the next iteration
+                        ilayer = 0
+                     ENDIF
                      !
-                     DO jiter = 1, jpk          ! vertical iteration
-                        !
-                        ! search of ikup : the first static instability from the sea surface
-                        !
-                        ik = 0
-220                     CONTINUE
-                        ik = ik + 1
-                        IF( ik >= ikbot ) GO TO 200
-                        zwx(ji,ik) = zrhop(ji,jj,ik) - zrhop(ji,jj,ik+1)
-                        IF( zwx(ji,ik) <= 0.e0 ) GO TO 220
-                        ikup = ik
-                        ! the density profil is instable below ikup
-                        ! ikdown : bottom of the instable portion of the density profil
-                        ! search of ikdown and vertical mixing from ikup to ikdown
-                        !
-                        ze3tot= fse3t(ji,jj,ikup)
-                        zta   = tsa  (ji,jj,ikup,jp_tem)
-                        zsa   = tsa  (ji,jj,ikup,jp_sal)
-                        zraua = zrhop(ji,jj,ikup)
-                        !
-                        DO jkdown = ikup+1, ikbot-1
-                           IF( zraua <= zrhop(ji,jj,jkdown) ) THEN
-                              ikdown = jkdown
-                              GO TO 240
-                           ENDIF
-                           ze3dwn =  fse3t(ji,jj,jkdown)
-                           ze3tot =  ze3tot + ze3dwn
-                           zta   = ( zta*(ze3tot-ze3dwn) + tsa(ji,jj,jkdown,jp_tem)*ze3dwn )/ze3tot
-                           zsa   = ( zsa*(ze3tot-ze3dwn) + tsa(ji,jj,jkdown,jp_sal)*ze3dwn )/ze3tot
-                           zraua = ( zraua*(ze3tot-ze3dwn) + zrhop(ji,jj,jkdown)*ze3dwn )/ze3tot
-                           inpci = inpci+1
-                        END DO
-                        ikdown = ikbot-1
-240                     CONTINUE
-                        !
-                        DO jkp = ikup, ikdown-1
-                           tsa  (ji,jj,jkp,jp_tem) = zta
-                           tsa  (ji,jj,jkp,jp_sal) = zsa
-                           zrhop(ji,jj,jkp       ) = zraua
-                        END DO
-                        IF (ikdown == ikbot-1 .AND. zraua >= zrhop(ji,jj,ikdown) ) THEN
-                           tsa  (ji,jj,jkp,jp_tem) = zta
-                           tsa  (ji,jj,jkp,jp_sal) = zsa
-                           zrhop(ji,jj,ikdown    ) = zraua
-                        ENDIF
-                     END DO
-                  ENDIF
-200               CONTINUE
-               END DO
-               ! <<-- no more static instability on slab jj
-            ENDIF
-            !                                             ! ===============
-         END DO                                           !   End of slab
-         !                                                ! ===============
-         ! 
-         IF( l_trdtra )   THEN         ! save the Non penetrative mixing trends for diagnostic
-            ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
-            ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
+                     IF( ik >= ikbot ) THEN
+                        IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '    --- exiting jiter loop ---'
+                        l_column_treated = .TRUE.
+                     ENDIF
+                     !
+                  END DO ! DO WHILE ( .NOT. l_column_treated )
+
+                  !! Updating tsa:
+                  tsa(ji,jj,:,jp_tem) = zvts(:,jp_tem)
+                  tsa(ji,jj,:,jp_sal) = zvts(:,jp_sal)
+
+                  !! lolo:  Should we update something else????
+                  !! => like alpha and beta?
+
+                  IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) ''
+
+               ENDIF ! IF( tmask(ji,jj,3) == 1 ) THEN
+
+            END DO ! ji
+         END DO ! jj
+         !
+         IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic
+            z1_r2dt = 1._wp / (2._wp * rdt)
+            ztrdt(:,:,:) = ( tsa(:,:,:,jp_tem) - ztrdt(:,:,:) ) * z1_r2dt
+            ztrds(:,:,:) = ( tsa(:,:,:,jp_sal) - ztrds(:,:,:) ) * z1_r2dt
             CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt )
             CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrds )
             CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
          ENDIF
-      
-         ! Lateral boundary conditions on ( ta, sa )   ( Unchanged sign)
-         ! ------------------------------============
+         !
          CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ;   CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )
-      
-
-         !  2. non penetrative convective scheme statistics
-         !  -----------------------------------------------
-         IF( nn_npcp /= 0 .AND. MOD( kt, nn_npcp ) == 0 ) THEN
-            IF(lwp) WRITE(numout,*)' kt=',kt, ' number of statically instable',   &
-               &                   ' water column : ',inpcc, ' number of iteration : ',inpci
-         ENDIF
-         !
-      ENDIF
-      !
-      CALL wrk_dealloc(jpi, jpj, jpk, zrhop )
-      CALL wrk_dealloc(jpi, jpk, zwx, zwy, zwz )
+         !
+         IF(lwp) THEN
+            WRITE(numout,*) 'LOLO: exiting tra_npc, kt =', kt
+            WRITE(numout,*)' => number of statically instable water column : ',inpcc
+            WRITE(numout,*) '' ; WRITE(numout,*) ''
+         ENDIF
+         !
+         CALL wrk_dealloc(jpi, jpj, jpk, zn2 )
+         CALL wrk_dealloc(jpi, jpj, jpk, 2, zab )
+         CALL wrk_dealloc(jpk, zvn2 )
+         CALL wrk_dealloc(jpk, 2, zvts, zvab )
+         !
+      ENDIF   ! IF( MOD( kt, nn_npc ) == 0 ) THEN
       !
       IF( nn_timing == 1 )  CALL timing_stop('tra_npc')
