Index: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90
===================================================================
--- branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90	(revision 6746)
+++ branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90	(revision 6971)
@@ -31,6 +31,7 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sio3eq   ! chemistry of Si
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   fekeq    ! chemistry of Fe
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   chemc    ! Solubilities of O2 and CO2
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemc    ! Solubilities of O2 and CO2
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemo2   ! Solubilities of O2 and CO2
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tempis   ! In situ temperature
 
    REAL(wp), PUBLIC ::   atcox  = 0.20946         ! units atm
@@ -39,38 +40,9 @@
    REAL(wp) ::   o2atm  = 1. / ( 1000. * 0.20946 )  
 
-   REAL(wp) ::   akcc1  = -171.9065       ! coeff. for apparent solubility equilibrium
-   REAL(wp) ::   akcc2  =   -0.077993     ! Millero et al. 1995 from Mucci 1983
-   REAL(wp) ::   akcc3  = 2839.319        
-   REAL(wp) ::   akcc4  =   71.595        
-   REAL(wp) ::   akcc5  =   -0.77712      
-   REAL(wp) ::   akcc6  =    0.00284263   
-   REAL(wp) ::   akcc7  =  178.34        
-   REAL(wp) ::   akcc8  =   -0.07711     
-   REAL(wp) ::   akcc9  =    0.0041249   
-
-   REAL(wp) ::   rgas   = 83.143         ! universal gas constants
+   REAL(wp) ::   rgas   = 83.14472       ! universal gas constants
    REAL(wp) ::   oxyco  = 1. / 22.4144   ! converts from liters of an ideal gas to moles
 
    REAL(wp) ::   bor1   = 0.00023        ! borat constants
    REAL(wp) ::   bor2   = 1. / 10.82
-
-   REAL(wp) ::   ca0    = -162.8301      ! WEISS & PRICE 1980, units mol/(kg atm)
-   REAL(wp) ::   ca1    =  218.2968
-   REAL(wp) ::   ca2    =   90.9241
-   REAL(wp) ::   ca3    =   -1.47696
-   REAL(wp) ::   ca4    =    0.025695
-   REAL(wp) ::   ca5    =   -0.025225
-   REAL(wp) ::   ca6    =    0.0049867
-
-   REAL(wp) ::   c10    = -3670.7        ! Coeff. for 1. dissoc. of carbonic acid (Edmond and Gieskes, 1970)   
-   REAL(wp) ::   c11    =    62.008     
-   REAL(wp) ::   c12    =    -9.7944    
-   REAL(wp) ::   c13    =     0.0118     
-   REAL(wp) ::   c14    =    -0.000116
-
-   REAL(wp) ::   c20    = -1394.7       ! coeff. for 2. dissoc. of carbonic acid (Millero, 1995)   
-   REAL(wp) ::   c21    =    -4.777   
-   REAL(wp) ::   c22    =     0.0184   
-   REAL(wp) ::   c23    =    -0.000118
 
    REAL(wp) ::   st1    =      0.14     ! constants for calculate concentrations for sulfate
@@ -146,6 +118,6 @@
       REAL(wp) ::   ztgg , ztgg2, ztgg3 , ztgg4 , ztgg5
       REAL(wp) ::   zpres, ztc  , zcl   , zcpexp, zoxy  , zcpexp2
-      REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1
-      REAL(wp) ::   zis  , zis2 , zsal15, zisqrt
+      REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1, zc1, zplat
+      REAL(wp) ::   zis  , zis2 , zsal15, zisqrt, za1  , za2
       REAL(wp) ::   zckb , zck1 , zck2  , zckw  , zak1 , zak2  , zakb , zaksp0, zakw
       REAL(wp) ::   zst  , zft  , zcks  , zckf  , zaksp1
@@ -154,4 +126,20 @@
       IF( nn_timing == 1 )  CALL timing_start('p4z_che')
       !
+      ! Computations of chemical constants require in situ temperature
+      ! Here a quite simple formulation is used to convert 
+      ! potential temperature to in situ temperature. The errors is less than 
+      ! 0.04°C relative to an exact computation
+      ! ---------------------------------------------------------------------
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zpres = fsdept(ji,jj,jk) / 1000.
+               za1 = 0.04 * ( 1.0 + 0.185 * tsn(ji,jj,jk,jp_tem) + 0.035 * (tsn(ji,jj,jk,jp_sal) - 35.0) )
+               za2 = 0.0075 * ( 1.0 - tsn(ji,jj,jk,jp_tem) / 30.0 )
+               tempis(ji,jj,jk) = tsn(ji,jj,jk,jp_tem) - za1 * zpres + za2 * zpres**2
+            END DO
+         END DO
+      END DO
+      !
       ! CHEMICAL CONSTANTS - SURFACE LAYER
       ! ----------------------------------
