Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90	(revision 5462)
@@ -9,5 +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.
+   !!            3.6  ! 2015-05  (L. Brodeau) new algorithm based on local Brunt-Vaisala freq.
    !!----------------------------------------------------------------------
 
@@ -64,8 +64,9 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   inpcc        ! number of statically instable water column
-      INTEGER  ::   jiter, ikbot, ik, ikup, ikdown, ilayer, ikm   ! local integers
+      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...
@@ -75,9 +76,7 @@
       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 = .FALSE.
-      !!LB debug.
+      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"
       !!----------------------------------------------------------------------
       !
@@ -97,22 +96,13 @@
          ENDIF
 
-         !LB debug:
-         IF( lwp .AND. l_LB_debug ) THEN
-            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
+            ! 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
-         !LBdebug.
-
-         CALL eos_rab( tsa, zab )         ! after alpha and beta
-         CALL bn2    ( tsa, zab, zn2 )    ! after Brunt-Vaisala
+         
+         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
@@ -134,19 +124,9 @@
                      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
+                     lp_monitor_point = (narea == nncpu).AND.lp_monitor_point                      
                   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)
+                  ikp = 1                  ! because N2 is irrelevant at the surface level (will start at ikp=2)
                   ilayer = 0
                   jiter  = 0
@@ -163,34 +143,48 @@
                      DO WHILE ( .NOT. l_bottom_reached )
 
-                        ik = ik + 1
+                        ikp = ikp + 1
                        
-                        !! Checking level ik for instability
+                        !! Testing level ikp 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
-                          
+                        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
@@ -199,7 +193,5 @@
                            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
+                           DO jk = ikup, ikbot      ! Inside the instable (and overlying neutral) portion of the column
                               !
                               zdz       = fse3t(ji,jj,jk)
@@ -209,16 +201,14 @@
                               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
+                              !                              
+                              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
                           
-                           !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:
+                           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
@@ -226,8 +216,10 @@
                            zbeta = zsum_beta/zsum_z
 
-                           IF(l_LB_debug .AND. lp_monitor_point) THEN
+                           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 Alpha in that portion =', zalfa
+                              WRITE(numout,*) '  => Mean Alfa  in that portion =', zalfa
                               WRITE(numout,*) '  => Mean Beta  in that portion =', zbeta
                            ENDIF
@@ -240,57 +232,64 @@
                               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.
+                           
+                           
+                           !! 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( ik /= ikbot )   STOP 'ERROR: tranpc.F90 => PROBLEM #1'
+                     IF( ikp /= ikbot )   CALL ctl_stop( 'tra_npc :  PROBLEM #3')
                     
-                     ! ******* At this stage ik == ikbot ! *******
+                     ! ******* At this stage ikp == 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!
+                     IF( ilayer > 0 ) THEN      !! least an unstable layer has been found
                         !
-                        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
+                        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,*) ' '
+                           WRITE(numout,*)
                         ENDIF
-
-                        ik     = 1  ! starting again at the surface for the next iteration
+                        !
+                        ikp    = 1     ! starting again at the surface for the next iteration
                         ilayer = 0
                      ENDIF
                      !
-                     IF( ik >= ikbot ) THEN
-                        IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '    --- exiting jiter loop ---'
-                        l_column_treated = .TRUE.
-                     ENDIF
+                     IF( ikp >= ikbot )   l_column_treated = .TRUE.
                      !
                   END DO ! DO WHILE ( .NOT. l_column_treated )
@@ -300,8 +299,8 @@
                   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,*) ''
+                  !! 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
@@ -321,8 +320,7 @@
          CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ;   CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )
          !
-         IF(lwp) THEN
-            WRITE(numout,*) 'LOLO: exiting tra_npc, kt =', kt
-            WRITE(numout,*)' => number of statically instable water column : ',inpcc
-            WRITE(numout,*) '' ; WRITE(numout,*) ''
+         IF( lwp .AND. l_LB_debug ) THEN
+            WRITE(numout,*) 'Exiting tra_npc , kt = ',kt,', => numb. of statically instable water-columns: ', inpcc
+            WRITE(numout,*)
          ENDIF
          !
