Index: branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 4737)
+++ branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 5837)
@@ -64,5 +64,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_prod( kt , jnt )
+   SUBROUTINE p4z_prod( kt , knt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_prod  ***
@@ -74,5 +74,5 @@
       !!---------------------------------------------------------------------
       !
-      INTEGER, INTENT(in) :: kt, jnt
+      INTEGER, INTENT(in) :: kt, knt
       !
       INTEGER  ::   ji, jj, jk
@@ -83,8 +83,8 @@
       REAL(wp) ::   zpislopen  , zpislope2n
       REAL(wp) ::   zrum, zcodel, zargu, zval
-      REAL(wp) ::   zrfact2
+      REAL(wp) ::   zfact
       CHARACTER (len=25) :: charout
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zmixnano, zmixdiat, zstrn
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpislopead, zpislopead2, zprdia, zprbio, zprdch, zprnch, zysopt   
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zmixnano, zmixdiat, zstrn, zw2d
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpislopead, zpislopead2, zprdia, zprbio, zprdch, zprnch, zysopt, zw3d   
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zprorca, zprorcad, zprofed, zprofen, zprochln, zprochld, zpronew, zpronewd
       !!---------------------------------------------------------------------
@@ -129,19 +129,17 @@
       END DO
 
-      IF( ln_newprod ) THEN
-         ! Impact of the day duration on phytoplankton growth
-         DO jk = 1, jpkm1
-            DO jj = 1 ,jpj
-               DO ji = 1, jpi
-                  IF( etot(ji,jj,jk) > 1.E-3 ) THEN
-                     zval = MAX( 1., zstrn(ji,jj) )
-                     zval = 1.5 * zval / ( 12. + zval )
-                     zprbio(ji,jj,jk) = prmax(ji,jj,jk) * zval
-                     zprdia(ji,jj,jk) = zprbio(ji,jj,jk)
-                  ENDIF
-               END DO
-            END DO
-         END DO
-      ENDIF
+      ! Impact of the day duration on phytoplankton growth
+      DO jk = 1, jpkm1
+         DO jj = 1 ,jpj
+            DO ji = 1, jpi
+               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
+                  zval = MAX( 1., zstrn(ji,jj) )
+                  zval = 1.5 * zval / ( 12. + zval )
+                  zprbio(ji,jj,jk) = prmax(ji,jj,jk) * zval
+                  zprdia(ji,jj,jk) = zprbio(ji,jj,jk)
+               ENDIF
+            END DO
+         END DO
+      END DO
 
       ! Maximum light intensity
@@ -157,17 +155,17 @@
                DO ji = 1, jpi
                   ! Computation of the P-I slope for nanos and diatoms
-                  IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                       ztn         = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )
                       zadap       = xadap * ztn / ( 2.+ ztn )
-                      zconctemp   = MAX( 0.e0 , trn(ji,jj,jk,jpdia) - xsizedia )
-                      zconctemp2  = trn(ji,jj,jk,jpdia) - zconctemp
+                      zconctemp   = MAX( 0.e0 , trb(ji,jj,jk,jpdia) - xsizedia )
+                      zconctemp2  = trb(ji,jj,jk,jpdia) - zconctemp
                       znanotot    = enano(ji,jj,jk) * zstrn(ji,jj)
                       zdiattot    = ediat(ji,jj,jk) * zstrn(ji,jj)
                       !
                       zpislopead (ji,jj,jk) = pislope * ( 1.+ zadap  * EXP( -znanotot ) )  &
-                         &                   * trn(ji,jj,jk,jpnch) /( trn(ji,jj,jk,jpphy) * 12. + rtrn)
+                         &                   * trb(ji,jj,jk,jpnch) /( trb(ji,jj,jk,jpphy) * 12. + rtrn)
                       !
-                      zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp) / ( trn(ji,jj,jk,jpdia) + rtrn )   &
-                         &                   * trn(ji,jj,jk,jpdch) /( trn(ji,jj,jk,jpdia) * 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
@@ -196,30 +194,32 @@
 
                   ! Computation of the P-I slope for nanos and diatoms
-                  IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                       ztn         = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )
                       zadap       = ztn / ( 2.+ ztn )
