Index: /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 8532)
+++ /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 8533)
@@ -207,4 +207,9 @@
                   &                                / ( oxymin + trb(ji,jj,jk,jpoxy) )  )
                nitrfac(ji,jj,jk) = MIN( 1., nitrfac(ji,jj,jk) )
+               !
+               ! denitrification factor computed from NO3 levels
+               nitrfac2(ji,jj,jk) = MAX( 0.e0,       ( 1.E-6 - trb(ji,jj,jk,jpno3) )  &
+                  &                                / ( 1.E-6 + trb(ji,jj,jk,jpno3) ) )
+               nitrfac2(ji,jj,jk) = MIN( 1., nitrfac2(ji,jj,jk) )
             END DO
          END DO
Index: /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90	(revision 8532)
+++ /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90	(revision 8533)
@@ -95,5 +95,5 @@
                !  no real reason except that it seems to be more stable and may mimic predation
                !  ---------------------------------------------------------------
-               ztortz2   = mzrat2 * 1.e6 * zfact * trb(ji,jj,jk,jpmes)
+               ztortz2   = mzrat2 * 1.e6 * zfact * trb(ji,jj,jk,jpmes)  * (1. - nitrfac(ji,jj,jk) )
                !
                zcompadi  = MAX( ( trb(ji,jj,jk,jpdia) - xthresh2dia ), 0.e0 )
@@ -125,8 +125,10 @@
                !  ----------------------------------
                zgrazffeg = grazflux  * xstep * wsbio4(ji,jj,jk)      &
-               &           * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpgoc) * trb(ji,jj,jk,jpmes)
+               &           * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpgoc) * trb(ji,jj,jk,jpmes) &
+               &           * (1. - nitrfac(ji,jj,jk))
                zgrazfffg = zgrazffeg * trb(ji,jj,jk,jpbfe) / (trb(ji,jj,jk,jpgoc) + rtrn)
                zgrazffep = grazflux  * xstep *  wsbio3(ji,jj,jk)     &
-               &           * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpmes)
+               &           * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpmes) &
+               &           * (1. - nitrfac(ji,jj,jk))
                zgrazfffp = zgrazffep * trb(ji,jj,jk,jpsfe) / (trb(ji,jj,jk,jppoc) + rtrn)
               !
Index: /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90	(revision 8532)
+++ /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90	(revision 8533)
@@ -93,5 +93,5 @@
                !  no real reason except that it seems to be more stable and may mimic predation.
                !  ---------------------------------------------------------------
-               ztortz = mzrat * 1.e6 * zfact * trb(ji,jj,jk,jpzoo)
+               ztortz = mzrat * 1.e6 * zfact * trb(ji,jj,jk,jpzoo) * (1. - nitrfac(ji,jj,jk))
 
                zcompadi  = MIN( MAX( ( trb(ji,jj,jk,jpdia) - xthreshdia ), 0.e0 ), xsizedia )
@@ -105,5 +105,5 @@
                zdenom    = zfoodlim / ( xkgraz + zfoodlim )
                zdenom2   = zdenom / ( zfood + rtrn )
-               zgraze    = grazrat * xstep * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpzoo) 
+               zgraze    = grazrat * xstep * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpzoo) * (1. - nitrfac(ji,jj,jk))
 
                zgrazp    = zgraze  * xpref2p * zcompaph  * zdenom2 
Index: /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 8532)
+++ /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 8533)
@@ -65,4 +65,5 @@
       REAL(wp) ::   zsatur, zsatur2, znusil, znusil2, zdep, zdepmin, zfactdep
       REAL(wp) ::   zbactfer, zolimit, zonitr, zrfact2
+      REAL(wp) ::   zammonic, zoxyrem
       REAL(wp) ::   zosil, ztem, zdenitnh4, zolimic, zolimin, zolimip, zdenitrn, zdenitrp
       CHARACTER (len=25) :: charout
@@ -118,17 +119,20 @@
                   ! Ammonification in suboxic waters with denitrification
                   ! -------------------------------------------------------
