Index: NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90
===================================================================
--- NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 10076)
+++ NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 10115)
@@ -190,7 +190,12 @@
             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) )    &
+               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) )
+               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
@@ -264,4 +269,5 @@
       !
       nitrfac (:,:,:) = 0._wp
+      nitrfac2(:,:,:) = 0._wp
       !
    END SUBROUTINE p4z_lim_init
Index: NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90
===================================================================
--- NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90	(revision 10076)
+++ NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90	(revision 10115)
@@ -131,5 +131,5 @@
 
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
-         IF( iom_use( "PH"     ) ) CALL iom_put( "PH"    , -1. * LOG10( hi(:,:,:) )          * tmask(:,:,:) )
+         IF( iom_use( "PH"     ) ) CALL iom_put( "PH"    , -1. * LOG10( MAX( rtrn, hi(:,:,:) ) )  * 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(:,:,:) )
Index: NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90
===================================================================
--- NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90	(revision 10076)
+++ NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90	(revision 10115)
@@ -111,5 +111,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 )
@@ -126,5 +126,6 @@
                zdenom    = zfoodlim / ( xkgraz2 + zfoodlim )
                zdenom2   = zdenom / ( zfood + rtrn )
-               zgraze2   = grazrat2 * zstep * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpmes) 
+               zgraze2   = grazrat2 * zstep * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpmes) * (1. - nitrfac(ji,jj,jk) )  
+ 
 
                zgrazd    = zgraze2  * xprefc   * zcompadi  * zdenom2 
@@ -142,9 +143,11 @@
 # if ! defined key_kriest
                zgrazffeg = grazflux  * zstep * 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)
 # endif
                zgrazffep = grazflux  * zstep *  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: NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90
===================================================================
--- NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90	(revision 10076)
+++ NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90	(revision 10115)
@@ -103,5 +103,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 )
@@ -115,5 +115,5 @@
                zdenom    = zfoodlim / ( xkgraz + zfoodlim )
                zdenom2   = zdenom / ( zfood + rtrn )
-               zgraze    = grazrat * zstep * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpzoo) 
+               zgraze    = grazrat * zstep * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpzoo) * (1. - nitrfac(ji,jj,jk)) 
 
                zgrazp    = zgraze  * xpref2p * zcompaph  * zdenom2 
Index: NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90
===================================================================
--- NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 10076)
+++ NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 10115)
@@ -36,5 +36,5 @@
 
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_par      ! structure of input par
-   INTEGER , PARAMETER :: nbtimes = 365  !: maximum number of times record in a file
+   INTEGER , PARAMETER :: nbtimes = 366  !: maximum number of times record in a file
    INTEGER  :: ntimes_par                ! number of time steps in a file
    REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:) :: par_varsw    !: PAR fraction of shortwave
Index: NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90
===================================================================
--- NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 10076)
+++ NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 10115)
@@ -67,5 +67,5 @@
       !
       INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zremip, zremik, zsiremin 
+      REAL(wp) ::   zremip, zremik, zsiremin, zammonic 
       REAL(wp) ::   zsatur, zsatur2, znusil, znusil2, zdep, zdepmin, zfactdep
       REAL(wp) ::   zbactfer, zorem, zorem2, zofer, zolimit
@@ -78,4 +78,5 @@
       REAL(wp), POINTER, DIMENSION(:,:  ) :: ztempbac
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepbac, zolimi, zdepprod, zw3d
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zoxyrem
       !!---------------------------------------------------------------------
       !
@@ -84,5 +85,5 @@
       ! Allocate temporary workspace
       CALL wrk_alloc( jpi, jpj,      ztempbac                  )
-      CALL wrk_alloc( jpi, jpj, jpk, zdepbac, zdepprod, zolimi )
+      CALL wrk_alloc( jpi, jpj, jpk, zdepbac, zdepprod, zolimi, zoxyrem )
 
       ! Initialisation of temprary arrys
@@ -129,9 +130,11 @@
                ! 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(ji,jj,jk) = 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(ji,jj,jk) = MAX( 0.e0, zoxyrem(ji,jj,jk) )
                !
             END DO
@@ -157,5 +160,6 @@
                tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) + zonitr - rdenita * denitnh4(ji,jj,jk)
                tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - o2nit * zonitr
-               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2 * rno3 * zonitr + rno3 * ( rdenita - 1. ) * denitnh4(ji,jj,jk)
+               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2 * rno3 * zonitr &
+                  &                + rno3 * ( rdenita - 1. ) * denitnh4(ji,jj,jk)
             END DO
          END DO
@@ -291,11 +295,12 @@
 
       DO jk = 1, jpkm1
-         tra(:,:,jk,jppo4) = tra(:,:,jk,jppo4) + zolimi (:,:,jk) + denitr(:,:,jk)
-         tra(:,:,jk,jpnh4) = tra(:,:,jk,jpnh4) + zolimi (:,:,jk) + denitr(:,:,jk)
+         tra(:,:,jk,jppo4) = tra(:,:,jk,jppo4) + zolimi (:,:,jk) + denitr(:,:,jk) + zoxyrem(:,:,jk)
+         tra(:,:,jk,jpnh4) = tra(:,:,jk,jpnh4) + zolimi (:,:,jk) + denitr(:,:,jk) + zoxyrem(:,:,jk)
          tra(:,:,jk,jpno3) = tra(:,:,jk,jpno3) - denitr (:,:,jk) * rdenit