-                      zconctemp   = MAX( 0.e0 , trn(ji,jj,jk,jpdia) - xsizedia )
-                      zconctemp2  = trn(ji,jj,jk,jpdia) - zconctemp
+                      zconctemp   = MAX( 0.e0 , trb(ji,jj,jk,jpdia) - xsizedia )
+                      zconctemp2  = trb(ji,jj,jk,jpdia) - zconctemp
+                      znanotot    = enano(ji,jj,jk) * zstrn(ji,jj)
+                      zdiattot    = ediat(ji,jj,jk) * zstrn(ji,jj)
                       !
-                      zpislopead (ji,jj,jk) = pislope  * ( 1.+ zadap  * EXP( -0.21 * enano(ji,jj,jk) ) )
-                      zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp)  / ( trn(ji,jj,jk,jpdia) + rtrn )
-
-                      zpislopen =  zpislopead(ji,jj,jk) * trn(ji,jj,jk,jpnch)                &
-                        &          / ( trn(ji,jj,jk,jpphy) * 12.                  + rtrn )   &
+                      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) * trn(ji,jj,jk,jpdch)                &
-                        &          / ( trn(ji,jj,jk,jpdia) * 12.                  + 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 )
 
                       ! Computation of production function for Carbon
                       !  ---------------------------------------------
-                      zprbio(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen  * enano(ji,jj,jk) ) )
-                      zprdia(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) ) )
+                      zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen  * znanotot ) )
+                      zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpislope2n * zdiattot ) )
 
                       !  Computation of production function for Chlorophyll
                       !--------------------------------------------------
-                      zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen  * enano(ji,jj,jk) * zstrn(ji,jj) ) )
-                      zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) * zstrn(ji,jj) ) )
+                      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) ) )
                   ENDIF
                END DO
@@ -252,5 +252,5 @@
             DO ji = 1, jpi
 
-                IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+                IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                    !    Si/C of diatoms
                    !    ------------------------
@@ -258,8 +258,8 @@
                    !    Si/C is arbitrariliy increased for very high Si concentrations
                    !    to mimic the very high ratios observed in the Southern Ocean (silpot2)
-                  zlim  = trn(ji,jj,jk,jpsil) / ( trn(ji,jj,jk,jpsil) + xksi1 )
+                  zlim  = trb(ji,jj,jk,jpsil) / ( trb(ji,jj,jk,jpsil) + xksi1 )
                   zsilim = MIN( zprdia(ji,jj,jk) / ( prmax(ji,jj,jk) + rtrn ), xlimsi(ji,jj,jk) )
                   zsilfac = 4.4 * EXP( -4.23 * zsilim ) * MAX( 0.e0, MIN( 1., 2.2 * ( zlim - 0.5 ) )  ) + 1.e0
-                  zsiborn = trn(ji,jj,jk,jpsil) * trn(ji,jj,jk,jpsil) * trn(ji,jj,jk,jpsil)
+                  zsiborn = trb(ji,jj,jk,jpsil) * trb(ji,jj,jk,jpsil) * trb(ji,jj,jk,jpsil)
                   IF (gphit(ji,jj) < -30 ) THEN
                     zsilfac2 = 1. + 2. * zsiborn / ( zsiborn + xksi2**3 )
@@ -302,10 +302,10 @@
 !CDIR NOVERRCHK
             DO ji = 1, jpi
-               IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                   !  production terms for nanophyto.
-                  zprorca(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * trn(ji,jj,jk,jpphy) * rfact2
+                  zprorca(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * trb(ji,jj,jk,jpphy) * rfact2
                   zpronew(ji,jj,jk) = zprorca(ji,jj,jk) * xnanono3(ji,jj,jk) / ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) + rtrn )
                   !
-                  zratio = trn(ji,jj,jk,jpnfe) / ( trn(ji,jj,jk,jpphy) + rtrn )
+                  zratio = trb(ji,jj,jk,jpnfe) / ( trb(ji,jj,jk,jpphy) + rtrn )
                   zratio = zratio / fecnm 
                   zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) ) 
