Index: branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90	(revision 4750)
+++ branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90	(revision 5965)
@@ -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.6  ! 2015-05  (L. Brodeau) new algorithm based on local Brunt-Vaisala freq.
    !!----------------------------------------------------------------------
 
@@ -14,10 +15,12 @@
    !!   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 trdmod_oce      ! ocean active tracer trends
+   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 +32,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.6 , NEMO Consortium (2014)
+   !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -44,44 +48,47 @@
       !!                  ***  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, ikp, ikup, ikdown, ilayer, ik_low   ! 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), PARAMETER :: zn2_zero = 1.e-14_wp       ! acceptance criteria for neutrality (N2==0)
+      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
+      !
+      LOGICAL, PARAMETER :: l_LB_debug = .FALSE. ! set to true if you want to follow what is
+      INTEGER :: ilc1, jlc1, klc1, nncpu         ! actually happening in a water column at point "ilc1, jlc1"
+      LOGICAL :: lp_monitor_point = .FALSE.      ! in CPU domain "nncpu"
       !!----------------------------------------------------------------------
       !
       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,239 @@
          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
+         IF( l_LB_debug ) THEN
+            ! Location of 1 known convection site to follow what's happening in the water column
+            ilc1 = 45 ;  jlc1 = 3 ; !  ORCA2 4x4, Antarctic coast, more than 2 unstable portions in the  water column...           
+            nncpu = 1  ;            ! the CPU domain contains the convection spot
+            klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...          
+         ENDIF
+         
+         CALL eos_rab( tsa, zab )         ! after alpha and beta (given on T-points)
+         CALL bn2    ( tsa, zab, zn2 )    ! after Brunt-Vaisala  (given on W-points)
+        
+         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                      
+                  ENDIF                                  !LB debug  end
+
+                  ikbot = mbkt(ji,jj)   ! ikbot: ocean bottom T-level
+                  ikp = 1                  ! because N2 is irrelevant at the surface level (will start at ikp=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 )
+
+                        ikp = ikp + 1
+                       
+                        !! Testing level ikp for instability
+                        !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+                        IF( zvn2(ikp) <  -zn2_zero ) THEN ! Instability found!
+
+                           ilayer = ilayer + 1    ! yet another instable portion of the water column found....
+
+                           IF( lp_monitor_point ) THEN 
+                              WRITE(numout,*)
+                              IF( ilayer == 1 .AND. jiter == 1 ) THEN   ! first time a column is spoted with an instability
+                                 WRITE(numout,*)
+                                 WRITE(numout,*) 'Time step = ',kt,' !!!'
+                              ENDIF
+                              WRITE(numout,*)  ' * Iteration #',jiter,': found instable portion #',ilayer,   &
+                                 &                                    ' in column! Starting at ikp =', ikp
+                              WRITE(numout,*)  ' *** N2 for point (i,j) = ',ji,' , ',jj
+                              DO jk = 1, klc1
+                                 WRITE(numout,*) jk, zvn2(jk)
+                              END DO
+                              WRITE(numout,*)
+                           ENDIF
+                           
+
+                           IF( jiter == 1 )   inpcc = inpcc + 1 
+
+                           IF( lp_monitor_point )   WRITE(numout, *) 'Negative N2 at ikp =',ikp,' for layer #', ilayer
+
+                           !! ikup is the uppermost point where mixing will start:
+                           ikup = ikp - 1 ! ikup is always "at most at ikp-1", less if neutral levels overlying
+                           
+                           !! If the points above ikp-1 have N2 == 0 they must also be mixed:
+                           IF( ikp > 2 ) THEN
+                              DO jk = ikp-1, 2, -1
+                                 IF( ABS(zvn2(jk)) < zn2_zero ) THEN
+                                    ikup = ikup - 1  ! 1 more upper level has N2=0 and must be added for the mixing
+                                 ELSE
+                                    EXIT
+                                 ENDIF
+                              END DO
+                           ENDIF
+                           
+                           IF( ikup < 1 )   CALL ctl_stop( 'tra_npc :  PROBLEM #1')
+
+                           zsum_temp = 0._wp
+                           zsum_sali = 0._wp
+                           zsum_alfa = 0._wp
+                           zsum_beta = 0._wp
+                           zsum_z    = 0._wp
+                                                    
+                           DO jk = ikup, ikbot      ! Inside the instable (and overlying neutral) portion of the column
+                              !
+                              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
+                              !                              
+                              IF( jk == ikbot ) EXIT ! avoid array-index overshoot in case ikbot = jpk, cause we're calling jk+1 next line
+                              !! EXIT when we have reached the last layer that is instable (N2<0) or neutral (N2=0):
+                              IF( zvn2(jk+1) > zn2_zero ) EXIT
+                           END DO
+                          
+                           ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative or neutral N2
+                           IF( ikup == ikdown )   CALL ctl_stop( 'tra_npc :  PROBLEM #2')
+
+                           ! Mixing Temperature, salinity, alpha and beta from ikup to ikdown included:
+                           zta   = zsum_temp/zsum_z
+                           zsa   = zsum_sali/zsum_z
+                           zalfa = zsum_alfa/zsum_z
+                           zbeta = zsum_beta/zsum_z
+
+                           IF( lp_monitor_point ) THEN
+                              WRITE(numout,*) 'MIXED T, S, alfa and beta between ikup =',ikup,   &
+                                 &            ' and ikdown =',ikdown,', in layer #',ilayer
+                              WRITE(numout,*) '  => Mean temp. in that portion =', zta
+                              WRITE(numout,*) '  => Mean sali. in that portion =', zsa
+                              WRITE(numout,*) '  => Mean Alfa  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
+                           
+                           
+                           !! Updating N2 in the relvant portion of the water column
+                           !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion
+                           !! => Need to re-compute N2! will use Alpha and Beta!
+                           
+                           ikup   = MAX(2,ikup)         ! ikup can never be 1 !
+                           ik_low = MIN(ikdown+1,ikbot) ! we must go 1 point deeper than ikdown!
+                           
+                           DO jk = ikup, ik_low              ! we must go 1 point deeper than ikdown!
+
+                              !! Interpolating alfa and beta at W point:
+                              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
+
+                              !! N2 at W point, doing exactly as in eosbn2.F90:
+                              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)
+
+                              !! OR, faster  => just considering the vertical gradient of density
+                              !! as only the signa maters...
+                              !zvn2(jk) = (  zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     &
+                              !     &      - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )
+
+                           END DO
+                        
+                           ikp = MIN(ikdown+1,ikbot)
+                           
+
+                        ENDIF  !IF( zvn2(ikp) < 0. )
+
+
+                        IF( ikp == ikbot ) l_bottom_reached = .TRUE.
+                        !
+                     END DO ! DO WHILE ( .NOT. l_bottom_reached )
+
+                     IF( ikp /= ikbot )   CALL ctl_stop( 'tra_npc :  PROBLEM #3')
+                    
+                     ! ******* At this stage ikp == ikbot ! *******
+                    
+                     IF( ilayer > 0 ) THEN      !! least an unstable layer has been found
+                        !
+                        IF( lp_monitor_point ) THEN
+                           WRITE(numout,*)
+                           WRITE(numout,*) 'After ',jiter,' iteration(s), we neutralized ',ilayer,' instable layer(s)'
+                           WRITE(numout,*) '   ==> N2 at i,j=',ji,',',jj,' now looks like this:'
+                           DO jk = 1, klc1
+                              WRITE(numout,*) jk, zvn2(jk)
+                           END DO
+                           WRITE(numout,*)
+                        ENDIF
+                        !
+                        ikp    = 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(:,:,:)
-            CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_npc, ztrdt )
-            CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_npc, ztrds )
+                     IF( ikp >= ikbot )   l_column_treated = .TRUE.
+                     !
+                  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)
+
+                  !! LB:  Potentially some other global variable beside theta and S can be treated here
+                  !!      like BGC tracers.
+
+                  IF( 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
+         !
+         IF( lwp .AND. l_LB_debug ) THEN
+            WRITE(numout,*) 'Exiting tra_npc , kt = ',kt,', => numb. of statically instable water-columns: ', inpcc
+            WRITE(numout,*)
          ENDIF
          !
-      ENDIF
-      !
-      CALL wrk_dealloc(jpi, jpj, jpk, zrhop )
-      CALL wrk_dealloc(jpi, jpk, zwx, zwy, zwz )
+         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')