@@ -161,5 +149,5 @@
          DO ji = 1, jpi
             !                             ! SET ABSOLUTE TEMPERATURE
-            ztkel = tsn(ji,jj,1,jp_tem) + 273.15
+            ztkel = tempis(ji,jj,1) + 273.15
             zt    = ztkel * 0.01
             zt2   = zt * zt
@@ -169,7 +157,10 @@
             !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980)
             !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS
-            zcek1 = ca0 + ca1 / zt + ca2 * zlogt + ca3 * zt2 + zsal * ( ca4 + ca5 * zt + ca6 * zt2 )
+            zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel    &
+            &       + 0.0047036e-4*ztkel**2)
             !                             ! SET SOLUBILITIES OF O2 AND CO2 
-            chemc(ji,jj) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000.  ! mol/(L uatm)
+            chemc(ji,jj,1) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000. ! mol/(kg uatm)
+            chemc(ji,jj,2) = -1636.75 + 12.0408*ztkel - 0.0327957*ztkel**2 + 0.0000316528*ztkel**3
+            chemc(ji,jj,3) = 57.7 - 0.118*ztkel
             !
          END DO
@@ -184,8 +175,8 @@
 !CDIR NOVERRCHK
             DO ji = 1, jpi
-              ztkel = tsn(ji,jj,jk,jp_tem) + 273.15
+              ztkel = tempis(ji,jj,jk) + 273.15
               zsal  = tsn(ji,jj,jk,jp_sal) + ( 1.- tmask(ji,jj,jk) ) * 35.
               zsal2 = zsal * zsal
-              ztgg  = LOG( ( 298.15 - tsn(ji,jj,jk,jp_tem) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature
+              ztgg  = LOG( ( 298.15 - tempis(ji,jj,jk) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature
               ztgg2 = ztgg  * ztgg
               ztgg3 = ztgg2 * ztgg
@@ -210,9 +201,12 @@
             DO ji = 1, jpi
 
-               ! SET PRESSION
-               zpres   = 1.025e-1 * fsdept(ji,jj,jk)
+               ! SET PRESSION ACCORDING TO SAUNDER (1980)
+               zplat   = SIN ( ABS(gphit(ji,jj)*3.141592654/180.) )
+               zc1 = 5.92E-3 + zplat**2 * 5.25E-3
+               zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*fsdept(ji,jj,jk)))) / 4.42E-6
+               zpres = zpres / 10.0
 
                ! SET ABSOLUTE TEMPERATURE
-               ztkel   = tsn(ji,jj,jk,jp_tem) + 273.15
+               ztkel   = tempis(ji,jj,jk) + 273.15
                zsal    = tsn(ji,jj,jk,jp_sal) + ( 1.-tmask(ji,jj,jk) ) * 35.
                zsqrt  = SQRT( zsal )
@@ -223,5 +217,5 @@
                zis2   = zis * zis
                zisqrt = SQRT( zis )
-               ztc     = tsn(ji,jj,jk,jp_tem) + ( 1.- tmask(ji,jj,jk) ) * 20.
+               ztc     = tempis(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20.
 
                ! CHLORINITY (WOOSTER ET AL., 1969)
@@ -256,6 +250,10 @@
 
 
-               zck1    = c10 * ztr + c11 + c12 * zlogt + c13 * zsal + c14 * zsal * zsal
-               zck2    = c20 * ztr + c21 + c22 * zsal   + c23 * zsal**2
+               ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO 
+               ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale
+               zck1    = -1.0*(3633.86*ztr - 61.2172 + 9.6777*zlogt  &
+                  - 0.011555*zsal + 0.0001152*zsal*zsal)
+               zck2    = -1.0*(471.78*ztr + 25.9290 - 3.16967*zlogt      &
+                  - 0.01781*zsal + 0.0001122*zsal*zsal)
 
                ! PKW (H2O) (DICKSON AND RILEY, 1979)
@@ -266,6 +264,7 @@
                ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER
                !       (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983)
-               zaksp0  = akcc1 + akcc2 * ztkel + akcc3 * ztr + akcc4 * LOG10( ztkel )   &
-                  &   + ( akcc5 + akcc6 * ztkel + akcc7 * ztr ) * zsqrt + akcc8 * zsal + akcc9 * zsal15
+               zaksp0  = -171.9065 -0.077993*ztkel + 2839.319*ztr + 71.595*LOG10( ztkel )   &
+                  &      + (-0.77712 + 0.00284263*ztkel + 178.34*ztr) * zsqrt  &
+                  &      - 0.07711*zsal + 0.0041249*zsal15
 
                ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
@@ -337,6 +336,6 @@
       !!                     ***  ROUTINE p4z_che_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj), chemo2(jpi,jpj,jpk),   &
-      &         STAT=p4z_che_alloc )
+      ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk),   &
+      &         tempis(jpi,jpj,jpk), STAT=p4z_che_alloc )
       !
       IF( p4z_che_alloc /= 0 )   CALL ctl_warn('p4z_che_alloc : failed to allocate arrays.')