@@ -313,10 +313,10 @@
                   &             * ( 4. - 4.5 * xlimnfe(ji,jj,jk) / ( xlimnfe(ji,jj,jk) + 0.5 ) )    &
                   &             * biron(ji,jj,jk) / ( biron(ji,jj,jk) + concnfe(ji,jj,jk) )  &
-                  &             * zmax * trn(ji,jj,jk,jpphy) * rfact2
+                  &             * zmax * trb(ji,jj,jk,jpphy) * rfact2
                   !  production terms for diatomees
-                  zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * trn(ji,jj,jk,jpdia) * rfact2
+                  zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * trb(ji,jj,jk,jpdia) * rfact2
                   zpronewd(ji,jj,jk) = zprorcad(ji,jj,jk) * xdiatno3(ji,jj,jk) / ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) + rtrn )
                   !
-                  zratio = trn(ji,jj,jk,jpdfe) / ( trn(ji,jj,jk,jpdia) + rtrn )
+                  zratio = trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) + rtrn )
                   zratio = zratio / fecdm 
                   zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) ) 
@@ -324,5 +324,5 @@
                   &             * ( 4. - 4.5 * xlimdfe(ji,jj,jk) / ( xlimdfe(ji,jj,jk) + 0.5 ) )    &
                   &             * biron(ji,jj,jk) / ( biron(ji,jj,jk) + concdfe(ji,jj,jk) )  &
-                  &             * zmax * trn(ji,jj,jk,jpdia) * rfact2
+                  &             * zmax * trb(ji,jj,jk,jpdia) * rfact2
                ENDIF
             END DO
@@ -341,5 +341,5 @@
                      zprdch(ji,jj,jk) = zprdch(ji,jj,jk) * zmixdiat(ji,jj)
                   ENDIF
-                  IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                      !  production terms for nanophyto. ( chlorophyll )
                      znanotot = enano(ji,jj,jk) * zstrn(ji,jj)
@@ -365,17 +365,17 @@
 !CDIR NOVERRCHK
                DO ji = 1, jpi
-                  IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+                  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) * trn(ji,jj,jk,jpphy) * xlimphy(ji,jj,jk)
+                     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) * trn(ji,jj,jk,jpnch) * znanotot +rtrn )
+                     &                    / ( zpislopead(ji,jj,jk) * trb(ji,jj,jk,jpnch) * 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) * trn(ji,jj,jk,jpdia) * xlimdia(ji,jj,jk)
+                     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) * trn(ji,jj,jk,jpdch) * zdiattot +rtrn )
+                     &                    / ( zpislopead2(ji,jj,jk) * trb(ji,jj,jk,jpdch) * zdiattot +rtrn )
                   ENDIF
                END DO
@@ -412,51 +412,143 @@
      END DO
 