-                  denitr(ji,jj,jk)  = MIN(  ( trb(ji,jj,jk,jpno3) - rtrn ) / rdenit,   &
-                     &                     zremik * nitrfac(ji,jj,jk) * trb(ji,jj,jk,jpdoc)  )
+                  zammonic = zremik * nitrfac(ji,jj,jk) * trb(ji,jj,jk,jpdoc)
+                  denitr(ji,jj,jk)  = zammonic * ( 1. - nitrfac2(ji,jj,jk) )
+                  zoxyrem           = zammonic *        nitrfac2(ji,jj,jk)
                   !
                   zolimi (ji,jj,jk) = MAX( 0.e0, zolimi (ji,jj,jk) )
                   denitr (ji,jj,jk) = MAX( 0.e0, denitr (ji,jj,jk) )
+                  zoxyrem           = MAX( 0.e0, zoxyrem  )
+
                   !
-                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zolimi (ji,jj,jk) + denitr(ji,jj,jk)
-                  tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zolimi (ji,jj,jk) + denitr(ji,jj,jk)
+                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zolimi (ji,jj,jk) + denitr(ji,jj,jk) + zoxyrem
+                  tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zolimi (ji,jj,jk) + denitr(ji,jj,jk) + zoxyrem
                   tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) - denitr (ji,jj,jk) * rdenit
-                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zolimi (ji,jj,jk) - denitr(ji,jj,jk)
+                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zolimi (ji,jj,jk) - denitr(ji,jj,jk) - zoxyrem
                   tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - zolimi (ji,jj,jk) * o2ut
-                  tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zolimi (ji,jj,jk) + denitr(ji,jj,jk)
-                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * ( zolimi(ji,jj,jk)    &
+                  tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zolimi (ji,jj,jk) + denitr(ji,jj,jk) + zoxyrem
+                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * ( zolimi(ji,jj,jk) + zoxyrem    &
                   &                     + ( rdenit + 1.) * denitr(ji,jj,jk) )
                END DO
Index: /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 8532)
+++ /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 8533)
@@ -54,5 +54,5 @@
       REAL(wp) ::  zwflux, zfminus, zfplus
       REAL(wp) ::  zlim, zfact, zfactcal
-      REAL(wp) ::  zo2, zno3, zflx, zpdenit, z1pdenit, zdenitt, zolimit
+      REAL(wp) ::  zo2, zno3, zflx, zpdenit, z1pdenit, zolimit
       REAL(wp) ::  zsiloss, zcaloss, zws3, zws4, zwsc, zdep
       REAL(wp) ::  zwstpoc, zwstpon, zwstpop
@@ -319,6 +319,6 @@
                tra(ji,jj,ikt,jptal) =  tra(ji,jj,ikt,jptal) + zcaloss * zrivalk * 2.0
                tra(ji,jj,ikt,jpdic) =  tra(ji,jj,ikt,jpdic) + zcaloss * zrivalk
-               zsedcal(ji,jj) = (1.0 - zrivalk) * zcaloss / zdep
-               zsedsi (ji,jj) = (1.0 - zrivsil) * zsiloss / zdep
+               zsedcal(ji,jj) = (1.0 - zrivalk) * zcaloss * e3t_n(ji,jj,ikt) 
+               zsedsi (ji,jj) = (1.0 - zrivsil) * zsiloss * e3t_n(ji,jj,ikt) 
             END DO
          END DO
@@ -365,6 +365,6 @@
 
       IF( .NOT.lk_sed ) THEN
-         ! The 0.5 factor in zpdenit and zdenitt is to avoid negative NO3 concentration after both denitrification
-         ! in the sediments and just above the sediments. Not very clever, but simpliest option.
+         ! The 0.5 factor in zpdenit is to avoid negative NO3 concentration after
+         ! denitrification in the sediments. Not very clever, but simpliest option.
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -378,14 +378,13 @@
                z1pdenit = zwstpoc * zrivno3 - zpdenit
                zolimit = MIN( ( trb(ji,jj,ikt,jpoxy) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) )
-               zdenitt = MIN(  0.5 * ( trb(ji,jj,ikt,jpno3) - rtrn ) / rdenit, z1pdenit * nitrfac(ji,jj,ikt) )
-               tra(ji,jj,ikt,jpdoc) = tra(ji,jj,ikt,jpdoc) + z1pdenit - zolimit - zdenitt
-               tra(ji,jj,ikt,jppo4) = tra(ji,jj,ikt,jppo4) + zpdenit + zolimit + zdenitt
-               tra(ji,jj,ikt,jpnh4) = tra(ji,jj,ikt,jpnh4) + zpdenit + zolimit + zdenitt
-               tra(ji,jj,ikt,jpno3) = tra(ji,jj,ikt,jpno3) - rdenit * (zpdenit + zdenitt)
+               tra(ji,jj,ikt,jpdoc) = tra(ji,jj,ikt,jpdoc) + z1pdenit - zolimit
+               tra(ji,jj,ikt,jppo4) = tra(ji,jj,ikt,jppo4) + zpdenit + zolimit
+               tra(ji,jj,ikt,jpnh4) = tra(ji,jj,ikt,jpnh4) + zpdenit + zolimit
+               tra(ji,jj,ikt,jpno3) = tra(ji,jj,ikt,jpno3) - rdenit * zpdenit
                tra(ji,jj,ikt,jpoxy) = tra(ji,jj,ikt,jpoxy) - zolimit * o2ut