Index: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90
===================================================================
--- branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90	(revision 6746)
+++ branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90	(revision 6971)
@@ -86,13 +86,14 @@
       REAL(wp) ::   ztc, ztc2, ztc3, ztc4, zws, zkgwan
       REAL(wp) ::   zfld, zflu, zfld16, zflu16, zfact
+      REAL(wp) ::   zvapsw, zsal, zfco2, zxc2, xCO2approx, ztkel, zfugcoeff
       REAL(wp) ::   zph, zah2, zbot, zdic, zalk, zsch_o2, zalka, zsch_co2
       REAL(wp) ::   zyr_dec, zdco2dt
       CHARACTER (len=25) :: charout
-      REAL(wp), POINTER, DIMENSION(:,:) :: zkgco2, zkgo2, zh2co3, zoflx, zw2d 
+      REAL(wp), POINTER, DIMENSION(:,:) :: zkgco2, zkgo2, zh2co3, zoflx, zw2d, zpco2atm 
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('p4z_flx')
       !
-      CALL wrk_alloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx )
+      CALL wrk_alloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx, zpco2atm )
       !
 
@@ -183,15 +184,24 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
+            ztkel  = tsn(ji,jj,1,jp_tem) + 273.15
+            zsal   = tsn(ji,jj,1,jp_sal) + ( 1.- tmask(ji,jj,1) ) * 35.
+            zvapsw = EXP(24.4543 - 67.4509*(100.0/ztkel) - 4.8489*LOG(ztkel/100) - 0.000544*zsal)
+            zpco2atm(ji,jj) = satmco2(ji,jj) * ( patm(ji,jj) - zvapsw )
+            zxc2 = (1.0 - zpco2atm(ji,jj) * 1E-6 )**2
+            zfugcoeff = EXP(patm(ji,jj) * (chemc(ji,jj,2) + 2.0 * zxc2 * chemc(ji,jj,3) )   &
+            &           / (82.05736 * ztkel))
+            zfco2 = zpco2atm(ji,jj) * zfugcoeff
+
             ! Compute CO2 flux for the sea and air
-            zfld = satmco2(ji,jj) * patm(ji,jj) * tmask(ji,jj,1) * chemc(ji,jj) * zkgco2(ji,jj)   ! (mol/L) * (m/s)
-            zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj)                                   ! (mol/L) (m/s) ?
+            zfld = zfco2 * chemc(ji,jj,1) * zkgco2(ji,jj)  ! (mol/L) * (m/s)
+            zflu = zh2co3(ji,jj) * zkgco2(ji,jj)                                   ! (mol/L) (m/s) ?
             oce_co2(ji,jj) = ( zfld - zflu ) * rfact2 * e1e2t(ji,jj) * tmask(ji,jj,1) * 1000.
             ! compute the trend
-            tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) * rfact2 / fse3t(ji,jj,1)
+            tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) * rfact2 / fse3t(ji,jj,1) * tmask(ji,jj,1)
 
             ! Compute O2 flux 
