Index: branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 4819)
+++ branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 4820)
@@ -233,9 +233,9 @@
                   zn3 = EOS013*zt   &
                      &   + EOS103*zs+EOS003
-                  !
+                     !
                   zn2 = (EOS022*zt   &
                      &   + EOS112*zs+EOS012)*zt   &
                      &   + (EOS202*zs+EOS102)*zs+EOS002
-                  !
+                     !
                   zn1 = (((EOS041*zt   &
                      &   + EOS131*zs+EOS031)*zt   &
@@ -243,5 +243,5 @@
                      &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
                      &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
-                  !
+                     !
                   zn0 = (((((EOS060*zt   &
                      &   + EOS150*zs+EOS050)*zt   &
@@ -251,6 +251,6 @@
                      &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
                      &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
-                  !
-                  zn   = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+                     !
+                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
                   !
                   prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm  ! density anomaly (masked)
@@ -273,5 +273,5 @@
                      &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
                      &  - rn_nu * zt * zs
-                  !                                 
+                     !                                 
                   prd(ji,jj,jk) = zn * r1_rau0 * ztm                ! density anomaly (masked)
                END DO
@@ -329,9 +329,9 @@
                   zn3 = EOS013*zt   &
                      &   + EOS103*zs+EOS003
-                  !
+                     !
                   zn2 = (EOS022*zt   &
                      &   + EOS112*zs+EOS012)*zt   &
                      &   + (EOS202*zs+EOS102)*zs+EOS002
-                  !
+                     !
                   zn1 = (((EOS041*zt   &
                      &   + EOS131*zs+EOS031)*zt   &
@@ -339,5 +339,5 @@
                      &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
                      &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
-                  !
+                     !
                   zn0 = (((((EOS060*zt   &
                      &   + EOS150*zs+EOS050)*zt   &
@@ -347,6 +347,6 @@
                      &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
                      &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
-                  !
-                  zn   = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+                     !
+                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
                   !
                   prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface
@@ -426,9 +426,9 @@
                zn3 = EOS013*zt   &
                   &   + EOS103*zs+EOS003
-               !
+                  !
                zn2 = (EOS022*zt   &
                   &   + EOS112*zs+EOS012)*zt   &
                   &   + (EOS202*zs+EOS102)*zs+EOS002
-               !
+                  !
                zn1 = (((EOS041*zt   &
                   &   + EOS131*zs+EOS031)*zt   &
@@ -436,5 +436,5 @@
                   &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
                   &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
-               !
+                  !
                zn0 = (((((EOS060*zt   &
                   &   + EOS150*zs+EOS050)*zt   &
@@ -444,6 +444,6 @@
                   &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
                   &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
-               !
-               zn   = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+                  !
+               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
                !
                prd(ji,jj) = zn * r1_rau0 - 1._wp               ! unmasked in situ density anomaly
@@ -464,5 +464,5 @@
                   &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
                   &  - rn_nu * zt * zs
-               !
+                  !
                prd(ji,jj) = zn * r1_rau0               ! unmasked in situ density anomaly
                !
@@ -515,13 +515,12 @@
                   zn3 = ALP003
                   !
-                  zn2 = (ALP012)*zt   &
-                     &   + ALP102*zs+ALP002
-                  !
-                  zn1 = (((ALP031)*zt   &
+                  zn2 = ALP012*zt + ALP102*zs+ALP002
+                  !
+                  zn1 = ((ALP031*zt   &
                      &   + ALP121*zs+ALP021)*zt   &
                      &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
                      &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
-                  !
-                  zn0 = (((((ALP050)*zt   &
+                     !
+                  zn0 = ((((ALP050*zt   &
                      &   + ALP140*zs+ALP040)*zt   &
                      &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
@@ -529,6 +528,6 @@
                      &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
                      &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
-                  !
-                  zn   = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+                     !
+                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
                   !
                   pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm
@@ -537,13 +536,12 @@
                   zn3 = BET003
                   !
-                  zn2 = (BET012)*zt   &
-                     &   + BET102*zs+BET002
-                  !
-                  zn1 = (((BET031)*zt   &
+                  zn2 = BET012*zt + BET102*zs+BET002
+                  !
+                  zn1 = ((BET031*zt   &
                      &   + BET121*zs+BET021)*zt   &
                      &   + (BET211*zs+BET111)*zs+BET011)*zt   &
                      &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
-                  !
-                  zn0 = (((((BET050)*zt   &
+                     !
+                  zn0 = ((((BET050*zt   &
                      &   + BET140*zs+BET040)*zt   &
                      &   + (BET230*zs+BET130)*zs+BET030)*zt   &
@@ -551,6 +549,6 @@
                      &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
                      &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
-                  !
-                  zn   = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+                     !
+                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
                   !
                   pab(ji,jj,jk,jp_sal) = zn / zs * r1_rau0 * ztm
@@ -627,13 +625,12 @@
                zn3 = ALP003
                !
-               zn2 = (ALP012)*zt   &
-                  &   + ALP102*zs+ALP002
-               !
-               zn1 = (((ALP031)*zt   &
+               zn2 = ALP012*zt + ALP102*zs+ALP002
+               !
+               zn1 = ((ALP031*zt   &
                   &   + ALP121*zs+ALP021)*zt   &
                   &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
                   &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
-               !
-               zn0 = (((((ALP050)*zt   &
+                  !
+               zn0 = ((((ALP050*zt   &
                   &   + ALP140*zs+ALP040)*zt   &
                   &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
@@ -641,6 +638,6 @@
                   &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
                   &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
-               !
-               zn   = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+                  !
+               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
                !
                pab(ji,jj,jp_tem) = zn * r1_rau0
@@ -649,13 +646,12 @@
                zn3 = BET003
                !
-               zn2 = (BET012)*zt   &
-                  &   + BET102*zs+BET002
-               !
-               zn1 = (((BET031)*zt   &
+               zn2 = BET012*zt + BET102*zs+BET002
+               !
+               zn1 = ((BET031*zt   &
                   &   + BET121*zs+BET021)*zt   &
                   &   + (BET211*zs+BET111)*zs+BET011)*zt   &
                   &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
-               !
-               zn0 = (((((BET050)*zt   &
+                  !
+               zn0 = ((((BET050*zt   &
                   &   + BET140*zs+BET040)*zt   &
                   &   + (BET230*zs+BET130)*zs+BET030)*zt   &
@@ -663,6 +659,6 @@
                   &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
                   &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
-               !
-               zn   = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+                  !
+               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
                !
                pab(ji,jj,jp_sal) = zn / zs * r1_rau0
@@ -731,5 +727,5 @@
                zrw =   ( fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk) )   &
                   &  / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) ) 
-               !
+                  !
                zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 
                zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw
@@ -775,30 +771,30 @@
       !
       zdeltaS = 5._wp
-      z1_S0   = 35._wp/35.16504_wp/40._wp
+      z1_S0   = 0.875_wp/35.16504_wp
       z1_T0   = 1._wp/40._wp
       !
       DO jj = 1, jpj
-          DO ji = 1, jpi
-              !
-              zt  = ctmp   (ji,jj) * z1_T0
-              zs  = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 )
-              ztm = tmask(ji,jj,1)
-              !
-              zn = ((((-2.1385727895e-01_wp*zt   &
-                 &   - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt   &
-                 &   + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt   &
-                 &   + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt   &
-                 &   + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs   &
-                 &      +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt   &
-                 &   + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs   &
-                 &      -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp
-              !
-              zd = (2.0035003456_wp*zt   &
-                 &   -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt   &
-                 &   + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp
-              !
-              ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm
-              !
-          END DO
+         DO ji = 1, jpi
+            !
+            zt  = ctmp   (ji,jj) * z1_T0
+            zs  = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 )
+            ztm = tmask(ji,jj,1)
+            !
+            zn = ((((-2.1385727895e-01_wp*zt   &
+               &   - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt   &
+               &   + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt   &
+               &   + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt   &
+               &   + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs   &
+               &      +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt   &
+               &   + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs   &
+               &      -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp
+               !
+            zd = (2.0035003456_wp*zt   &
+               &   -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt   &
+               &   + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp
+               !
+            ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm
+               !
+         END DO
       END DO
       !
@@ -822,9 +818,8 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   psal   ! salinity   [psu]
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ), OPTIONAL ::   pdep   ! depth      [m]
-      ! Leave result array automatic rather than making explicitly allocated
-      REAL(wp), DIMENSION(jpi,jpj) ::   ptf   ! freezing temperature [Celcius]
-      !
-      INTEGER  ::   ji, jj               ! dummy loop indices
-      REAL(wp) ::   zt , zsr             ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj)                          ::   ptf   ! freezing temperature [Celcius]
+      !
+      INTEGER  ::   ji, jj   ! dummy loop indices
+      REAL(wp) ::   zt, zs   ! local scalars
       !!----------------------------------------------------------------------
       !
@@ -835,13 +830,12 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zsr= SQRT( ABS( psal(ji,jj) ) ) * 0.1_wp   ! square root salinity x 0.1
-               ptf(ji,jj) = 0.017947064327968_wp &
-                    &    + zsr*zsr*(-6.07609909992982_wp + zsr*(4.883198653547851_wp &
-                    &    + zsr*(-11.8808160123054_wp + zsr*(13.34658511480257_wp &
-                    &    + zsr*(-8.722761043208607_wp + 2.082038908808201_wp*zsr)))))
+               zs= SQRT( ABS( psal(ji,jj) ) * r1_S0 )           ! square root salinity
+               ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
+                  &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp
             END DO
          END DO
-         !
-         IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.69e-4 * pdep(:,:)
+         ptf(:,:) = ptf(:,:) * psal(:,:)
+         !
+         IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:)
          !
       CASE ( 0 )                     !==  PT,SP (UNESCO formulation)  ==!
@@ -849,5 +843,5 @@
          ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   &
             &                     - 2.154996e-4_wp *       psal(:,:)   ) * psal(:,:)
-         !
+            !
          IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:)
          !
@@ -909,9 +903,9 @@
                   zn2 = (PEN012)*zt   &
                      &   + PEN102*zs+PEN002
-                  !
+                     !
                   zn1 = ((PEN021)*zt   &
                      &   + PEN111*zs+PEN011)*zt   &
                      &   + (PEN201*zs+PEN101)*zs+PEN001
-                  !
+                     !
                   zn0 = ((((PEN040)*zt   &
                      &   + PEN130*zs+PEN030)*zt   &
@@ -919,6 +913,6 @@
                      &   + ((PEN310*zs+PEN210)*zs+PEN110)*zs+PEN010)*zt   &
                      &   + (((PEN400*zs+PEN300)*zs+PEN200)*zs+PEN100)*zs+PEN000
-                  !
-                  zn   = ( zn2 * zh + zn1 ) * zh + zn0
+                     !
+                  zn  = ( zn2 * zh + zn1 ) * zh + zn0
                   !
                   ppen(ji,jj,jk)  = zn * zh * r1_rau0 * ztm
@@ -929,11 +923,11 @@
                   zn1 = (APE011)*zt   &
                      &   + APE101*zs+APE001
-                  !
+                     !
                   zn0 = (((APE030)*zt   &
                      &   + APE120*zs+APE020)*zt   &
                      &   + (APE210*zs+APE110)*zs+APE010)*zt   &
                      &   + ((APE300*zs+APE200)*zs+APE100)*zs+APE000
-                  !
-                  zn   = ( zn2 * zh + zn1 ) * zh + zn0
+                     !
+                  zn  = ( zn2 * zh + zn1 ) * zh + zn0
                   !                              
                   pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rau0 * ztm
@@ -944,11 +938,11 @@
                   zn1 = (BPE011)*zt   &
                      &   + BPE101*zs+BPE001
-                  !
+                     !
                   zn0 = (((BPE030)*zt   &
                      &   + BPE120*zs+BPE020)*zt   &
                      &   + (BPE210*zs+BPE110)*zs+BPE010)*zt   &
                      &   + ((BPE300*zs+BPE200)*zs+BPE100)*zs+BPE000
-                  !
-                  zn   = ( zn2 * zh + zn1 ) * zh + zn0
+                     !
+                  zn  = ( zn2 * zh + zn1 ) * zh + zn0
                   !                              
                   pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rau0 * ztm
@@ -1013,5 +1007,5 @@
       WRITE( numond, nameos )
       !
-      rau0        = 1026._wp                  !: volumic mass of reference     [kg/m3]
+      rau0        = 1026._wp                 !: volumic mass of reference     [kg/m3]
       rcp         = 3991.86795711963_wp      !: heat capacity     [J/K]
       !
@@ -1035,5 +1029,5 @@
          !
          rdeltaS = 32._wp
-         r1_S0  = 35._wp/35.16504_wp/40._wp
+         r1_S0  = 0.875_wp/35.16504_wp
          r1_T0  = 1._wp/40._wp
          r1_Z0  = 1.e-4_wp
@@ -1128,39 +1122,39 @@
          ALP003 = -9.4924551138e-03_wp
          !
-         ALP000 = -6.5025362670e-01_wp
-         ALP100 = 1.6320471316_wp
-         ALP200 = -2.0442606277_wp
-         ALP300 = 1.4222011580_wp
-         ALP400 = -4.4204535284e-01_wp
-         ALP500 = 4.7983755487e-02_wp
-         ALP010 = 1.8537085209_wp
-         ALP110 = -3.0774129064_wp
-         ALP210 = 3.0181275751_wp
-         ALP310 = -1.4565010626_wp
-         ALP410 = 2.7361846370e-01_wp
-         ALP020 = -1.6246342147_wp
-         ALP120 = 2.5086831352_wp
-         ALP220 = -1.4787808849_wp
-         ALP320 = 2.3807209899e-01_wp
-         ALP030 = 8.3627885467e-01_wp
-         ALP130 = -1.1311538584_wp
-         ALP230 = 5.3563304045e-01_wp
-         ALP040 = -6.7560904739e-02_wp
-         ALP140 = -6.0212475204e-02_wp
-         ALP050 = 2.8625353333e-02_wp
-         ALP001 = 3.3340752782e-01_wp
-         ALP101 = 1.1217528644e-01_wp
-         ALP201 = -1.2510649515e-01_wp
-         ALP301 = 1.6349760916e-02_wp
-         ALP011 = -3.3540239802e-01_wp
-         ALP111 = -1.7531540640e-01_wp
-         ALP211 = 9.3976864981e-02_wp
-         ALP021 = 1.8487252150e-01_wp
-         ALP121 = 4.1307825959e-02_wp
-         ALP031 = -5.5927935970e-02_wp
-         ALP002 = -5.1410778748e-02_wp
-         ALP102 = 5.3278413794e-03_wp
-         ALP012 = 6.2099915132e-02_wp
-         ALP003 = -9.4924551138e-03_wp
+         BET000 = 1.0783203594e+01_wp
+         BET100 = -4.4452095908e+01_wp
+         BET200 = 7.6048755820e+01_wp
+         BET300 = -6.3944280668e+01_wp
+         BET400 = 2.6890441098e+01_wp
+         BET500 = -4.5221697773_wp
+         BET010 = -8.1219372432e-01_wp
+         BET110 = 2.0346663041_wp
+         BET210 = -2.1232895170_wp
+         BET310 = 8.7994140485e-01_wp
+         BET410 = -1.1939638360e-01_wp
+         BET020 = 7.6574242289e-01_wp
+         BET120 = -1.5019813020_wp
+         BET220 = 1.0872489522_wp
+         BET320 = -2.7233429080e-01_wp
+         BET030 = -4.1615152308e-01_wp
+         BET130 = 4.9061350869e-01_wp
+         BET230 = -1.1847737788e-01_wp
+         BET040 = 1.4073062708e-01_wp
+         BET140 = -1.3327978879e-01_wp
+         BET050 = 5.9929880134e-03_wp
+         BET001 = -5.2937873009e-01_wp
+         BET101 = 1.2634116779_wp
+         BET201 = -1.1547328025_wp
+         BET301 = 3.2870876279e-01_wp
+         BET011 = -5.5824407214e-02_wp
+         BET111 = 1.2451933313e-01_wp
+         BET211 = -2.4409539932e-02_wp
+         BET021 = 4.3623149752e-02_wp
+         BET121 = -4.6767901790e-02_wp
+         BET031 = -6.8523260060e-03_wp
+         BET002 = -6.1618945251e-02_wp
+         BET102 = 6.2255521644e-02_wp
+         BET012 = -2.6514181169e-03_wp
+         BET003 = -2.3025968587e-04_wp
          !
          PEN000 = -9.8409626043_wp
Index: branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90
===================================================================
--- branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 4819)
+++ branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 4820)
@@ -19,5 +19,5 @@
    USE zdf_oce        ! ocean vertical physics
    USE trd_oce        ! trends: ocean variables
-   USE trdtrc     ! ocean passive mixed layer tracers trends 
+   USE trdtrc         ! ocean passive mixed layer tracers trends 
    USE trdglo         ! trends: global domain averaged
    USE trdpen         ! trends: Potential ENergy
Index: branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
===================================================================
--- branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 4819)
+++ branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 4820)
@@ -80,5 +80,5 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   iikn, iiki   ! local integer
-      REAL(wp) ::   zmult        ! local scalar
+      REAL(wp) ::   zN2_c        ! local scalar
       INTEGER, POINTER, DIMENSION(:,:) ::   imld   ! 2D workspace
       !!----------------------------------------------------------------------
@@ -97,28 +97,22 @@
 
       ! w-level of the mixing and mixed layers
-      nmln     (:,:)  = nlb10                ! Initialization to the number of w ocean point
-      hmlp     (:,:)  = 0._wp                ! used as a dummy variable, integrating vertically N2
-      zmult = rho_c * r1_rau0 * grav
+      nmln(:,:)  = nlb10               ! Initialization to the number of w ocean point
+      hmlp(:,:)  = 0._wp               ! here hmlp used as a dummy variable, integrating vertically N^2
+      zN2_c = grav * rho_c * r1_rau0   ! convert density criteria into N^2 criteria
       DO jk = nlb10, jpkm1
-         DO jj = 1, jpj
+         DO jj = 1, jpj                ! Mixed layer level: w-level 
             DO ji = 1, jpi
-               hmlp(ji,jj) = hmlp(ji,jj) + rn2b(ji,jj,jk) * fse3w(ji,jj,jk)
-               IF( hmlp(ji,jj) < zmult )   nmln(ji,jj) = jk + 1    ! Mixed layer
+               hmlp(ji,jj) = hmlp(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * fse3w(ji,jj,jk)
+               IF( hmlp(ji,jj) < zN2_c )   nmln(ji,jj) = MIN( jk , mbkt(ji,jj) ) + 1   ! Mixed layer level
             END DO
          END DO
       END DO
       !
-      DO jj = 1, jpj                      ! bottom k-index of u- (v-) level
-         DO ji = 1, jpi
-            nmln(ji,jj) = MIN(  nmln(ji,jj) , mbkt(ji,jj) + 1  )
-         END DO
-      END DO
-      !! end fabien
       ! w-level of the turbocline
       imld(:,:) = mbkt(:,:) + 1        ! Initialization to the number of w ocean point
-      DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10
+      DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10 
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( avt (ji,jj,jk) < avt_c                     )   imld(ji,jj) = jk      ! Turbocline 
+               IF( avt (ji,jj,jk) < avt_c )   imld(ji,jj) = jk      ! Turbocline 
             END DO
          END DO
Index: branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90
===================================================================
--- branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 4819)
+++ branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 4820)
@@ -11,5 +11,5 @@
    !!   'key_pisces'                                       PISCES bio-model
    !!----------------------------------------------------------------------
-   !!   p4zsms        :  Time loop of passive tracers sms
+   !!   p4zsms         :  Time loop of passive tracers sms
    !!----------------------------------------------------------------------
    USE oce_trc         !  shared variables between ocean and passive tracers
@@ -25,6 +25,6 @@
    USE p4zint          !  time interpolation
    USE iom             !  I/O manager
-   USE trd_oce      !  Ocean trends variables
-   USE trdtrc      !  TOP trends variables
+   USE trd_oce         !  Ocean trends variables
+   USE trdtrc          !  TOP trends variables
    USE sedmodel        !  Sediment model
    USE prtctl_trc      !  print control for debugging
@@ -33,6 +33,6 @@
    PRIVATE
 
-   PUBLIC   p4z_sms_init    ! called in p4zsms.F90
-   PUBLIC   p4z_sms    ! called in p4zsms.F90
+   PUBLIC   p4z_sms_init   ! called in p4zsms.F90
+   PUBLIC   p4z_sms        ! called in p4zsms.F90
 
    REAL(wp) :: alkbudget, no3budget, silbudget, ferbudget