-               tra(ji,jj,ikt,jptal) = tra(ji,jj,ikt,jptal) + rno3 * (zolimit + (1.+rdenit) * (zpdenit + zdenitt) )
-               tra(ji,jj,ikt,jpdic) = tra(ji,jj,ikt,jpdic) + zpdenit + zolimit + zdenitt
+               tra(ji,jj,ikt,jptal) = tra(ji,jj,ikt,jptal) + rno3 * (zolimit + (1.+rdenit) * zpdenit )
+               tra(ji,jj,ikt,jpdic) = tra(ji,jj,ikt,jpdic) + zpdenit + zolimit 
                sdenit(ji,jj) = rdenit * zpdenit * e3t_n(ji,jj,ikt)
-               zsedc(ji,jj)   = (1. - zrivno3) * zwstpoc / zdep
+               zsedc(ji,jj)   = (1. - zrivno3) * zwstpoc * e3t_n(ji,jj,ikt)
                IF( ln_p5z ) THEN
                   zwstpop              = trb(ji,jj,ikt,jpgop) * zws4 + trb(ji,jj,ikt,jppop) * zws3
@@ -494,17 +493,17 @@
       IF( lk_iomput ) THEN
          IF( knt == nrdttrc ) THEN
-            zfact = 1.e+3 * rfact2r * rno3  !  conversion from molC/l/kt  to molN/m3/s
-            IF( iom_use("Nfix"   ) ) CALL iom_put( "Nfix", nitrpot(:,:,:) * nitrfix * zfact * tmask(:,:,:) )  ! nitrogen fixation 
+            zfact = 1.e+3 * rfact2r !  conversion from molC/l/kt  to molN/m3/s
+            IF( iom_use("Nfix"   ) ) CALL iom_put( "Nfix", nitrpot(:,:,:) * nitrfix * rno3 * zfact * tmask(:,:,:) )  ! nitrogen fixation 
             IF( iom_use("INTNFIX") ) THEN   ! nitrogen fixation rate in ocean ( vertically integrated )
                zwork1(:,:) = 0.
                DO jk = 1, jpkm1
-                 zwork1(:,:) = zwork1(:,:) + nitrpot(:,:,jk) * nitrfix * zfact * e3t_n(:,:,jk) * tmask(:,:,jk)
+                 zwork1(:,:) = zwork1(:,:) + nitrpot(:,:,jk) * nitrfix * rno3 * zfact * e3t_n(:,:,jk) * tmask(:,:,jk)
                ENDDO
                CALL iom_put( "INTNFIX" , zwork1 ) 
             ENDIF
-            IF( iom_use("SedCal" ) ) CALL iom_put( "SedCal", zsedcal(:,:) * 1.e+3 )
-            IF( iom_use("SedSi" ) )  CALL iom_put( "SedSi",  zsedsi (:,:) * 1.e+3 )
-            IF( iom_use("SedC" ) )   CALL iom_put( "SedC",   zsedc  (:,:) * 1.e+3 )
-            IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", sdenit (:,:) * 1.e+3 * rno3 )
+            IF( iom_use("SedCal" ) ) CALL iom_put( "SedCal", zsedcal(:,:) * zfact )
+            IF( iom_use("SedSi" ) )  CALL iom_put( "SedSi",  zsedsi (:,:) * zfact )
+            IF( iom_use("SedC" ) )   CALL iom_put( "SedC",   zsedc  (:,:) * zfact )
+            IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", sdenit (:,:) * zfact * rno3 )
          ENDIF
       ENDIF
Index: /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 8532)
+++ /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 8533)
@@ -431,11 +431,11 @@
 
       IF( kt == nittrc000 ) THEN 