-         tra(:,:,jk,jpdoc) = tra(:,:,jk,jpdoc) - zolimi (:,:,jk) - denitr(:,:,jk)
+         tra(:,:,jk,jpdoc) = tra(:,:,jk,jpdoc) - zolimi (:,:,jk) - denitr(:,:,jk) - zoxyrem(:,:,jk)
          tra(:,:,jk,jpoxy) = tra(:,:,jk,jpoxy) - zolimi (:,:,jk) * o2ut
-         tra(:,:,jk,jpdic) = tra(:,:,jk,jpdic) + zolimi (:,:,jk) + denitr(:,:,jk)
-         tra(:,:,jk,jptal) = tra(:,:,jk,jptal) + rno3 * ( zolimi(:,:,jk) + ( rdenit + 1.) * denitr(:,:,jk) )
+         tra(:,:,jk,jpdic) = tra(:,:,jk,jpdic) + zolimi (:,:,jk) + denitr(:,:,jk) + zoxyrem(:,:,jk)
+         tra(:,:,jk,jptal) = tra(:,:,jk,jptal) + rno3 * ( zolimi(:,:,jk) + zoxyrem(:,:,jk) &
+              &                                        + ( rdenit + 1.) * denitr(:,:,jk) )
       END DO
 
@@ -323,5 +328,5 @@
       !
       CALL wrk_dealloc( jpi, jpj,      ztempbac                  )
-      CALL wrk_dealloc( jpi, jpj, jpk, zdepbac, zdepprod, zolimi )
+      CALL wrk_dealloc( jpi, jpj, jpk, zdepbac, zdepprod, zolimi, zoxyrem )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_rem')
Index: NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90
===================================================================
--- NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 10076)
+++ NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 10115)
@@ -64,5 +64,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, zwstpoc
       REAL(wp) ::  ztrfer, ztrpo4, zwdust, zlight
@@ -323,6 +323,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 * fse3t(ji,jj,ikt) 
+            zsedsi (ji,jj) = (1.0 - zrivsil) * zsiloss * fse3t(ji,jj,ikt) 
 #endif
          END DO
@@ -354,19 +354,19 @@
 
 #if ! defined key_sed
-            ! 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.
-            zpdenit  = MIN( 0.5 * ( trb(ji,jj,ikt,jpno3) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 )
+            ! The 0.5 factor in zpdenit is to avoid negative NO3 concentration after
+            ! denitrification in the sediments. Not very clever, but simpliest option.
+            zpdenit  = MIN( 0.5 * ( trb(ji,jj,ikt,jpno3) - rtrn ) / rdenit, &
+                            zdenit2d(ji,jj) * zwstpoc * zrivno3 * (1. - nitrfac2(ji,jj,ikt) ) )
             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
-            sdenit(ji,jj) = rdenit * zpdenit / zdep
-            zsedc(ji,jj)   = (1. - zrivno3) * zwstpoc / zdep
+            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 * fse3t(ji,jj,ikt) 
+            zsedc(ji,jj)   = (1. - zrivno3) * zwstpoc * fse3t(ji,jj,ikt) 
 #endif
          END DO
@@ -415,21 +415,23 @@
       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 molC/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 * fse3t(:,:,jk) * tmask(:,:,jk)
+                 zwork1(:,:) = zwork1(:,:) + nitrpot(:,:,jk) * nitrfix * rno3 * zfact * fse3t(:,:,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 (:,:) * rno3 * zfact )
          ENDIF
       ELSE
-         IF( ln_diatrc )  &
-            &  trc2d(:,:,jp_pcs0_2d + 12) = nitrpot(:,:,1) * nitrfix * rno3 * 1.e+3 * rfact2r * fse3t(:,:,1) * tmask(:,:,1)
+         IF( ln_diatrc ) THEN 
+            zfact = 1.e+3 * rfact2r  !  conversion from molC/l/kt  to molC/m3/s
+            trc2d(:,:,jp_pcs0_2d + 12) = nitrpot(:,:,1) * nitrfix * rno3 * zfact * fse3t(:,:,1) * tmask(:,:,1)
+         ENDIF
       ENDIF
       !
Index: NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90
===================================================================
--- NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 10076)
+++ NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 10115)
@@ -450,11 +450,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)
@@ -529,15 +529,11 @@
       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
+         zrdenittot = glob_sum ( denitr(:,:,:) * rdenit * xnegtr(:,:,:) * cvol(:,:,:) ) ! denitrification in the water column
+         zsdenittot = glob_sum ( sdenit(:,:) * e1e2t(:,:) * tmask(:,:,1) )              ! denitrification in the sediments
+         CALL iom_put( "tdenit"  , ( zrdenittot + zsdenittot ) * xfact3 )               ! Total denitrification in molN/m3 
       ENDIF
 