-     ! Total primary production per year
-     tpp = tpp + glob_sum( ( zprorca(:,:,:) + zprorcad(:,:,:) ) * cvol(:,:,:) )
-
-     IF( ln_diatrc ) THEN
-         !
-         zrfact2 = 1.e3 * rfact2r  ! conversion from mol/L/timestep into mol/m3/s
-         IF( lk_iomput ) THEN
-           IF( jnt == nrdttrc ) THEN
-              CALL iom_put( "PPPHY"  , zprorca (:,:,:) * zrfact2 * tmask(:,:,:) )  ! primary production by nanophyto
-              CALL iom_put( "PPPHY2" , zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) )  ! primary production by diatom
-              CALL iom_put( "PPNEWN" , zpronew (:,:,:) * zrfact2 * tmask(:,:,:) )  ! new primary production by nanophyto
-              CALL iom_put( "PPNEWD" , zpronewd(:,:,:) * zrfact2 * tmask(:,:,:) )  ! new primary production by diatom
-              CALL iom_put( "PBSi"   , zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) * zysopt(:,:,:) ) ! biogenic silica production
-              CALL iom_put( "PFeD"   , zprofed (:,:,:) * zrfact2 * tmask(:,:,:) )  ! biogenic iron production by diatom
-              CALL iom_put( "PFeN"   , zprofen (:,:,:) * zrfact2 * tmask(:,:,:) )  ! biogenic iron production by nanophyto
-              CALL iom_put( "Mumax"  , prmax(:,:,:) * tmask(:,:,:) )  ! Maximum growth rate
-              CALL iom_put( "MuN"    , zprbio(:,:,:) * xlimphy(:,:,:) * tmask(:,:,:) )  ! Realized growth rate for nanophyto
-              CALL iom_put( "MuD"    , zprdia(:,:,:) * xlimdia(:,:,:) * tmask(:,:,:) )  ! Realized growth rate for diatoms
-              CALL iom_put( "LNlight", zprbio (:,:,:) / (prmax(:,:,:) + rtrn) * tmask(:,:,:) )  ! light limitation term
-              CALL iom_put( "LDlight", zprdia (:,:,:) / (prmax(:,:,:) + rtrn) * tmask(:,:,:) )  ! light limitation term
-           ENDIF
-         ELSE
-              trc3d(:,:,:,jp_pcs0_3d + 4)  = zprorca (:,:,:) * zrfact2 * tmask(:,:,:)
-              trc3d(:,:,:,jp_pcs0_3d + 5)  = zprorcad(:,:,:) * zrfact2 * tmask(:,:,:)
-              trc3d(:,:,:,jp_pcs0_3d + 6)  = zpronew (:,:,:) * zrfact2 * tmask(:,:,:)
-              trc3d(:,:,:,jp_pcs0_3d + 7)  = zpronewd(:,:,:) * zrfact2 * tmask(:,:,:)
-              trc3d(:,:,:,jp_pcs0_3d + 8)  = zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) * zysopt(:,:,:)
-              trc3d(:,:,:,jp_pcs0_3d + 9)  = zprofed (:,:,:) * zrfact2 * tmask(:,:,:)
+
+    ! Total primary production per year
+    IF( iom_use( "tintpp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  &
+         & tpp = glob_sum( ( zprorca(:,:,:) + zprorcad(:,:,:) ) * cvol(:,:,:) )
+
+    IF( lk_iomput ) THEN
+       IF( knt == nrdttrc ) THEN
+          CALL wrk_alloc( jpi, jpj,      zw2d )
+          CALL wrk_alloc( jpi, jpj, jpk, zw3d )
+          zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
+          !
+          IF( iom_use( "PPPHY" ) .OR. iom_use( "PPPHY2" ) )  THEN
+              zw3d(:,:,:) = zprorca (:,:,:) * zfact * tmask(:,:,:)  ! primary production by nanophyto
+              CALL iom_put( "PPPHY"  , zw3d )
+              !
+              zw3d(:,:,:) = zprorcad(:,:,:) * zfact * tmask(:,:,:)  ! primary production by diatomes
+              CALL iom_put( "PPPHY2"  , zw3d )
+          ENDIF
+          IF( iom_use( "PPNEWN" ) .OR. iom_use( "PPNEWD" ) )  THEN
+              zw3d(:,:,:) = zpronew (:,:,:) * zfact * tmask(:,:,:)  ! new primary production by nanophyto
+              CALL iom_put( "PPNEWN"  , zw3d )
+              !
+              zw3d(:,:,:) = zpronewd(:,:,:) * zfact * tmask(:,:,:)  ! new primary production by diatomes
+              CALL iom_put( "PPNEWD"  , zw3d )
+          ENDIF
+          IF( iom_use( "PBSi" ) )  THEN
+              zw3d(:,:,:) = zprorcad(:,:,:) * zfact * tmask(:,:,:) * zysopt(:,:,:) ! biogenic silica production
+              CALL iom_put( "PBSi"  , zw3d )
+          ENDIF
+          IF( iom_use( "PFeN" ) .OR. iom_use( "PFeD" ) )  THEN
+              zw3d(:,:,:) = zprofen(:,:,:) * zfact * tmask(:,:,:)  ! biogenic iron production by nanophyto
+              CALL iom_put( "PFeN"  , zw3d )
+              !
+              zw3d(:,:,:) = zprofed(:,:,:) * zfact * tmask(:,:,:)  ! biogenic iron production by  diatomes
+              CALL iom_put( "PFeD"  , zw3d )
+          ENDIF
+          IF( iom_use( "Mumax" ) )  THEN
+              zw3d(:,:,:) = prmax(:,:,:) * tmask(:,:,:)   ! Maximum growth rate
+              CALL iom_put( "Mumax"  , zw3d )
+          ENDIF
+          IF( iom_use( "MuN" ) .OR. iom_use( "MuD" ) )  THEN
+              zw3d(:,:,:) = zprbio(:,:,:) * xlimphy(:,:,:) * tmask(:,:,:)  ! Realized growth rate for nanophyto
+              CALL iom_put( "MuN"  , zw3d )
+              !
+              zw3d(:,:,:) =  zprdia(:,:,:) * xlimdia(:,:,:) * tmask(:,:,:)  ! Realized growth rate for diatoms
+              CALL iom_put( "MuD"  , zw3d )
+          ENDIF
+          IF( iom_use( "LNlight" ) .OR. iom_use( "LDlight" ) )  THEN
+              zw3d(:,:,:) = zprbio (:,:,:) / (prmax(:,:,:) + rtrn) * tmask(:,:,:) ! light limitation term
+              CALL iom_put( "LNlight"  , zw3d )
+              !
+              zw3d(:,:,:) =  zprdia (:,:,:) / (prmax(:,:,:) + rtrn) * tmask(:,:,:)  ! light limitation term
+              CALL iom_put( "LDlight"  , zw3d )
+          ENDIF
+          IF( iom_use( "TPP" ) )  THEN
+              zw3d(:,:,:) = ( zprorca(:,:,:) + zprorcad(:,:,:) ) * zfact * tmask(:,:,:)  ! total primary production
+              CALL iom_put( "TPP"  , zw3d )
+          ENDIF
+          IF( iom_use( "TPNEW" ) )  THEN
+              zw3d(:,:,:) = ( zpronew(:,:,:) + zpronewd(:,:,:) ) * zfact * tmask(:,:,:)  ! total new production
+              CALL iom_put( "TPNEW"  , zw3d )
+          ENDIF
+          IF( iom_use( "TPBFE" ) )  THEN
+              zw3d(:,:,:) = ( zprofen(:,:,:) + zprofed(:,:,:) ) * zfact * tmask(:,:,:)  ! total biogenic iron production
+              CALL iom_put( "TPBFE"  , zw3d )
+          ENDIF
+          IF( iom_use( "INTPPPHY" ) .OR. iom_use( "INTPPPHY2" ) ) THEN  
+             zw2d(:,:) = 0.
+             DO jk = 1, jpkm1
+               zw2d(:,:) = zw2d(:,:) + zprorca (:,:,jk) * fse3t(:,:,jk) * zfact * tmask(:,:,jk)  ! vert. integrated  primary produc. by nano
+             ENDDO
+             CALL iom_put( "INTPPPHY" , zw2d )
+             !
+             zw2d(:,:) = 0.
+             DO jk = 1, jpkm1
+                zw2d(:,:) = zw2d(:,:) + zprorcad(:,:,jk) * fse3t(:,:,jk) * zfact * tmask(:,:,jk) ! vert. integrated  primary produc. by diatom
+             ENDDO
+             CALL iom_put( "INTPPPHY2" , zw2d )
+          ENDIF
+          IF( iom_use( "INTPP" ) ) THEN   
+             zw2d(:,:) = 0.
+             DO jk = 1, jpkm1
+                zw2d(:,:) = zw2d(:,:) + ( zprorca(:,:,jk) + zprorcad(:,:,jk) ) * fse3t(:,:,jk) * zfact * tmask(:,:,jk) ! vert. integrated pp
+             ENDDO
+             CALL iom_put( "INTPP" , zw2d )
+          ENDIF
+          IF( iom_use( "INTPNEW" ) ) THEN    
+             zw2d(:,:) = 0.
+             DO jk = 1, jpkm1
+                zw2d(:,:) = zw2d(:,:) + ( zpronew(:,:,jk) + zpronewd(:,:,jk) ) * fse3t(:,:,jk) * zfact * tmask(:,:,jk)  ! vert. integrated new prod
+             ENDDO
+             CALL iom_put( "INTPNEW" , zw2d )
+          ENDIF
+          IF( iom_use( "INTPBFE" ) ) THEN           !   total biogenic iron production  ( vertically integrated )
+             zw2d(:,:) = 0.
+             DO jk = 1, jpkm1
+                zw2d(:,:) = zw2d(:,:) + ( zprofen(:,:,jk) + zprofed(:,:,jk) ) * fse3t(:,:,jk) * zfact * tmask(:,:,jk) ! vert integr. bfe prod
+             ENDDO
+            CALL iom_put( "INTPBFE" , zw2d )
+          ENDIF
+          IF( iom_use( "INTPBSI" ) ) THEN           !   total biogenic silica production  ( vertically integrated )
+             zw2d(:,:) = 0.
+             DO jk = 1, jpkm1
+                zw2d(:,:) = zw2d(:,:) + zprorcad(:,:,jk) * zysopt(:,:,jk) * fse3t(:,:,jk) * zfact * tmask(:,:,jk)  ! vert integr. bsi prod
+             ENDDO
+             CALL iom_put( "INTPBSI" , zw2d )
+          ENDIF
+          IF( iom_use( "tintpp" ) )  CALL iom_put( "tintpp" , tpp * zfact )  !  global total integrated primary production molC/s
+          !
+          CALL wrk_dealloc( jpi, jpj,      zw2d )
+          CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
+       ENDIF
+     ELSE
+        IF( ln_diatrc ) THEN
+           zfact = 1.e+3 * rfact2r
+           trc3d(:,:,:,jp_pcs0_3d + 4)  = zprorca (:,:,:) * zfact * tmask(:,:,:)
+           trc3d(:,:,:,jp_pcs0_3d + 5)  = zprorcad(:,:,:) * zfact * tmask(:,:,:)
+           trc3d(:,:,:,jp_pcs0_3d + 6)  = zpronew (:,:,:) * zfact * tmask(:,:,:)
+           trc3d(:,:,:,jp_pcs0_3d + 7)  = zpronewd(:,:,:) * zfact * tmask(:,:,:)
+           trc3d(:,:,:,jp_pcs0_3d + 8)  = zprorcad(:,:,:) * zfact * tmask(:,:,:) * zysopt(:,:,:)
+           trc3d(:,:,:,jp_pcs0_3d + 9)  = zprofed (:,:,:) * zfact * tmask(:,:,:)
 #  if ! defined key_kriest
-              trc3d(:,:,:,jp_pcs0_3d + 10) = zprofen (:,:,:) * zrfact2 * tmask(:,:,:)
+           trc3d(:,:,:,jp_pcs0_3d + 10) = zprofen (:,:,:) * zfact * tmask(:,:,:)
 #  endif
-         ENDIF
-         !
-      ENDIF
-
-      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+        ENDIF
+     ENDIF
+
+     IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('prod')")
          CALL prt_ctl_trc_info(charout)
          CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi, jpj,      zmixnano, zmixdiat, zstrn                                                  )
-      CALL wrk_dealloc( jpi, jpj, jpk, zpislopead, zpislopead2, zprdia, zprbio, zprdch, zprnch, zysopt            ) 
-      CALL wrk_dealloc( jpi, jpj, jpk, zprorca, zprorcad, zprofed, zprofen, zprochln, zprochld, zpronew, zpronewd )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('p4z_prod')
-      !
+     ENDIF
+     !
+     CALL wrk_dealloc( jpi, jpj,      zmixnano, zmixdiat, zstrn                                                  )
+     CALL wrk_dealloc( jpi, jpj, jpk, zpislopead, zpislopead2, zprdia, zprbio, zprdch, zprnch, zysopt            ) 
+     CALL wrk_dealloc( jpi, jpj, jpk, zprorca, zprorcad, zprofed, zprofen, zprochln, zprochld, zpronew, zpronewd )
+     !
+     IF( nn_timing == 1 )  CALL timing_stop('p4z_prod')
+     !
    END SUBROUTINE p4z_prod
 