+         xfact1 = rfact2r * 12. / 1.e15 * ryyss    ! conversion molC/kt --> PgC/yr
+         xfact2 = 1.e+3 * rno3 * 14. / 1.e12 * ryyss   ! conversion molC/l/s ----> TgN/m3/yr
+         xfact3 = 1.e+3 * rfact2r * rno3   ! conversion molC/l/kt ----> molN/m3/s
          IF( ln_check_mass .AND. lwp) THEN      !   Open budget file of NO3, ALK, Si, Fer
             CALL ctl_opn( numco2, 'carbon.budget'  , 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
             CALL ctl_opn( numnut, 'nutrient.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
             CALL ctl_opn( numnit, 'nitrogen.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
-            xfact1 = rfact2r * 12. / 1.e15 * ryyss    ! conversion molC/kt --> PgC/yr
-            xfact2 = 1.e+3 * rno3 * 14. / 1.e12 * ryyss   ! conversion molC/l/s ----> TgN/m3/yr
-            xfact3 = 1.e+3 * rfact2r * rno3   ! conversion molC/l/kt ----> molN/m3/s
             cltxt='time-step   Alkalinity        Nitrate        Phosphorus         Silicate           Iron'
             IF( lwp ) WRITE(numnut,*)  TRIM(cltxt)
@@ -517,17 +517,13 @@
       IF( iom_use( "tnfix" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
          znitrpottot  = glob_sum ( nitrpot(:,:,:) * nitrfix * cvol(:,:,:) )
-         CALL iom_put( "tnfix"  , znitrpottot * 1.e+3 * rno3 )  ! Global  nitrogen fixation molC/l  to molN/m3 
+         CALL iom_put( "tnfix"  , znitrpottot * xfact3 )  ! Global  nitrogen fixation molC/l  to molN/m3 
       ENDIF
       !
       IF( iom_use( "tdenit" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
-         zrdenittot   = glob_sum ( denitr(:,:,:) * rdenit * xnegtr(:,:,:) * cvol(:,:,:) )
-         CALL iom_put( "tdenit"  , zrdenittot * 1.e+3 * rno3 )  ! Total denitrification molC/l to molN/m3 
-      ENDIF
-      !
-      IF( iom_use( "Sdenit" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
-         zsdenittot   = glob_sum ( sdenit(:,:) * e1e2t(:,:) )
-         CALL iom_put( "Sdenit", sdenit(:,:) * xfact3 * tmask(:,:,1) )  ! Nitrate reduction in the sediments
-      ENDIF
-
+         zrdenittot = glob_sum ( denitr(:,:,:) * rdenit * xnegtr(:,:,:) * cvol(:,:,:) )
+         zsdenittot = glob_sum ( sdenit(:,:) * e1e2t(:,:) * tmask(:,:,1) )
+         CALL iom_put( "tdenit" , ( zrdenittot + zsdenittot ) * xfact3 )  ! Total denitrification molC/l to molN/m3 
+      ENDIF
+      !
       IF( ln_check_mass .AND. kt == nitend ) THEN   ! Compute the budget of NO3, ALK, Si, Fer
          t_atm_co2_flx  = t_atm_co2_flx / glob_sum( e1e2t(:,:) )
Index: /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90	(revision 8532)
+++ /trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90	(revision 8533)
@@ -97,4 +97,5 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xfracal    !: ??
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   nitrfac    !: ??
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   nitrfac2   !: ??
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   orem       !: ??
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xdiss      !: ??
@@ -159,7 +160,7 @@
          !
          !*  SMS for the organic matter
-         ALLOCATE( xfracal (jpi,jpj,jpk), nitrfac(jpi,jpj,jpk) ,    &
-            &      orem    (jpi,jpj,jpk),                           &
-            &      prodcal(jpi,jpj,jpk),  xdiss   (jpi,jpj,jpk),    &
+         ALLOCATE( xfracal (jpi,jpj,jpk), orem(jpi,jpj,jpk)    ,    &
+            &      nitrfac(jpi,jpj,jpk), nitrfac2(jpi,jpj,jpk) ,    &
+            &      prodcal(jpi,jpj,jpk) , xdiss   (jpi,jpj,jpk),    &
             &      prodpoc(jpi,jpj,jpk) , conspoc(jpi,jpj,jpk) ,    &
             &      prodgoc(jpi,jpj,jpk) , consgoc(jpi,jpj,jpk) ,  STAT=ierr(4) )