-            zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * tmask(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s)
-            zflu16 = trb(ji,jj,1,jpoxy) * tmask(ji,jj,1) * zkgo2(ji,jj)
-            zoflx(ji,jj) = zfld16 - zflu16
+            zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s)
+            zflu16 = trb(ji,jj,1,jpoxy) * zkgo2(ji,jj)
+            zoflx(ji,jj) = ( zfld16 - zflu16 ) * tmask(ji,jj,1)
             tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + zoflx(ji,jj) * rfact2 / fse3t(ji,jj,1)
          END DO
@@ -224,5 +234,5 @@
          ENDIF
          IF( iom_use( "Dpco2" ) ) THEN
-           zw2d(:,:) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:) + rtrn ) ) * tmask(:,:,1)
+           zw2d(:,:) = ( zpco2atm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1)
            CALL iom_put( "Dpco2" ,  zw2d )
          ENDIF
@@ -240,9 +250,9 @@
             trc2d(:,:,jp_pcs0_2d + 1) = zoflx(:,:) * 1000 * tmask(:,:,1) 
             trc2d(:,:,jp_pcs0_2d + 2) = zkgco2(:,:) * tmask(:,:,1) 
-            trc2d(:,:,jp_pcs0_2d + 3) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:) + rtrn ) ) * tmask(:,:,1) 
-         ENDIF
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx )
+            trc2d(:,:,jp_pcs0_2d + 3) = ( zpco2atm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1)
+         ENDIF
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx, zpco2atm )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_flx')
Index: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90
===================================================================
--- branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 6746)
+++ branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 6971)
@@ -44,4 +44,5 @@
    REAL(wp), PUBLIC ::  xkdoc       !:  2nd half-sat. of DOC remineralization  
    REAL(wp), PUBLIC ::  concbfe     !:  Fe half saturation for bacteria 
+   REAL(wp), PUBLIC ::  oxymin      !:  half saturation constant for anoxia
    REAL(wp), PUBLIC ::  qnfelim     !:  optimal Fe quota for nanophyto
    REAL(wp), PUBLIC ::  qdfelim     !:  optimal Fe quota for diatoms
@@ -121,5 +122,5 @@
                zlim1    = xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk)
                zlim2    = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + concbnh4 )
-               zlim3    = trb(ji,jj,jk,jpfer) / ( concbfe + trb(ji,jj,jk,jpfer) )
+               zlim3    = biron(ji,jj,jk)     / ( concbfe + biron(ji,jj,jk) )
                zlim4    = trb(ji,jj,jk,jpdoc) / ( xkdoc   + trb(ji,jj,jk,jpdoc) )
                xlimbacl(ji,jj,jk) = MIN( zlim1, zlim2, zlim3 )
@@ -187,4 +188,14 @@
       END DO
       !
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ! denitrification factor computed from O2 levels
+               nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - trb(ji,jj,jk,jpoxy) )    &
+                  &                                / ( oxymin + trb(ji,jj,jk,jpoxy) )  )
+               nitrfac(ji,jj,jk) = MIN( 1., nitrfac(ji,jj,jk) )
+            END DO
+         END DO
+      END DO
       !
       IF( lk_iomput .AND. knt == nrdttrc ) THEN        ! save output diagnostics
@@ -216,5 +227,5 @@
       NAMELIST/nampislim/ concnno3, concdno3, concnnh4, concdnh4, concnfer, concdfer, concbfe,   &
          &                concbno3, concbnh4, xsizedia, xsizephy, xsizern, xsizerd,          & 
-         &                xksi1, xksi2, xkdoc, qnfelim, qdfelim, caco3r
+         &                xksi1, xksi2, xkdoc, qnfelim, qdfelim, caco3r, oxymin
       INTEGER :: ios                 ! Local integer output status for namelist read
 
@@ -249,8 +260,11 @@
          WRITE(numout,*) '    Minimum size criteria for nanophyto      xsizephy  = ', xsizephy
          WRITE(numout,*) '    Fe half saturation for bacteria          concbfe   = ', concbfe
+         WRITE(numout,*) '    halk saturation constant for anoxia       oxymin   =' , oxymin
          WRITE(numout,*) '    optimal Fe quota for nano.               qnfelim   = ', qnfelim
          WRITE(numout,*) '    Optimal Fe quota for diatoms             qdfelim   = ', qdfelim
       ENDIF
-
+      !
+      nitrfac (:,:,:) = 0._wp
+      !
    END SUBROUTINE p4z_lim_init
 
Index: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90
===================================================================
--- branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90	(revision 6746)
+++ branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90	(revision 6971)
@@ -65,10 +65,10 @@
       REAL(wp) ::   zomegaca, zexcess, zexcess0
       CHARACTER (len=25) :: charout
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zcaldiss   
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zco3sat, zcaldiss   
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('p4z_lys')
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zco3, zcaldiss )
+      CALL wrk_alloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss )
       !
       zco3    (:,:,:) = 0.
@@ -120,5 +120,6 @@
                zcalcon  = calcon * ( tsn(ji,jj,jk,jp_sal) / 35._wp )
                zfact    = rhop(ji,jj,jk) / 1000._wp
-               zomegaca = ( zcalcon * zco3(ji,jj,jk) * zfact ) / aksp(ji,jj,jk) 
+               zomegaca = ( zcalcon * zco3(ji,jj,jk) ) / ( aksp(ji,jj,jk) * zfact + rtrn )
+               zco3sat(ji,jj,jk) = aksp(ji,jj,jk) * zfact / ( zcalcon + rtrn )
 
                ! SET DEGREE OF UNDER-/SUPERSATURATION
@@ -149,12 +150,12 @@
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
          IF( iom_use( "PH"     ) ) CALL iom_put( "PH"    , -1. * LOG10( hi(:,:,:) )          * tmask(:,:,:) )
-         IF( iom_use( "CO3"    ) ) CALL iom_put( "CO3"   , zco3(:,:,:) * 1.e+3               * tmask(:,:,:) )
-         IF( iom_use( "CO3sat" ) ) CALL iom_put( "CO3sat", aksp(:,:,:) * 1.e+3 / calcon      * tmask(:,:,:) )
-         IF( iom_use( "DCAL"   ) ) CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * 1.e+3 * rfact2r   * tmask(:,:,:) )
+         IF( iom_use( "CO3"    ) ) CALL iom_put( "CO3"   , zco3(:,:,:)    * 1.e+3            * tmask(:,:,:) )
+         IF( iom_use( "CO3sat" ) ) CALL iom_put( "CO3sat", zco3sat(:,:,:) * 1.e+3            * tmask(:,:,:) )
+         IF( iom_use( "DCAL"   ) ) CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * 1.e+3 * rfact2r * tmask(:,:,:) )
       ELSE
          IF( ln_diatrc ) THEN
             trc3d(:,:,:,jp_pcs0_3d    ) = -1. * LOG10( hi(:,:,:) ) * tmask(:,:,:)
             trc3d(:,:,:,jp_pcs0_3d + 1) = zco3(:,:,:)              * tmask(:,:,:)
-            trc3d(:,:,:,jp_pcs0_3d + 2) = aksp(:,:,:) / calcon     * tmask(:,:,:)
+            trc3d(:,:,:,jp_pcs0_3d + 2) = zco3sat(:,:,:)           * tmask(:,:,:)
          ENDIF
       ENDIF
@@ -166,5 +167,5 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zco3, zcaldiss )
+      CALL wrk_dealloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_lys')
Index: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90
===================================================================
--- branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 6746)
+++ branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 6971)
@@ -114,10 +114,8 @@
       !                                        !  --------------------------------------
       IF( l_trcdm2dc ) THEN                     !  diurnal cycle
-         !                                       ! 1% of qsr to compute euphotic layer
-         zqsr100(:,:) = 0.01 * qsr_mean(:,:)     !  daily mean qsr
          !
          zqsr_corr(:,:) = qsr_mean(:,:) / ( 1. - fr_i(:,:) + rtrn )
          !
-         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 ) 
+         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 ) 
          !
          DO jk = 1, nksrp      
@@ -136,10 +134,8 @@
          !
       ELSE
-         ! 1% of qsr to compute euphotic layer
-         zqsr100(:,:) = 0.01 * qsr(:,:)
          !
          zqsr_corr(:,:) = qsr(:,:) / ( 1. - fr_i(:,:) + rtrn )
          !
-         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 ) 
+         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 ) 
          !
          DO jk = 1, nksrp      
@@ -169,5 +165,5 @@
          DO jj = 1, jpj
            DO ji = 1, jpi
-              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.43 * zqsr100(ji,jj) )  THEN
+              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= zqsr100(ji,jj) )  THEN
                  neln(ji,jj) = jk+1                    ! Euphotic level : 1rst T-level strictly below Euphotic layer
                  !                                     ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint
@@ -242,5 +238,5 @@
    END SUBROUTINE p4z_opt
 
-   SUBROUTINE p4z_opt_par( kt, pqsr, pe1, pe2, pe3, pe0 ) 
+   SUBROUTINE p4z_opt_par( kt, pqsr, pe1, pe2, pe3, pe0, pqsr100 ) 
       !!----------------------------------------------------------------------
       !!                  ***  routine p4z_opt_par  ***
@@ -251,8 +247,9 @@
       !!----------------------------------------------------------------------
       !! * arguments
-      INTEGER, INTENT(in)                                       ::  kt            !   ocean time-step
-      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in)              ::  pqsr          !   shortwave
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::  pe1 , pe2 , pe3   !  PAR ( R-G-B)
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL ::  pe0  
+      INTEGER, INTENT(in)                                        ::  kt            !   ocean time-step
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in)               ::  pqsr          !   shortwave
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)            ::  pe1 , pe2 , pe3   !  PAR ( R-G-B)
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL  ::  pe0 
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(out)  , OPTIONAL  ::  pqsr100  
       !! * local variables
       INTEGER    ::   ji, jj, jk     ! dummy loop indices
@@ -264,4 +261,7 @@
       ELSE                  ;  zqsr(:,:) = xparsw         * pqsr(:,:)
       ENDIF
+
+      !  Light at the euphotic depth 
+      IF( PRESENT( pqsr100 ) )  pqsr100(:,:) = 0.01 * 3. * zqsr(:,:)
       !
       IF( PRESENT( pe0 ) ) THEN     !  W-level
Index: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90
===================================================================
--- branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 6746)
+++ branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 6971)
@@ -136,5 +136,5 @@
                   zval = MAX( 1., zstrn(ji,jj) )
                   zval = 1.5 * zval / ( 12. + zval )
-                  zprbio(ji,jj,jk) = prmax(ji,jj,jk) * zval * ( 1. - fr_i(ji,jj) )
+                  zprbio(ji,jj,jk) = prmax(ji,jj,jk) * zval
                   zprdia(ji,jj,jk) = zprbio(ji,jj,jk)
                ENDIF
@@ -202,17 +202,13 @@
                       zdiattot    = ediat(ji,jj,jk) * zstrn(ji,jj)
                       !
-                      zpislopead (ji,jj,jk) = pislope  * ( 1.+ zadap  * EXP( -znanotot ) )
-                      zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp)  / ( trb(ji,jj,jk,jpdia) + rtrn )
-
-                      zpislopen =  zpislopead(ji,jj,jk) * trb(ji,jj,jk,jpnch)                &
-                        &          / ( trb(ji,jj,jk,jpphy) * 12.                  + rtrn )   &
-                        &          / ( prmax(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn )
-
-                      zpislope2n = zpislopead2(ji,jj,jk) * trb(ji,jj,jk,jpdch)                &
-                        &          / ( trb(ji,jj,jk,jpdia) * 12.                  + rtrn )   &
-                        &          / ( prmax(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn )
+                      zpislopead (ji,jj,jk) = pislope  * ( 1.+ zadap  * EXP( -znanotot ) )           &
+                         &                   * trb(ji,jj,jk,jpnch) /( trb(ji,jj,jk,jpphy) * 12. + rtrn)
+                      zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp)  / ( trb(ji,jj,jk,jpdia) + rtrn )   &
+                         &                   * trb(ji,jj,jk,jpdch) /( trb(ji,jj,jk,jpdia) * 12. + rtrn)
 
                       ! Computation of production function for Carbon
                       !  ---------------------------------------------
+                      zpislopen  =  zpislopead(ji,jj,jk)  / ( prmax(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn )
+                      zpislope2n =  zpislopead2(ji,jj,jk) / ( prmax(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn )
                       zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen  * znanotot ) )
                       zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpislope2n * zdiattot ) )
@@ -220,6 +216,6 @@
                       !  Computation of production function for Chlorophyll
                       !--------------------------------------------------
-                      zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen  * enano(ji,jj,jk) ) )
-                      zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) ) )
+                      zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen  * znanotot ) )
+                      zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislope2n * zdiattot ) )
                   ENDIF
                END DO
@@ -227,6 +223,5 @@
          END DO
       ENDIF
-
-
+      
       !  Computation of a proxy of the N/C ratio
       !  ---------------------------------------
@@ -278,10 +273,12 @@
             zmxltst = MAX( 0.e0, hmld(ji,jj) - heup(ji,jj) )
             zmxlday = zmxltst * zmxltst * r1_rday
-            zmixnano(ji,jj) = 1. - zmxlday / ( 2. + zmxlday )
-            zmixdiat(ji,jj) = 1. - zmxlday / ( 4. + zmxlday )
+            zmixnano(ji,jj) = 1. - zmxlday / ( 1. + zmxlday )
+            zmixdiat(ji,jj) = 1. - zmxlday / ( 2. + zmxlday )
          END DO
       END DO
  
-      !  Mixed-layer effect on production                                                                               
+      !  Mixed-layer effect on production 
+      !  Sea-ice effect on production
+
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -291,4 +288,6 @@
                   zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * zmixdiat(ji,jj)
                ENDIF
+                  zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
+                  zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
             END DO
          END DO
@@ -330,57 +329,31 @@
       END DO
 
-      IF( ln_newprod ) THEN
-!CDIR NOVERRCHK
-         DO jk = 1, jpkm1
-!CDIR NOVERRCHK
-            DO jj = 1, jpj
-!CDIR NOVERRCHK
-               DO ji = 1, jpi
-                  IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
-                     zprnch(ji,jj,jk) = zprnch(ji,jj,jk) * zmixnano(ji,jj)
-                     zprdch(ji,jj,jk) = zprdch(ji,jj,jk) * zmixdiat(ji,jj)
-                  ENDIF
-                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-                     !  production terms for nanophyto. ( chlorophyll )
-                     znanotot = enano(ji,jj,jk) * zstrn(ji,jj)
-                     zprod    = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * xlimphy(ji,jj,jk)
-                     zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk)
-                     zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 12. * zprod / &
-                                        & (  zpislopead(ji,jj,jk) * znanotot +rtrn)
-                     !  production terms for diatomees ( chlorophyll )
-                     zdiattot = ediat(ji,jj,jk) * zstrn(ji,jj)
-                     zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * xlimdia(ji,jj,jk)
-                     zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk)
-                     zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 12. * zprod / &
-                                        & ( zpislopead2(ji,jj,jk) * zdiattot +rtrn )
-                  ENDIF
-               END DO
-            END DO
-         END DO
-      ELSE
-!CDIR NOVERRCHK
-         DO jk = 1, jpkm1
-!CDIR NOVERRCHK
-            DO jj = 1, jpj
-!CDIR NOVERRCHK
-               DO ji = 1, jpi
-                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-                     !  production terms for nanophyto. ( chlorophyll )
-                     znanotot = enano(ji,jj,jk)
-                     zprod = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * trb(ji,jj,jk,jpphy) * xlimphy(ji,jj,jk)
-                     zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk)
-                     zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 144. * zprod            &
-                     &                    / ( zpislopead(ji,jj,jk) * trb(ji,jj,jk,jpnch) * znanotot +rtrn )
-                     !  production terms for diatomees ( chlorophyll )
-                     zdiattot = ediat(ji,jj,jk)
-                     zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * trb(ji,jj,jk,jpdia) * xlimdia(ji,jj,jk)
-                     zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk)
-                     zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 144. * zprod             &
-                     &                    / ( zpislopead2(ji,jj,jk) * trb(ji,jj,jk,jpdch) * zdiattot +rtrn )
-                  ENDIF
-               END DO
-            END DO
-         END DO
-      ENDIF
+!CDIR NOVERRCHK
+      DO jk = 1, jpkm1
+!CDIR NOVERRCHK
+         DO jj = 1, jpj
+!CDIR NOVERRCHK
+            DO ji = 1, jpi
+               IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
+                  zprnch(ji,jj,jk) = zprnch(ji,jj,jk) * zmixnano(ji,jj)
+                  zprdch(ji,jj,jk) = zprdch(ji,jj,jk) * zmixdiat(ji,jj)
+               ENDIF
+               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
+                  !  production terms for nanophyto. ( chlorophyll )
+                  znanotot = enano(ji,jj,jk) * zstrn(ji,jj)
+                  zprod    = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * xlimphy(ji,jj,jk)
+                  zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk)
+                  zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 12. * zprod / &
+                                     & (  zpislopead(ji,jj,jk) * znanotot +rtrn)
+                  !  production terms for diatomees ( chlorophyll )
+                  zdiattot = ediat(ji,jj,jk) * zstrn(ji,jj)
+                  zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * xlimdia(ji,jj,jk)
+                  zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk)
+                  zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 12. * zprod / &
+                                     & ( zpislopead2(ji,jj,jk) * zdiattot +rtrn )
+               ENDIF
+            END DO
+         END DO
+      END DO
 
       !   Update the arrays TRA which contain the biological sources and sinks
Index: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90
===================================================================
--- branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 6746)
+++ branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 6971)
@@ -44,6 +44,4 @@
    REAL(wp), PUBLIC ::  xsiremlab  !: fast remineralisation rate of POC 
    REAL(wp), PUBLIC ::  xsilab     !: fraction of labile biogenic silica 
-   REAL(wp), PUBLIC ::  oxymin     !: halk saturation constant for anoxia 
-
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   denitr     !: denitrification array
@@ -111,15 +109,4 @@
                   zdepprod(ji,jj,jk) = zdepmin**0.273
                ENDIF
-            END DO
-         END DO
-      END DO
-
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ! denitrification factor computed from O2 levels
-               nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - trb(ji,jj,jk,jpoxy) )    &
-                  &                                / ( oxymin + trb(ji,jj,jk,jpoxy) )  )
-               nitrfac(ji,jj,jk) = MIN( 1., nitrfac(ji,jj,jk) )
             END DO
          END DO
@@ -357,6 +344,5 @@
       !!
       !!----------------------------------------------------------------------
-      NAMELIST/nampisrem/ xremik, xremip, nitrif, xsirem, xsiremlab, xsilab,   &
-      &                   oxymin
+      NAMELIST/nampisrem/ xremik, xremip, nitrif, xsirem, xsiremlab, xsilab
       INTEGER :: ios                 ! Local integer output status for namelist read
 
@@ -380,8 +366,6 @@
          WRITE(numout,*) '    fraction of labile biogenic silica        xsilab    =', xsilab
          WRITE(numout,*) '    NH4 nitrification rate                    nitrif    =', nitrif
-         WRITE(numout,*) '    halk saturation constant for anoxia       oxymin    =', oxymin
       ENDIF
       !
-      nitrfac (:,:,:) = 0._wp
       denitr  (:,:,:) = 0._wp
       denitnh4(:,:,:) = 0._wp
Index: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90
===================================================================
--- branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 6746)
+++ branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 6971)
@@ -159,10 +159,11 @@
       IF( ln_ndepo ) THEN
          IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_ndep > 1 ) ) THEN
-            CALL fld_read( kt, 1, sf_ndepo )
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  nitdep(ji,jj) = sf_ndepo(1)%fnow(ji,jj,1) / rno3 / ( 14E6 * ryyss * fse3t(ji,jj,1) + rtrn )
-               END DO
-            END DO
+             zcoef = rno3 * 14E6 * ryyss
+             CALL fld_read( kt, 1, sf_ndepo )
+             nitdep(:,:) = sf_ndepo(1)%fnow(:,:,1) / zcoef / fse3t(:,:,1) 
+         ENDIF
+         IF( lk_vvl ) THEN
+           zcoef = rno3 * 14E6 * ryyss
+           nitdep(:,:) = sf_ndepo(1)%fnow(:,:,1) / zcoef / fse3t(:,:,1) 
          ENDIF
       ENDIF
@@ -455,5 +456,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zexpide   = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) )
+                  zexpide   = MIN( 8.,( gdept_0(ji,jj,jk) / 500. )**(-1.5) )
                   zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2
                   zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 )
@@ -465,5 +466,5 @@
          ironsed(:,:,jpk) = 0._wp
          DO jk = 1, jpkm1
-            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( fse3t(:,:,jk) * rday )
+            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( e3t_0(:,:,jk) * rday )
          END DO
          DEALLOCATE( zcmask)
@@ -483,5 +484,7 @@
          CALL iom_close( numhydro )
          !
-         hydrofe(:,:,:) = ( hydrofe(:,:,:) * hratio ) / ( cvol(:,:,:) * ryyss + rtrn ) / 1000._wp
+         DO jk = 1, jpk
+            hydrofe(:,:,jk) = ( hydrofe(:,:,jk) * hratio ) / ( e1e2t(:,:) * e3t_0(:,:,jk) * ryyss + rtrn ) / 1000._wp
+         ENDDO
          !
       ENDIF
