Index: trunk/NEMOGCM/CONFIG/ORCA2_LIM_PISCES/EXP00/iodef.xml
===================================================================
--- trunk/NEMOGCM/CONFIG/ORCA2_LIM_PISCES/EXP00/iodef.xml	(revision 4995)
+++ trunk/NEMOGCM/CONFIG/ORCA2_LIM_PISCES/EXP00/iodef.xml	(revision 4996)
@@ -46,5 +46,17 @@
 	  <field field_ref="svoce"        name="vos"     long_name="sea_surface_y_velocity"    />
 	</file>
-	
+
+        <file id="file4" name_suffix="_bioscalar" description="pisces sms variables" >
+          <field field_ref="tdenit"   name="tdenit"   unit="TgN/yr" operation="instant" > tdenit * 14. * 86400. * 365. / 1e12 </field>
+          <field field_ref="tnfix"    name="tnfix"    unit="TgN/yr" operation="instant" > tnfix * 14. * 86400. * 365. / 1e12 </field>
+          <field field_ref="tcflx"    name="tcflx"    unit="PgC/yr" operation="instant" > tcflx * -1. * 12. * 86400. * 365. / 1e15 </field>
+          <field field_ref="tcflxcum" name="tcflxcum" unit="PgC"    operation="instant" > tcflxcum * -1. * 12. / 1e15 </field>
+          <field field_ref="tcexp"    name="tcexp"    unit="PgC/yr" operation="instant" > tcexp * 12. * 86400. * 365. / 1e15 </field>
+          <field field_ref="tintpp"   name="tintpp"   unit="PgC/yr" operation="instant" > tintpp * 12. * 86400. * 365. / 1e15 </field>
+          <field field_ref="pno3tot"  name="pno3tot"  unit="umolN"  > pno3tot * 16. / 122. * 1e6 </field>
+          <field field_ref="psiltot"  name="psiltot"  unit="umolC"  > psiltot * 1e6  </field>
+          <field field_ref="palktot"  name="palktot"  unit="umolC"  > palktot * 1e6  </field>
+        <field field_ref="pfertot"  name="pfertot"  unit="nmolFe" > pfertot * 1e9  </field>
+                                                                                                                             </file>
       </file_group>
 
@@ -53,5 +65,5 @@
       <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE.">  <!-- 5d files -->   
 
-	<file id="file4" name_suffix="_grid_T" description="ocean T grid variables" >
+	<file id="file5" name_suffix="_grid_T" description="ocean T grid variables" >
 	  <field field_ref="toce"         name="thetao"   long_name="sea_water_potential_temperature"               />
 	  <field field_ref="soce"         name="so"       long_name="sea_water_salinity"                            />
@@ -72,5 +84,5 @@
 	</file>
 	
-	<file id="file5" name_suffix="_grid_U" description="ocean U grid variables" >
+	<file id="file6" name_suffix="_grid_U" description="ocean U grid variables" >
 	  <field field_ref="uoce"         name="uo"      long_name="sea_water_x_velocity"      />
 	  <field field_ref="suoce"        name="uos"     long_name="sea_surface_x_velocity"    />
@@ -81,5 +93,5 @@
 	</file>
 	
-	<file id="file6" name_suffix="_grid_V" description="ocean V grid variables" >
+	<file id="file7" name_suffix="_grid_V" description="ocean V grid variables" >
 	  <field field_ref="voce"         name="vo"      long_name="sea_water_y_velocity"      />
 	  <field field_ref="svoce"        name="vos"     long_name="sea_surface_y_velocity"    />
@@ -90,10 +102,10 @@
 	</file>
 	
-	<file id="file7" name_suffix="_grid_W" description="ocean W grid variables" >
+	<file id="file8" name_suffix="_grid_W" description="ocean W grid variables" >
 	  <field field_ref="woce"         name="wo"      long_name="ocean vertical velocity"         />
 	  <field field_ref="avt"          name="difvho"  long_name="ocean_vertical_heat_diffusivity" />
 	</file>
 	
-	<file id="file8" name_suffix="_icemod" description="ice variables" >
+	<file id="file9" name_suffix="_icemod" description="ice variables" >
 	  <field field_ref="ice_pres"                     />
 	  <field field_ref="snowthic_cea" name="snd"     long_name="surface_snow_thickness"   />
@@ -114,5 +126,5 @@
       <file_group id="1m" output_freq="1mo" output_level="10" enabled=".TRUE."> <!-- real monthly files -->
 
-	<file id="file9" name_suffix="_ptrc_T" description="pisces sms variables" >
+	<file id="file10" name_suffix="_ptrc_T" description="pisces sms variables" >
           <field field_ref="DIC"      />
           <field field_ref="Alkalini" />
@@ -126,5 +138,5 @@
 	</file>
 	
-	<file id="file10" name_suffix="_diad_T" description="additional pisces diagnostics" >
+	<file id="file11" name_suffix="_diad_T" description="additional pisces diagnostics" >
           <field field_ref="Cflx"     />
           <field field_ref="Dpco2"    />
@@ -140,5 +152,5 @@
       <file_group id="1y"  output_freq="1y" output_level="10" enabled=".TRUE."> <!-- real yearly files -->
 
-	<file id="file11" name_suffix="_ptrc_T" description="pisces sms variables" >
+	<file id="file12" name_suffix="_ptrc_T" description="pisces sms variables" >
           <field field_ref="DIC"      />
           <field field_ref="Alkalini" />
@@ -167,5 +179,5 @@
 	</file>
 
-	<file id="file12" name_suffix="_diad_T" description="additional pisces diagnostics" >
+	<file id="file13" name_suffix="_diad_T" description="additional pisces diagnostics" >
           <field field_ref="PH"       />
           <field field_ref="CO3"      />
Index: trunk/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/iodef.xml
===================================================================
--- trunk/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/iodef.xml	(revision 4995)
+++ trunk/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/iodef.xml	(revision 4996)
@@ -30,5 +30,19 @@
       <file_group id="6h" output_freq="6h"  output_level="10" enabled=".TRUE."/> <!-- 6h files -->
      
-      <file_group id="1d" output_freq="1d"  output_level="10" enabled=".TRUE."/> <!-- 1d files -->
+      <file_group id="1d" output_freq="1d"  output_level="10" enabled=".TRUE."> <!-- 1d files -->
+        <file id="file1" name_suffix="_bioscalar" description="pisces sms variables" >
+           <field field_ref="tdenit"   name="tdenit"   unit="TgN/yr" operation="instant" > tdenit * 14. * 86400. * 365. / 1e12 </field>
+           <field field_ref="tnfix"    name="tnfix"    unit="TgN/yr" operation="instant" > tnfix * 14. * 86400. * 365. / 1e12 </field>
+           <field field_ref="tcflx"    name="tcflx"    unit="PgC/yr" operation="instant" > tcflx * -1. * 12. * 86400. * 365. / 1e15 </field>
+           <field field_ref="tcflxcum" name="tcflxcum" unit="PgC"    operation="instant" > tcflxcum * -1. * 12. / 1e15 </field>
+           <field field_ref="tcexp"    name="tcexp"    unit="PgC/yr" operation="instant" > tcexp * 12. * 86400. * 365. / 1e15 </field>
+           <field field_ref="tintpp"   name="tintpp"   unit="PgC/yr" operation="instant" > tintpp * 12. * 86400. * 365. / 1e15 </field>
+           <field field_ref="pno3tot"  name="pno3tot"  unit="umolN"  > pno3tot * 16. / 122. * 1e6 </field>
+           <field field_ref="psiltot"  name="psiltot"  unit="umolC"  > psiltot * 1e6  </field>
+           <field field_ref="palktot"  name="palktot"  unit="umolC"  > palktot * 1e6  </field>
+           <field field_ref="pfertot"  name="pfertot"  unit="nmolFe" > pfertot * 1e9  </field>
+        </file>
+      </file_group>
+
       <file_group id="3d" output_freq="3d"  output_level="10" enabled=".TRUE."/> <!-- 3d files -->    
       <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE."/>  <!-- 5d files -->   
@@ -36,5 +50,5 @@
       <file_group id="1m" output_freq="1mo" output_level="10" enabled=".TRUE."> <!-- real monthly files -->
 
-	<file id="file1" name_suffix="_ptrc_T" description="pisces sms variables" >
+	<file id="file2" name_suffix="_ptrc_T" description="pisces sms variables" >
           <field field_ref="DIC"      />
           <field field_ref="Alkalini" />
@@ -48,5 +62,5 @@
 	</file>
 	
-	<file id="file2" name_suffix="_diad_T" description="additional pisces diagnostics" >
+	<file id="file3" name_suffix="_diad_T" description="additional pisces diagnostics" >
           <field field_ref="Cflx"     />
           <field field_ref="Dpco2"    />
@@ -61,5 +75,5 @@
       <file_group id="1y"  output_freq="1y" output_level="10" enabled=".TRUE."> <!-- real yearly files -->
 
-	<file id="file3" name_suffix="_ptrc_T" description="pisces sms variables" >
+	<file id="file4" name_suffix="_ptrc_T" description="pisces sms variables" >
           <field field_ref="DIC"      />
           <field field_ref="Alkalini" />
@@ -88,5 +102,5 @@
 	</file>
 
-	<file id="file4" name_suffix="_diad_T" description="additional pisces diagnostics" >
+	<file id="file5" name_suffix="_diad_T" description="additional pisces diagnostics" >
           <field field_ref="PH"       />
           <field field_ref="CO3"      />
Index: trunk/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- trunk/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 4995)
+++ trunk/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 4996)
@@ -599,4 +599,8 @@
        <field id="EPSI100"     long_name="Export of Silicate at 100 m"             unit="mol/m2/s"                     />
        <field id="EPCAL100"    long_name="Export of Calcite at 100 m"              unit="mol/m2/s"                     />
+       <field id="EXPC"        long_name="Export of carbon"                        unit="molC/m2/s"  grid_ref="grid_T_3D" />
+       <field id="EXPFE"       long_name="Export of biogenic iron"                 unit="molFe/m2/s" grid_ref="grid_T_3D" />
+       <field id="EXPSI"       long_name="Export of Silicate"                      unit="molSi/m2/s" grid_ref="grid_T_3D" />
+       <field id="EXPCAL"      long_name="Export of Calcite"                       unit="molC/m2/s"  grid_ref="grid_T_3D" />
        <field id="Cflx"        long_name="DIC flux"                                unit="mol/m2/s"                     />
        <field id="Oflx"        long_name="Oxygen flux"                             unit="mol/m2/s"                     />
@@ -608,12 +612,27 @@
        <field id="Ironsed"     long_name="Iron deposition from sediment"           unit="mol/m2/s"  grid_ref="grid_T_3D"/>
 
+
        <!-- PISCES with Kriest parametisation : variables available with key_kriest -->
-       <field id="POCFlx"      long_name="Particulate organic C flux"              unit="mol/m2/s"   grid_ref="grid_T_3D" />
-       <field id="NumFlx"      long_name="Particle number flux"                    unit="nbr/m2/s"   grid_ref="grid_T_3D" />
-       <field id="SiFlx"       long_name="Biogenic Si flux"                        unit="mol/m2/s"   grid_ref="grid_T_3D" />
-       <field id="CaCO3Flx"    long_name="CaCO3 flux"                              unit="mol/m2/s"   grid_ref="grid_T_3D" />
-       <field id="xnum"        long_name="Number of particles in aggregats"        unit="-"          grid_ref="grid_T_3D" />
-       <field id="W1"          long_name="sinking speed of mass flux"              unit="m2/s"       grid_ref="grid_T_3D" />
-       <field id="W2"          long_name="sinking speed of number flux"            unit="m2/s"       grid_ref="grid_T_3D" />
+       <field id="EPN100"      long_name="Particulate number flux at 100 m"        unit="mol/m2/s"                     />
+       <field id="EXPN"        long_name="Particulate number flux"                 unit="mol/m2/s"   grid_ref="grid_T_3D" />
+       <field id="XNUM"        long_name="Number of particles in aggregats"        unit="-"          grid_ref="grid_T_3D" />
+       <field id="WSC"         long_name="sinking speed of mass flux"              unit="m2/s"       grid_ref="grid_T_3D" />
+       <field id="WSN"         long_name="sinking speed of number flux"            unit="m2/s"       grid_ref="grid_T_3D" />
+
+      <!-- dbio_T on T grid : variables available with key_diaar5 -->
+       <field id="TPP"         long_name="Total Primary production of phyto"                   unit="molC/m3/s"  grid_ref="grid_T_3D" />
+       <field id="TPNEW"       long_name="New Primary production of phyto"                     unit="molC/m3/s"  grid_ref="grid_T_3D" />
+       <field id="TPBFE"       long_name="Total biogenic iron production"                      unit="molFe/m3/s" grid_ref="grid_T_3D" />
+       <field id="INTDIC"      long_name="DIC content"                                         unit="kg/m2"   />
+       <field id="O2MIN"       long_name="Oxygen minimum concentration"                        unit="molC/L" />
+       <field id="ZO2MIN"      long_name="Depth of oxygen minimum concentration"               unit="m"       />
+       <field id="INTNFIX"     long_name="Nitrogen fixation rate : vert. integrated"           unit="molN/m2/s"  />
+       <field id="INTPPPHY"    long_name="Vertically integrated primary production by nanophy" unit="molC/m2/s"  />
+       <field id="INTPPPHY2"   long_name="Vertically integrated primary production by diatom"  unit="molC/m2/s"  />
+       <field id="INTPP"       long_name="Vertically integrated primary production by phyto"   unit="molC/m2/s"  />
+       <field id="INTPNEW"     long_name="Vertically integrated new primary production"        unit="molC/m2/s"  />
+       <field id="INTPBFE"     long_name="Vertically integrated of biogenic iron production"   unit="molFe/m2/s" />
+       <field id="INTPBSI"     long_name="Vertically integrated of biogenic Si production"     unit="molSi/m2/s" />
+       <field id="INTPCAL"     long_name="Vertically integrated of calcite production"         unit="molC/m2/s" />
 
        <!-- PISCES light : variables available with key_pisces_reduced -->
@@ -647,4 +666,17 @@
        <field id="qintC14b"    long_name="Cumulative air-sea flux of Bomb C14"      unit="mol/m2"     />
        <field id="fdecay"      long_name="Radiactive decay of Bomb C14"             unit="mol/m3"  grid_ref="grid_T_3D"  />
+     </field_group>
+
+    <field_group id="PISCES_scalar"  domain_ref="1point" >
+       <field id="pno3tot"         long_name="global mean nitrate concentration"                  unit="mol/l"   />
+       <field id="psiltot"         long_name="global mean silicate concentration"                 unit="mol/l"   />
+       <field id="palktot"         long_name="global mean alkalinity concentration"               unit="mol/l"   />
+       <field id="pfertot"         long_name="global mean iron concentration"                     unit="mol/l"   />
+       <field id="tcflx"           long_name="total Flux of Carbon out of the ocean"              unit="molC/s" />
+       <field id="tcflxcum"        long_name="cumulative total Flux of Carbon out of the ocean"   unit="molC/s" />
+       <field id="tcexp"           long_name="total Carbon export at 100m"                        unit="molC/s" />
+       <field id="tintpp"          long_name="global total integrated primary production"         unit="molC/s" />
+       <field id="tnfix"           long_name="global total nitrogen fixation"                     unit="molN/s" />
+       <field id="tdenit"          long_name="Total denitrification"                              unit="molN/s" />
      </field_group>
 
Index: trunk/NEMOGCM/NEMO/OFF_SRC/istate.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OFF_SRC/istate.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/OFF_SRC/istate.F90	(revision 4996)
@@ -12,5 +12,4 @@
    USE oce             ! ocean dynamics and active tracers 
    USE dom_oce         ! ocean space and time domain 
-   USE daymod          ! calendar
 
    IMPLICIT NONE
@@ -47,6 +46,4 @@
       rn2  (:,:,:) = 0.e0 
       !
-      CALL day_init                           ! model calendar (using both namelist and restart infos)
-
    END SUBROUTINE istate_init
 
Index: trunk/NEMOGCM/NEMO/TOP_SRC/C14b/trcsms_c14b.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/C14b/trcsms_c14b.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/C14b/trcsms_c14b.F90	(revision 4996)
@@ -290,14 +290,14 @@
       ENDIF
       !    
-      IF( ln_diatrc ) THEN
-         IF( lk_iomput ) THEN
-            CALL iom_put( "qtrC14b"  , qtr_c14  )
-            CALL iom_put( "qintC14b" , qint_c14 )
-            CALL iom_put( "fdecay"   , zdecay   )
-          ELSE
+      IF( lk_iomput ) THEN
+        CALL iom_put( "qtrC14b"  , qtr_c14  )
+        CALL iom_put( "qintC14b" , qint_c14 )
+        CALL iom_put( "fdecay"   , zdecay   )
+      ELSE
+         IF( ln_diatrc ) THEN
             trc2d(:,:  ,jp_c14b0_2d     ) = qtr_c14 (:,:)
             trc2d(:,:  ,jp_c14b0_2d + 1 ) = qint_c14(:,:)
             trc3d(:,:,:,jp_c14b0_3d     ) = zdecay  (:,:,:)
-          ENDIF
+         ENDIF
       ENDIF
 
Index: trunk/NEMOGCM/NEMO/TOP_SRC/C14b/trcwri_c14b.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/C14b/trcwri_c14b.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/C14b/trcwri_c14b.F90	(revision 4996)
@@ -37,9 +37,5 @@
       DO jn = jp_c14b0, jp_c14b1
          cltra = TRIM( ctrcnm(jn) )                  ! short title for tracer
-         IF( lk_vvl ) THEN
-            CALL iom_put( cltra, trn(:,:,:,jn) * fse3t_n(:,:,:) )
-         ELSE
-            CALL iom_put( cltra, trn(:,:,:,jn) )
-         ENDIF
+         CALL iom_put( cltra, trn(:,:,:,jn) )
       END DO
       !
Index: trunk/NEMOGCM/NEMO/TOP_SRC/CFC/trcsms_cfc.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/CFC/trcsms_cfc.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/CFC/trcsms_cfc.F90	(revision 4996)
@@ -185,18 +185,16 @@
             CALL iom_rstput( kt, nitrst, numrtw, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jn) )
          END DO
-      ENDIF
-      !                                              
-      IF( ln_diatrc ) THEN
-        !
-        IF( lk_iomput ) THEN
-           CALL iom_put( "qtrCFC11"  , qtr_cfc (:,:,1) )
-           CALL iom_put( "qintCFC11" , qint_cfc(:,:,1) )
-        ELSE
-           trc2d(:,:,jp_cfc0_2d    ) = qtr_cfc (:,:,1)
-           trc2d(:,:,jp_cfc0_2d + 1) = qint_cfc(:,:,1)
-        END IF
-        !
+      ENDIF                                            
+      !
+      IF( lk_iomput ) THEN
+         CALL iom_put( "qtrCFC11"  , qtr_cfc (:,:,1) )
+         CALL iom_put( "qintCFC11" , qint_cfc(:,:,1) )
+      ELSE
+         IF( ln_diatrc ) THEN
+            trc2d(:,:,jp_cfc0_2d    ) = qtr_cfc (:,:,1)
+            trc2d(:,:,jp_cfc0_2d + 1) = qint_cfc(:,:,1)
+         END IF
       END IF
- 
+      !
       IF( l_trdtrc ) THEN
           DO jn = jp_cfc0, jp_cfc1
Index: trunk/NEMOGCM/NEMO/TOP_SRC/CFC/trcwri_cfc.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/CFC/trcwri_cfc.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/CFC/trcwri_cfc.F90	(revision 4996)
@@ -37,9 +37,5 @@
       DO jn = jp_cfc0, jp_cfc1
          cltra = TRIM( ctrcnm(jn) )                  ! short title for tracer
-         IF( lk_vvl ) THEN
-            CALL iom_put( cltra, trn(:,:,:,jn) * fse3t_n(:,:,:) )
-         ELSE
-            CALL iom_put( cltra, trn(:,:,:,jn) )
-         ENDIF
+         CALL iom_put( cltra, trn(:,:,:,jn) )
       END DO
       !
Index: trunk/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcwri_my_trc.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcwri_my_trc.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcwri_my_trc.F90	(revision 4996)
@@ -37,9 +37,5 @@
       DO jn = jp_myt0, jp_myt1
          cltra = TRIM( ctrcnm(jn) )                  ! short title for tracer
-         IF( lk_vvl ) THEN
-            CALL iom_put( cltra, trn(:,:,:,jn) * fse3t_n(:,:,:) )
-         ELSE
-            CALL iom_put( cltra, trn(:,:,:,jn) )
-         ENDIF
+         CALL iom_put( cltra, trn(:,:,:,jn) )
       END DO
       !
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zbio.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zbio.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zbio.F90	(revision 4996)
@@ -110,5 +110,5 @@
       IF( nn_timing == 1 )  CALL timing_start('p2z_bio')
       !
-      IF( ln_diatrc ) THEN
+      IF( ln_diatrc .OR. lk_iomput ) THEN
          CALL wrk_alloc( jpi, jpj,     17, zw2d )
          CALL wrk_alloc( jpi, jpj, jpk, 3, zw3d )
@@ -122,5 +122,5 @@
 
       xksi(:,:) = 0.e0        ! zooplakton closure ( fbod)
-      IF( ln_diatrc ) THEN
+      IF( ln_diatrc .OR. lk_iomput ) THEN
          zw2d  (:,:,:) = 0.e0
          zw3d(:,:,:,:) = 0.e0
@@ -239,5 +239,5 @@
                   !  trend number 17 in p2zexp
                 ENDIF
-                IF( ln_diatrc ) THEN
+                IF( ln_diatrc .OR. lk_iomput ) THEN
                   ! convert fluxes in per day
                   ze3t = fse3t(ji,jj,jk) * 86400.
@@ -360,5 +360,5 @@
                   !  trend number 17 in p2zexp 
                 ENDIF
-                IF( ln_diatrc ) THEN
+                IF( ln_diatrc .OR. lk_iomput ) THEN
                   ! convert fluxes in per day
                   ze3t = fse3t(ji,jj,jk) * 86400.
@@ -390,6 +390,5 @@
       END DO
 
-      IF( ln_diatrc ) THEN
-         !
+      IF( ln_diatrc .OR. lk_iomput ) THEN
          DO jl = 1, 17 
             CALL lbc_lnk( zw2d(:,:,jl),'T', 1. )
@@ -398,28 +397,30 @@
             CALL lbc_lnk( zw3d(:,:,:,jl),'T', 1. )
          END DO
-         IF( lk_iomput ) THEN
+      ENDIF
+      IF( lk_iomput ) THEN
             ! Save diagnostics
-            CALL iom_put( "TNO3PHY", zw2d(:,:,1) )
-            CALL iom_put( "TNH4PHY", zw2d(:,:,2) )
-            CALL iom_put( "TPHYDOM", zw2d(:,:,3) )
-            CALL iom_put( "TPHYNH4", zw2d(:,:,4) )
-            CALL iom_put( "TPHYZOO", zw2d(:,:,5) )
-            CALL iom_put( "TPHYDET", zw2d(:,:,6) )
-            CALL iom_put( "TDETZOO", zw2d(:,:,7) )
-            CALL iom_put( "TZOODET", zw2d(:,:,8) )
-            CALL iom_put( "TZOOBOD", zw2d(:,:,9) )
-            CALL iom_put( "TZOONH4", zw2d(:,:,10) )
-            CALL iom_put( "TZOODOM", zw2d(:,:,11) )
-            CALL iom_put( "TNH4NO3", zw2d(:,:,12) )
-            CALL iom_put( "TDOMNH4", zw2d(:,:,13) )
-            CALL iom_put( "TDETNH4", zw2d(:,:,14) )
-            CALL iom_put( "TPHYTOT", zw2d(:,:,15) )
-            CALL iom_put( "TZOOTOT", zw2d(:,:,16) )
-            ! 
-            CALL iom_put( "FNO3PHY", zw3d(:,:,:,1) )
-            CALL iom_put( "FNH4PHY", zw3d(:,:,:,2) )
-            CALL iom_put( "FNH4NO3", zw3d(:,:,:,3) )
-            !
-         ELSE
+        CALL iom_put( "TNO3PHY", zw2d(:,:,1) )
+        CALL iom_put( "TNH4PHY", zw2d(:,:,2) )
+        CALL iom_put( "TPHYDOM", zw2d(:,:,3) )
+        CALL iom_put( "TPHYNH4", zw2d(:,:,4) )
+        CALL iom_put( "TPHYZOO", zw2d(:,:,5) )
+        CALL iom_put( "TPHYDET", zw2d(:,:,6) )
+        CALL iom_put( "TDETZOO", zw2d(:,:,7) )
+        CALL iom_put( "TZOODET", zw2d(:,:,8) )
+        CALL iom_put( "TZOOBOD", zw2d(:,:,9) )
+        CALL iom_put( "TZOONH4", zw2d(:,:,10) )
+        CALL iom_put( "TZOODOM", zw2d(:,:,11) )
+        CALL iom_put( "TNH4NO3", zw2d(:,:,12) )
+        CALL iom_put( "TDOMNH4", zw2d(:,:,13) )
+        CALL iom_put( "TDETNH4", zw2d(:,:,14) )
+        CALL iom_put( "TPHYTOT", zw2d(:,:,15) )
+        CALL iom_put( "TZOOTOT", zw2d(:,:,16) )
+         ! 
+        CALL iom_put( "FNO3PHY", zw3d(:,:,:,1) )
+        CALL iom_put( "FNH4PHY", zw3d(:,:,:,2) )
+        CALL iom_put( "FNH4NO3", zw3d(:,:,:,3) )
+         !
+       ELSE
+          IF( ln_diatrc ) THEN
             !
             trc2d(:,:,jp_pcs0_2d    ) = zw2d(:,:,1) 
@@ -467,5 +468,5 @@
       ENDIF
       !
-      IF( ln_diatrc ) THEN
+      IF( ln_diatrc .OR. lk_iomput ) THEN
          CALL wrk_dealloc( jpi, jpj,     17, zw2d )
          CALL wrk_dealloc( jpi, jpj, jpk, 3, zw3d )
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zexp.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zexp.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zexp.F90	(revision 4996)
@@ -127,8 +127,8 @@
  
       ! Oa & Ek: diagnostics depending on jpdia2d !          left as example
-      IF( ln_diatrc ) THEN
-         IF( lk_iomput ) THEN   ;   CALL iom_put( "SEDPOC" , sedpocn )
-         ELSE                   ;   trc2d(:,:,jp_pcs0_2d + 18) = sedpocn(:,:)
-         ENDIF
+      IF( lk_iomput ) THEN  
+         CALL iom_put( "SEDPOC" , sedpocn )
+      ELSE
+         IF( ln_diatrc )           trc2d(:,:,jp_pcs0_2d + 18) = sedpocn(:,:)
       ENDIF
 
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zsed.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zsed.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zsed.F90	(revision 4996)
@@ -108,16 +108,24 @@
       END DO
 
-      IF( ln_diatrc ) THEN 
-         CALL wrk_alloc( jpi, jpj, zw2d )
-         zw2d(:,:) =  ztra(:,:,1) * fse3t(:,:,1) * 86400.
-         DO jk = 2, jpkm1
-            zw2d(:,:) = zw2d(:,:) + ztra(:,:,jk) * fse3t(:,:,jk) * 86400.
-         END DO
-         IF( lk_iomput )  THEN
-           CALL iom_put( "TDETSED", zw2d )
-         ELSE
-           trc2d(:,:,jp_pcs0_2d + 7) = zw2d(:,:)
+      IF( lk_iomput )  THEN
+         IF( iom_use( "TDETSED" ) ) THEN
+            CALL wrk_alloc( jpi, jpj, zw2d )
+            zw2d(:,:) =  ztra(:,:,1) * fse3t(:,:,1) * 86400.
+            DO jk = 2, jpkm1
+               zw2d(:,:) = zw2d(:,:) + ztra(:,:,jk) * fse3t(:,:,jk) * 86400.
+            END DO
+            CALL iom_put( "TDETSED", zw2d )
+            CALL wrk_dealloc( jpi, jpj, zw2d )
          ENDIF
-         CALL wrk_dealloc( jpi, jpj, zw2d )
+      ELSE
+         IF( ln_diatrc ) THEN 
+            CALL wrk_alloc( jpi, jpj, zw2d )
+            zw2d(:,:) =  ztra(:,:,1) * fse3t(:,:,1) * 86400.
+            DO jk = 2, jpkm1
+               zw2d(:,:) = zw2d(:,:) + ztra(:,:,jk) * fse3t(:,:,jk) * 86400.
+            END DO
+            trc2d(:,:,jp_pcs0_2d + 7) = zw2d(:,:)
+            CALL wrk_dealloc( jpi, jpj, zw2d )
+         ENDIF
       ENDIF
       !
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 4996)
@@ -299,17 +299,15 @@
       !  Output of some diagnostics variables
       !     ---------------------------------
-      IF( ln_diatrc .AND. lk_iomput ) THEN
-         IF( jnt == nrdttrc ) THEN
-            CALL iom_put("Fe3"    , zFe3   (:,:,:)       * tmask(:,:,:) )   ! Fe3+
-            CALL iom_put("FeL1"   , zFeL1  (:,:,:)       * tmask(:,:,:) )   ! FeL1
-            CALL iom_put("TL1"    , zTL1   (:,:,:)       * tmask(:,:,:) )   ! TL1
-            CALL iom_put("Totlig" , ztotlig(:,:,:)       * tmask(:,:,:) )   ! TL
-            CALL iom_put("Biron"  , biron  (:,:,:) * 1e9 * tmask(:,:,:) )   ! biron
-            IF( ln_fechem ) THEN
-               CALL iom_put("Fe2" , zFe2   (:,:,:)       * tmask(:,:,:) )   ! Fe2+
-               CALL iom_put("FeL2", zFeL2  (:,:,:)       * tmask(:,:,:) )   ! FeL2
-               CALL iom_put("FeP" , zFeP   (:,:,:)       * tmask(:,:,:) )   ! FeP
-               CALL iom_put("TL2" , zTL2   (:,:,:)       * tmask(:,:,:) )   ! TL2
-            ENDIF
+      IF( lk_iomput .AND. jnt == nrdttrc ) THEN
+         IF( iom_use("Fe3")    )  CALL iom_put("Fe3"    , zFe3   (:,:,:)       * tmask(:,:,:) )   ! Fe3+
+         IF( iom_use("FeL1")   )  CALL iom_put("FeL1"   , zFeL1  (:,:,:)       * tmask(:,:,:) )   ! FeL1
+         IF( iom_use("TL1")    )  CALL iom_put("TL1"    , zTL1   (:,:,:)       * tmask(:,:,:) )   ! TL1
+         IF( iom_use("Totlig") )  CALL iom_put("Totlig" , ztotlig(:,:,:)       * tmask(:,:,:) )   ! TL
+         IF( iom_use("Biron")  )  CALL iom_put("Biron"  , biron  (:,:,:) * 1e9 * tmask(:,:,:) )   ! biron
+         IF( ln_fechem ) THEN
+            IF( iom_use("Fe2")  ) CALL iom_put("Fe2"    , zFe2   (:,:,:)       * tmask(:,:,:) )   ! Fe2+
+            IF( iom_use("FeL2") ) CALL iom_put("FeL2"   , zFeL2  (:,:,:)       * tmask(:,:,:) )   ! FeL2
+            IF( iom_use("FeP")  ) CALL iom_put("FeP"    , zFeP   (:,:,:)       * tmask(:,:,:) )   ! FeP
+            IF( iom_use("TL2")  ) CALL iom_put("TL2"    , zTL2   (:,:,:)       * tmask(:,:,:) )   ! TL2
          ENDIF
       ENDIF
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90	(revision 4996)
@@ -89,5 +89,5 @@
       REAL(wp) ::   zyr_dec, zdco2dt
       CHARACTER (len=25) :: charout
-      REAL(wp), POINTER, DIMENSION(:,:) :: zkgco2, zkgo2, zh2co3, zoflx 
+      REAL(wp), POINTER, DIMENSION(:,:) :: zkgco2, zkgo2, zh2co3, zoflx, zw2d 
       !!---------------------------------------------------------------------
       !
@@ -196,7 +196,9 @@
       END DO
 
-      t_oce_co2_flx = t_oce_co2_flx + glob_sum( oce_co2(:,:) )      ! Cumulative Total Flux of Carbon
-      t_atm_co2_flx = glob_sum( satmco2(:,:) * e1e2t(:,:) )         ! Total atmospheric pCO2
-
+      t_oce_co2_flx     = glob_sum( oce_co2(:,:) )                    !  Total Flux of Carbon
+      t_oce_co2_flx_cum = t_oce_co2_flx_cum + t_oce_co2_flx       !  Cumulative Total Flux of Carbon
+!      t_atm_co2_flx     = glob_sum( satmco2(:,:) * e1e2t(:,:) )       ! Total atmospheric pCO2
+      t_atm_co2_flx     =  atcco2      ! Total atmospheric pCO2
+ 
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('flx ')")
@@ -205,12 +207,32 @@
       ENDIF
 
-      IF( ln_diatrc ) THEN
-         IF( lk_iomput ) THEN
-            CALL iom_put( "Cflx" , oce_co2(:,:) / e1e2t(:,:) / rfact ) 
-            CALL iom_put( "Oflx" , zoflx(:,:) * 1000 * tmask(:,:,1)  )
-            CALL iom_put( "Kg"   , zkgco2(:,:) * tmask(:,:,1) )
-            CALL iom_put( "Dpco2", ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) )
-            CALL iom_put( "Dpo2" , ( atcox * patm(:,:) - trn(:,:,1,jpoxy) / ( chemc(:,:,2) + rtrn ) ) * tmask(:,:,1) )
-         ELSE
+      IF( lk_iomput ) THEN
+         CALL wrk_alloc( jpi, jpj, zw2d )  
+         IF( iom_use( "Cflx"  ) )  THEN
+            zw2d(:,:) = oce_co2(:,:) / e1e2t(:,:) / rfact
+            CALL iom_put( "Cflx"     , zw2d ) 
+         ENDIF
+         IF( iom_use( "Oflx"  ) )  THEN
+            zw2d(:,:) =  zoflx(:,:) * 1000 * tmask(:,:,1)
+            CALL iom_put( "Oflx" , zw2d )
+         ENDIF
+         IF( iom_use( "Kg"    ) )  THEN
+            zw2d(:,:) =  zkgco2(:,:) * tmask(:,:,1)
+            CALL iom_put( "Kg"   , zw2d )
+         ENDIF
+         IF( iom_use( "Dpco2" ) ) THEN
+           zw2d(:,:) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1)
+           CALL iom_put( "Dpco2" ,  zw2d )
+         ENDIF
+         IF( iom_use( "Dpo2" ) )  THEN
+           zw2d(:,:) = ( atcox * patm(:,:) - trn(:,:,1,jpoxy) / ( chemc(:,:,2) + rtrn ) ) * tmask(:,:,1)
+           CALL iom_put( "Dpo2"  , zw2d )
+         ENDIF
+         IF( iom_use( "tcflx" ) )  CALL iom_put( "tcflx"    , t_oce_co2_flx * rfact2r )   ! molC/s
+         CALL iom_put( "tcflxcum" , t_oce_co2_flx_cum )      ! molC
+         !
+         CALL wrk_dealloc( jpi, jpj, zw2d )
+      ELSE
+         IF( ln_diatrc ) THEN
             trc2d(:,:,jp_pcs0_2d    ) = oce_co2(:,:) / e1e2t(:,:) / rfact 
             trc2d(:,:,jp_pcs0_2d + 1) = zoflx(:,:) * 1000 * tmask(:,:,1) 
@@ -290,6 +312,6 @@
       !
       oce_co2(:,:)  = 0._wp                ! Initialization of Flux of Carbon
+      t_oce_co2_flx = 0._wp
       t_atm_co2_flx = 0._wp
-      t_oce_co2_flx = 0._wp
       !
       CALL p4z_patm( nit000 )
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 4996)
@@ -187,14 +187,12 @@
       END DO
       !
-      IF( ln_diatrc .AND. lk_iomput .AND. jnt == nrdttrc ) THEN        ! save output diagnostics
-        !
-        CALL iom_put( "xfracal", xfracal(:,:,:) * tmask(:,:,:) )  ! euphotic layer deptht
-        CALL iom_put( "LNnut"  , xlimphy(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
-        CALL iom_put( "LDnut"  , xlimdia(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
-        CALL iom_put( "LNFe"   , xlimnfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
-        CALL iom_put( "LDFe"   , xlimdfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
-        !
+      !
+      IF( lk_iomput .AND. jnt == nrdttrc ) THEN        ! save output diagnostics
+        IF( iom_use( "xfracal" ) ) CALL iom_put( "xfracal", xfracal(:,:,:) * tmask(:,:,:) )  ! euphotic layer deptht
+        IF( iom_use( "LNnut"   ) ) CALL iom_put( "LNnut"  , xlimphy(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
+        IF( iom_use( "LDnut"   ) ) CALL iom_put( "LDnut"  , xlimdia(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
+        IF( iom_use( "LNFe"    ) ) CALL iom_put( "LNFe"   , xlimnfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+        IF( iom_use( "LDFe"    ) ) CALL iom_put( "LDFe"   , xlimdfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
       ENDIF
-
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_lim')
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90	(revision 4996)
@@ -146,18 +146,14 @@
       END DO
       !
-      IF( ln_diatrc )  THEN
-         !
-         IF( lk_iomput ) THEN
-            zrfact2 = 1.e3 * rfact2r
-            CALL iom_put( "PH"    , -1. * LOG10( hi(:,:,:) )                * tmask(:,:,:) )
-            CALL iom_put( "CO3"   ,        zco3    (:,:,:) * 1e+3           * tmask(:,:,:) )
-            CALL iom_put( "CO3sat",        aksp    (:,:,:) * 1e+3 / calcon  * tmask(:,:,:) )
-            CALL iom_put( "DCAL"  ,        zcaldiss(:,:,:) * zrfact2        * tmask(:,:,:) )
-         ELSE
-            trc3d(:,:,:,jp_pcs0_3d    ) = -1. * LOG10( hi(:,:,:) ) * tmask(:,:,:)
-            trc3d(:,:,:,jp_pcs0_3d + 1) = zco3(:,:,:)              * tmask(:,:,:)
-            trc3d(:,:,:,jp_pcs0_3d + 2) = aksp(:,:,:) / calcon     * tmask(:,:,:)
-         ENDIF
-         !
+
+      IF( lk_iomput ) 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(:,:,:) )
+      ELSE
+         trc3d(:,:,:,jp_pcs0_3d    ) = -1. * LOG10( hi(:,:,:) ) * tmask(:,:,:)
+         trc3d(:,:,:,jp_pcs0_3d + 1) = zco3(:,:,:)              * tmask(:,:,:)
+         trc3d(:,:,:,jp_pcs0_3d + 2) = aksp(:,:,:) / calcon     * tmask(:,:,:)
       ENDIF
       !
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90	(revision 4996)
@@ -83,6 +83,5 @@
       REAL(wp) :: zgrazfffp, zgrazfffg, zgrazffep, zgrazffeg
       CHARACTER (len=25) :: charout
-      REAL(wp) :: zrfact2
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zgrazing
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zgrazing, zw3d
 
       !!---------------------------------------------------------------------
@@ -90,5 +89,5 @@
       IF( nn_timing == 1 )  CALL timing_start('p4z_meso')
       !
-      IF( ln_diatrc .AND. lk_iomput ) THEN
+      IF( lk_iomput ) THEN
          CALL wrk_alloc( jpi, jpj, jpk, zgrazing )
          zgrazing(:,:,:) = 0._wp
@@ -186,5 +185,5 @@
 
               ! Total grazing ( grazing by microzoo is already computed in p4zmicro )
-              IF( ln_diatrc .AND. lk_iomput )  zgrazing(ji,jj,jk) = zgraztot
+              IF( lk_iomput )  zgrazing(ji,jj,jk) = zgraztot
 
               !    Mesozooplankton efficiency
@@ -249,8 +248,15 @@
       END DO
       !
-      IF( ln_diatrc .AND. lk_iomput .AND. jnt == nrdttrc ) THEN
-         zrfact2 = 1.e3 * rfact2r
-         CALL iom_put( "GRAZ2", zgrazing(:,:,:) * zrfact2 * tmask(:,:,:) )  ! Total grazing of phyto by zooplankton
-         CALL iom_put( "PCAL" , prodcal(:,:,:)  * zrfact2 * tmask(:,:,:) )  ! Calcite production
+      IF( lk_iomput .AND. jnt == nrdttrc ) THEN
+         CALL wrk_alloc( jpi, jpj, jpk, zw3d )
+         IF( iom_use( "GRAZ2" ) ) THEN
+            zw3d(:,:,:) = zgrazing(:,:,:) * 1.e+3 * rfact2r * tmask(:,:,:)  !   Total grazing of phyto by zooplankton
+            CALL iom_put( "GRAZ2", zw3d )
+         ENDIF
+         IF( iom_use( "PCAL" ) ) THEN
+            zw3d(:,:,:) = prodcal(:,:,:) * 1.e+3 * rfact2r * tmask(:,:,:)   !  Calcite production
+            CALL iom_put( "PCAL", zw3d )  
+         ENDIF
+         CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
       ENDIF
       !
@@ -261,5 +267,5 @@
       ENDIF
       !
-      IF( ln_diatrc .AND. lk_iomput )  CALL wrk_dealloc( jpi, jpj, jpk, zgrazing )
+      IF( lk_iomput )  CALL wrk_dealloc( jpi, jpj, jpk, zgrazing )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_meso')
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90	(revision 4996)
@@ -79,6 +79,5 @@
       REAL(wp) :: zgrazp, zgrazm, zgrazsd
       REAL(wp) :: zgrazmf, zgrazsf, zgrazpf
-      REAL(wp) :: zrfact2
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zgrazing
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zgrazing, zw3d
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
@@ -86,5 +85,5 @@
       IF( nn_timing == 1 )  CALL timing_start('p4z_micro')
       !
-      IF( ln_diatrc .AND. lk_iomput )  CALL wrk_alloc( jpi, jpj, jpk, zgrazing )
+      IF( lk_iomput )  CALL wrk_alloc( jpi, jpj, jpk, zgrazing )
       !
       DO jk = 1, jpkm1
@@ -191,7 +190,11 @@
       END DO
       !
-      IF( ln_diatrc .AND. lk_iomput .AND. jnt == nrdttrc ) THEN
-         zrfact2 = 1.e3 * rfact2r
-         CALL iom_put( "GRAZ1" , zgrazing(:,:,:) * zrfact2 * tmask(:,:,:) )  ! Total grazing of phyto by zooplankton
+      IF( lk_iomput .AND. jnt == nrdttrc ) THEN
+         CALL wrk_alloc( jpi, jpj, jpk, zw3d )
+         IF( iom_use( "GRAZ1" ) ) THEN
+            zw3d(:,:,:) = zgrazing(:,:,:) * 1.e+3 * rfact2r * tmask(:,:,:)  !  Total grazing of phyto by zooplankton
+            CALL iom_put( "GRAZ1", zw3d )
+         ENDIF
+         CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
       ENDIF
       !
@@ -202,5 +205,5 @@
       ENDIF
       !
-      IF( ln_diatrc .AND. lk_iomput )  CALL wrk_dealloc( jpi, jpj, jpk, zgrazing )
+      IF( lk_iomput )  CALL wrk_dealloc( jpi, jpj, jpk, zgrazing )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_micro')
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 4996)
@@ -252,16 +252,14 @@
       END DO
 
-      IF( ln_diatrc ) THEN        ! save output diagnostics
-        !
-        IF( lk_iomput ) THEN
-           IF( jnt == nrdttrc ) THEN
-              CALL iom_put( "Heup", heup(:,:  ) * tmask(:,:,1) )  ! euphotic layer deptht
-              CALL iom_put( "PAR" , emoy(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
-           ENDIF
-        ELSE
-           trc2d(:,:,  jp_pcs0_2d + 10) = heup(:,:  ) * tmask(:,:,1)  
-           trc3d(:,:,:,jp_pcs0_3d + 3)  = etot(:,:,:) * tmask(:,:,:)
+      IF( lk_iomput ) THEN
+        IF( jnt == nrdttrc  ) THEN
+           IF( iom_use( "Heup" ) ) CALL iom_put( "Heup", heup(:,:  ) * tmask(:,:,1) )  ! euphotic layer deptht
+           IF( iom_use( "PAR"  ) ) CALL iom_put( "PAR" , emoy(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
         ENDIF
-        !
+      ELSE
+         IF( ln_diatrc ) THEN        ! save output diagnostics
+            trc2d(:,:,  jp_pcs0_2d + 10) = heup(:,:  ) * tmask(:,:,1)  
+            trc3d(:,:,:,jp_pcs0_3d + 3)  = etot(:,:,:) * tmask(:,:,:)
+         ENDIF
       ENDIF
       !
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 4996)
@@ -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
       !!---------------------------------------------------------------------
@@ -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. jnt == nrdttrc )  )  &
+         & tpp = glob_sum( ( zprorca(:,:,:) + zprorcad(:,:,:) ) * cvol(:,:,:) )
+
+    IF( lk_iomput ) THEN
+       IF( jnt == 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
 
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 4996)
@@ -78,8 +78,8 @@
       REAL(wp) ::   zofer2
 #endif
-      REAL(wp) ::   zonitr, zstep, zrfact2
+      REAL(wp) ::   zonitr, zstep, zfact
       CHARACTER (len=25) :: charout
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztempbac 
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepbac, zolimi, zdepprod
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztempbac
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepbac, zolimi, zdepprod, zw3d
       !!---------------------------------------------------------------------
       !
@@ -315,9 +315,19 @@
       END DO
 
-      IF( ln_diatrc .AND. lk_iomput .AND. jnt == nrdttrc ) THEN
-          zrfact2 = 1.e3 * rfact2r
-          CALL iom_put( "REMIN" , zolimi(:,:,:) * tmask(:,:,:) * zrfact2 )  ! Remineralisation rate
-          CALL iom_put( "DENIT" , denitr(:,:,:) * rdenit * rno3 * tmask(:,:,:) * zrfact2  )  ! Denitrification
-      ENDIF
+      IF( jnt == nrdttrc ) THEN
+          CALL wrk_alloc( jpi, jpj, jpk, zw3d )
+          zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
+          !
+          IF( iom_use( "REMIN" ) )  THEN
+              zw3d(:,:,:) = zolimi(:,:,:) * tmask(:,:,:) * zfact !  Remineralisation rate
+              CALL iom_put( "REMIN"  , zw3d )
+          ENDIF
+          IF( iom_use( "DENIT" ) )  THEN
+              zw3d(:,:,:) = denitr(:,:,:) * rdenit * rno3 * tmask(:,:,:) * zfact ! Denitrification
+              CALL iom_put( "DENIT"  , zw3d )
+          ENDIF
+          !
+          CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
+       ENDIF
 
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 4996)
@@ -80,5 +80,4 @@
    REAL(wp), PUBLIC :: rivdininput, rivdipinput, rivdsiinput
 
-   REAL(wp) :: ryyss                    !: number of seconds per year 
 
    !!* Substitution
@@ -207,6 +206,4 @@
       !
       IF( nn_timing == 1 )  CALL timing_start('p4z_sbc_init')
-      !
-      ryyss   = nyear_len(1) * rday    ! number of seconds per year and per month
       !
       !                            !* set file information
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 4996)
@@ -33,6 +33,4 @@
 
    !! * Module variables
-   REAL(wp) :: ryyss                    !: number of seconds per year 
-   REAL(wp) :: r1_ryyss                 !: inverse of ryyss
    REAL(wp) :: r1_rday                  !: inverse of rday
 
@@ -83,7 +81,5 @@
       !
       IF( kt == nittrc000 .AND. jnt == 1 )  THEN
-         ryyss    = nyear_len(1) * rday    ! number of seconds per year and per month
          r1_rday  = 1. / rday
-         r1_ryyss = 1. / ryyss
          IF( ln_check_mass .AND. lwp)  &
            &  CALL ctl_opn( numnit, 'nitrogen.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
@@ -119,7 +115,8 @@
          !
          trn(:,:,1,jpfer) = trn(:,:,1,jpfer) + zironice(:,:) 
-         !                                              
-         IF( ln_diatrc .AND. lk_iomput .AND. jnt == nrdttrc )   &
+         ! 
+         IF( lk_iomput .AND. jnt == nrdttrc .AND. iom_use( "Ironice" ) )   &
             &   CALL iom_put( "Ironice", zironice(:,:) * 1.e+3 * rfact2r * fse3t(:,:,1) * tmask(:,:,1) ) ! iron flux from ice
+         !
          CALL wrk_dealloc( jpi, jpj, zironice )
          !                                              
@@ -150,15 +147,15 @@
          trn(:,:,1,jpsil) = trn(:,:,1,jpsil) + zsidep  (:,:)
          trn(:,:,:,jpfer) = trn(:,:,:,jpfer) + zirondep(:,:,:) 
-         !                                              
-         IF( ln_diatrc ) THEN
-            zfact = 1.e+3 * rfact2r
-            IF( lk_iomput ) THEN
-               IF( jnt == nrdttrc ) THEN
-                  CALL iom_put( "Irondep", zirondep(:,:,1) * zfact * fse3t(:,:,1) * tmask(:,:,1) ) ! surface downward dust depo of iron
-                  CALL iom_put( "pdust"  , dust(:,:) / ( wdust * rday )  * tmask(:,:,1) ) ! dust concentration at surface
-               ENDIF
-            ELSE
-               trc2d(:,:,jp_pcs0_2d + 11) = zirondep(:,:,1) * zfact * fse3t(:,:,1) * tmask(:,:,1)
+         ! 
+         IF( lk_iomput ) THEN
+            IF( jnt == nrdttrc ) THEN
+                IF( iom_use( "Irondep" ) )   &
+                &  CALL iom_put( "Irondep", zirondep(:,:,1) * 1.e+3 * rfact2r * fse3t(:,:,1) * tmask(:,:,1) ) ! surface downward dust depo of iron
+                IF( iom_use( "pdust" ) )   &
+                &  CALL iom_put( "pdust"  , dust(:,:) / ( wdust * rday )  * tmask(:,:,1) ) ! dust concentration at surface
             ENDIF
+         ELSE                                    
+            IF( ln_diatrc )  &
+              &  trc2d(:,:,jp_pcs0_2d + 11) = zirondep(:,:,1) * 1.e+3 * rfact2r * fse3t(:,:,1) * tmask(:,:,1)
          ENDIF
          CALL wrk_dealloc( jpi, jpj,      zpdep, zsidep )
@@ -190,5 +187,5 @@
          trn(:,:,:,jpfer) = trn(:,:,:,jpfer) + ironsed(:,:,:) * rfact2
          !
-         IF( ln_diatrc .AND. lk_iomput .AND. jnt == nrdttrc )   &
+         IF( lk_iomput .AND. jnt == nrdttrc .AND. iom_use( "Ironsed" ) )   &
             &   CALL iom_put( "Ironsed", ironsed(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! iron inputs from sediments
       ENDIF
@@ -199,8 +196,7 @@
          trn(:,:,:,jpfer) = trn(:,:,:,jpfer) + hydrofe(:,:,:) * rfact2
          !
-         IF( ln_diatrc .AND. lk_iomput .AND. jnt == nrdttrc )   &
+         IF( lk_iomput .AND. jnt == nrdttrc .AND. iom_use( "HYDR" ) )   &
             &   CALL iom_put( "HYDR", hydrofe(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! hydrothermal iron input
       ENDIF
-
 
       ! OA: Warning, the following part is necessary, especially with Kriest
@@ -385,28 +381,35 @@
       END DO
 
- 
-      IF( ln_check_mass ) THEN
-        ! Global budget of N SMS : denitrification in the water column and in the sediment
-         !                          nitrogen fixation by the diazotrophs
-         ! --------------------------------------------------------------------------------
-         zrdenittot   = glob_sum ( denitr(:,:,:) * rdenit * xnegtr(:,:,:) * cvol(:,:,:) )
-         zsdenittot   = glob_sum ( zwork4(:,:)   * e1e2t(:,:) )
-         znitrpottot  = glob_sum ( znitrpot(:,:,:) * nitrfix              * cvol(:,:,:) )
-         IF( kt == nitend .AND. jnt == nrdttrc ) THEN
-            zfact = 1.e+3 * rfact2r * rno3 * ryyss * 14. / 1e12
-            IF(lwp) WRITE(numnit,9100) ndastp, znitrpottot * nitrfix * zfact, zrdenittot * zfact , zsdenittot * zfact
+      ! Global budget of N SMS : denitrification in the water column and in the sediment
+      !                          nitrogen fixation by the diazotrophs
+      ! --------------------------------------------------------------------------------
+      zrdenittot   = glob_sum ( denitr(:,:,:) * rdenit * xnegtr(:,:,:) * cvol(:,:,:) )
+      zsdenittot   = glob_sum ( zwork4(:,:)   * e1e2t(:,:) )
+      znitrpottot  = glob_sum ( znitrpot(:,:,:) * nitrfix              * cvol(:,:,:) )
+      zfact = 1.e+3 * rfact2r * rno3 * 14. / 1.e12 * ryyss   ! conversion molC/l/kt ----> TgN/m3/year
+      !
+      IF( ln_check_mass .AND. ( kt == nitend .AND. jnt == nrdttrc ) .AND. ( lwp )  )  &
+         &  WRITE(numnit,9100) ndastp, znitrpottot * zfact  , &
+         &                             zrdenittot  * zfact  , &
+         &                             zsdenittot  * zfact
+      !
+      IF( lk_iomput ) THEN
+         IF( jnt == 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"  , znitrpot(:,:,:) * nitrfix * zfact * tmask(:,:,:) )  ! nitrogen fixation 
+            IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", zwork4(:,:) * zfact * tmask(:,:,1) )  ! Nitrate reduction in the sediments
+            IF( iom_use("tnfix"  ) ) CALL iom_put( "tnfix"  , znitrpottot * zfact  )               ! Global  nitrogen fixation
+            IF( iom_use("tdenit" ) ) CALL iom_put( "tdenit" , zrdenittot  * zfact  )               ! Total denitrification
+            IF( iom_use("INTNFIX") ) THEN   ! nitrogen fixation rate in ocean ( vertically integrated )
+               zwork1(:,:) = 0.
+               DO jk = 1, jpkm1
+                 zwork1(:,:) = zwork1(:,:) + znitrpot(:,:,jk) * nitrfix * zfact * fse3t(:,:,jk) * tmask(:,:,jk)
+               ENDDO
+               CALL iom_put( "INTNFIX" , zwork1 ) 
+            ENDIF
          ENDIF
-       ENDIF
-      !
-      IF( ln_diatrc ) THEN
-         zfact = 1.e+3 * rfact2r
-         IF( lk_iomput ) THEN
-            IF( jnt == nrdttrc ) THEN
-               CALL iom_put( "Nfix"  , znitrpot(:,:,:) * nitrfix * rno3 * zfact * tmask(:,:,:) )  ! nitrogen fixation 
-               CALL iom_put( "Sdenit", zwork4(:,:)               * rno3 * zfact * tmask(:,:,1) )  ! Nitrate reduction in the sediments
-            ENDIF
-         ELSE
-            trc2d(:,:,jp_pcs0_2d + 12) = znitrpot(:,:,1) * nitrfix * zfact * fse3t(:,:,1) * tmask(:,:,1)
-         ENDIF
+      ELSE
+         IF( ln_diatrc )  &
+            &  trc2d(:,:,jp_pcs0_2d + 12) = znitrpot(:,:,1) * nitrfix * rno3 * 1.e+3 * rfact2r * fse3t(:,:,1) * tmask(:,:,1)
       ENDIF
       !
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90	(revision 4996)
@@ -41,5 +41,5 @@
 #endif
 
-   INTEGER  :: iksed  = 10
+   INTEGER  :: ik100
 
 #if  defined key_kriest
@@ -94,7 +94,7 @@
       REAL(wp) ::   zagg , zaggfe, zaggdoc, zaggdoc2, zaggdoc3
       REAL(wp) ::   zfact, zwsmax, zmax, zstep
-      REAL(wp) ::   zrfact2
-      INTEGER  ::   ik1
       CHARACTER (len=25) :: charout
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zw2d
       !!---------------------------------------------------------------------
       !
@@ -235,24 +235,61 @@
       END DO
 
-     ! Total primary production per year
-     t_oce_co2_exp = t_oce_co2_exp + glob_sum( ( sinking(:,:,iksed+1) + sinking2(:,:,iksed+1) ) * e1e2t(:,:) * tmask(:,:,1) )
+
+     ! Total carbon export per year
+     IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. jnt == nrdttrc )  )  &
+        &   t_oce_co2_exp = glob_sum( ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * e1e2t(:,:) * tmask(:,:,1) )
      !
-     IF( ln_diatrc ) THEN
-         zrfact2 = 1.e3 * rfact2r
-         ik1  = iksed + 1
-         IF( lk_iomput ) THEN
-           IF( jnt == nrdttrc ) THEN
-              CALL iom_put( "EPC100"  , ( sinking(:,:,ik1) + sinking2(:,:,ik1) ) * zrfact2 * tmask(:,:,1) ) ! Export of carbon at 100m
-              CALL iom_put( "EPFE100" , ( sinkfer(:,:,ik1) + sinkfer2(:,:,ik1) ) * zrfact2 * tmask(:,:,1) ) ! Export of iron at 100m
-              CALL iom_put( "EPCAL100",   sinkcal(:,:,ik1)                       * zrfact2 * tmask(:,:,1) ) ! Export of calcite  at 100m
-              CALL iom_put( "EPSI100" ,   sinksil(:,:,ik1)                       * zrfact2 * tmask(:,:,1) ) ! Export of biogenic silica at 100m
-           ENDIF
-         ELSE
-           trc2d(:,:,jp_pcs0_2d + 4) = sinking (:,:,ik1) * zrfact2 * tmask(:,:,1)
-           trc2d(:,:,jp_pcs0_2d + 5) = sinking2(:,:,ik1) * zrfact2 * tmask(:,:,1)
-           trc2d(:,:,jp_pcs0_2d + 6) = sinkfer (:,:,ik1) * zrfact2 * tmask(:,:,1)
-           trc2d(:,:,jp_pcs0_2d + 7) = sinkfer2(:,:,ik1) * zrfact2 * tmask(:,:,1)
-           trc2d(:,:,jp_pcs0_2d + 8) = sinksil (:,:,ik1) * zrfact2 * tmask(:,:,1)
-           trc2d(:,:,jp_pcs0_2d + 9) = sinkcal (:,:,ik1) * zrfact2 * tmask(:,:,1)
+     IF( lk_iomput ) THEN
+       IF( jnt == 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( "EPC100" ) )  THEN
+              zw2d(:,:) = ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * zfact * tmask(:,:,1) ! Export of carbon at 100m
+              CALL iom_put( "EPC100"  , zw2d )
+          ENDIF
+          IF( iom_use( "EPFE100" ) )  THEN
+              zw2d(:,:) = ( sinkfer(:,:,ik100) + sinkfer2(:,:,ik100) ) * zfact * tmask(:,:,1) ! Export of iron at 100m
+              CALL iom_put( "EPFE100"  , zw2d )
+          ENDIF
+          IF( iom_use( "EPCAL100" ) )  THEN
+              zw2d(:,:) = sinkcal(:,:,ik100) * zfact * tmask(:,:,1) ! Export of calcite at 100m
+              CALL iom_put( "EPCAL100"  , zw2d )
+          ENDIF
+          IF( iom_use( "EPSI100" ) )  THEN
+              zw2d(:,:) =  sinksil(:,:,ik100) * zfact * tmask(:,:,1) ! Export of bigenic silica at 100m
+              CALL iom_put( "EPSI100"  , zw2d )
+          ENDIF
+          IF( iom_use( "EXPC" ) )  THEN
+              zw3d(:,:,:) = ( sinking(:,:,:) + sinking2(:,:,:) ) * zfact * tmask(:,:,:) ! Export of carbon in the water column
+              CALL iom_put( "EXPC"  , zw3d )
+          ENDIF
+          IF( iom_use( "EXPFE" ) )  THEN
+              zw3d(:,:,:) = ( sinkfer(:,:,:) + sinkfer2(:,:,:) ) * zfact * tmask(:,:,:) ! Export of iron 
+              CALL iom_put( "EXPFE"  , zw3d )
+          ENDIF
+          IF( iom_use( "EXPCAL" ) )  THEN
+              zw3d(:,:,:) = sinkcal(:,:,:) * zfact * tmask(:,:,:) ! Export of calcite 
+              CALL iom_put( "EXPCAL"  , zw3d )
+          ENDIF
+          IF( iom_use( "EXPSI" ) )  THEN
+              zw3d(:,:,:) = sinksil(:,:,:) * zfact * tmask(:,:,:) ! Export of bigenic silica
+              CALL iom_put( "EXPSI"  , zw3d )
+          ENDIF
+          IF( iom_use( "tcexp" ) )  CALL iom_put( "tcexp" , t_oce_co2_exp * zfact )   ! molC/s
+          ! 
+          CALL wrk_dealloc( jpi, jpj,      zw2d )
+          CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
+        ENDIF
+      ELSE
+         IF( ln_diatrc ) THEN
+            zfact = 1.e3 * rfact2r
+            trc2d(:,:,jp_pcs0_2d + 4) = sinking (:,:,ik100) * zfact * tmask(:,:,1)
+            trc2d(:,:,jp_pcs0_2d + 5) = sinking2(:,:,ik100) * zfact * tmask(:,:,1)
+            trc2d(:,:,jp_pcs0_2d + 6) = sinkfer (:,:,ik100) * zfact * tmask(:,:,1)
+            trc2d(:,:,jp_pcs0_2d + 7) = sinkfer2(:,:,ik100) * zfact * tmask(:,:,1)
+            trc2d(:,:,jp_pcs0_2d + 8) = sinksil (:,:,ik100) * zfact * tmask(:,:,1)
+            trc2d(:,:,jp_pcs0_2d + 9) = sinkcal (:,:,ik100) * zfact * tmask(:,:,1)
          ENDIF
       ENDIF
@@ -272,5 +309,14 @@
       !!                  ***  ROUTINE p4z_sink_init  ***
       !!----------------------------------------------------------------------
-
+      INTEGER :: jk
+
+      ik100 = 10        !  last level where depth less than 100 m
+      DO jk = jpkm1, 1, -1
+         IF( gdept_1d(jk) > 100. )  ik100 = jk - 1
+      END DO
+      IF (lwp) WRITE(numout,*)
+      IF (lwp) WRITE(numout,*) ' Level corresponding to 100m depth ',  ik100 + 1
+      IF (lwp) WRITE(numout,*)
+      !
       t_oce_co2_exp = 0._wp
       !
@@ -300,8 +346,10 @@
       REAL(wp) :: zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5
       REAL(wp) :: zval1, zval2, zval3, zval4
-      REAL(wp) :: zrfact2
+      REAL(wp) :: zfact
       INTEGER  :: ik1
       CHARACTER (len=25) :: charout
       REAL(wp), POINTER, DIMENSION(:,:,:) :: znum3d 
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zw2d
       !!---------------------------------------------------------------------
       !
@@ -477,35 +525,77 @@
 
      ! Total primary production per year
-     t_oce_co2_exp = t_oce_co2_exp + glob_sum( ( sinking(:,:,:) ) * cvol(:,:,:) )
+     t_oce_co2_exp = t_oce_co2_exp + glob_sum( ( sinking(:,:,ik100) * e1e2t(:,:) * tmask(:,:,1) )
      !
-      IF( ln_diatrc ) THEN
-         !
-         ik1 = iksed + 1
-         zrfact2 = 1.e3 * rfact2r
-         IF( jnt == nrdttrc ) THEN
-           CALL iom_put( "POCFlx"  , sinking (:,:,:)      * zrfact2 * tmask(:,:,:) )  ! POC export
-           CALL iom_put( "NumFlx"  , sinking2 (:,:,:)     * zrfact2 * tmask(:,:,:) )  ! Num export
-           CALL iom_put( "SiFlx"   , sinksil (:,:,:)      * zrfact2 * tmask(:,:,:) )  ! Silica export
-           CALL iom_put( "CaCO3Flx", sinkcal (:,:,:)      * zrfact2 * tmask(:,:,:) )  ! Calcite export
-           CALL iom_put( "xnum"    , znum3d  (:,:,:)                * tmask(:,:,:) )  ! Number of particles in aggregats
-           CALL iom_put( "W1"      , wsbio3  (:,:,:)                * tmask(:,:,:) )  ! sinking speed of POC
-           CALL iom_put( "W2"      , wsbio4  (:,:,:)                * tmask(:,:,:) )  ! sinking speed of aggregats
+     IF( lk_iomput ) THEN
+        IF( jnt == 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( "EPC100" ) )  THEN
+              zw2d(:,:) = sinking(:,:,ik100) * zfact * tmask(:,:,1) ! Export of carbon at 100m
+              CALL iom_put( "EPC100"  , zw2d )
+          ENDIF
+          IF( iom_use( "EPN100" ) )  THEN
+              zw2d(:,:) = sinking2(:,:,ik100) * zfact * tmask(:,:,1) ! Export of number of aggregates ?
+              CALL iom_put( "EPN100"  , zw2d )
+          ENDIF
+          IF( iom_use( "EPCAL100" ) )  THEN
+              zw2d(:,:) = sinkcal(:,:,ik100) * zfact * tmask(:,:,1) ! Export of calcite at 100m
+              CALL iom_put( "EPCAL100"  , zw2d )
+          ENDIF
+          IF( iom_use( "EPSI100" ) )  THEN
+              zw2d(:,:) = sinksil(:,:,ik100) * zfact * tmask(:,:,1) ! Export of bigenic silica at 100m
+              CALL iom_put( "EPSI100"  , zw2d )
+          ENDIF
+          IF( iom_use( "EXPC" ) )  THEN
+              zw3d(:,:,:) = sinking(:,:,:) * zfact * tmask(:,:,:) ! Export of carbon in the water column
+              CALL iom_put( "EXPC"  , zw3d )
+          ENDIF
+          IF( iom_use( "EXPN" ) )  THEN
+              zw3d(:,:,:) = sinking(:,:,:) * zfact * tmask(:,:,:) ! Export of carbon in the water column
+              CALL iom_put( "EXPN"  , zw3d )
+          ENDIF
+          IF( iom_use( "EXPCAL" ) )  THEN
+              zw3d(:,:,:) = sinkcal(:,:,:) * zfact * tmask(:,:,:) ! Export of calcite 
+              CALL iom_put( "EXPCAL"  , zw3d )
+          ENDIF
+          IF( iom_use( "EXPSI" ) )  THEN
+              zw3d(:,:,:) = sinksil(:,:,:) * zfact * tmask(:,:,:) ! Export of bigenic silica
+              CALL iom_put( "EXPSI"  , zw3d )
+          ENDIF
+          IF( iom_use( "XNUM" ) )  THEN
+              zw3d(:,:,:) =  znum3d(:,:,:) * tmask(:,:,:) !  Number of particles on aggregats
+              CALL iom_put( "XNUM"  , zw3d )
+          ENDIF
+          IF( iom_use( "WSC" ) )  THEN
+              zw3d(:,:,:) = wsbio3(:,:,:) * tmask(:,:,:) ! Sinking speed of carbon particles
+              CALL iom_put( "WSC"  , zw3d )
+          ENDIF
+          IF( iom_use( "WSN" ) )  THEN
+              zw3d(:,:,:) = wsbio4(:,:,:) * tmask(:,:,:) ! Sinking speed of particles number
+              CALL iom_put( "WSN"  , zw3d )
+          ENDIF
+          !
+          CALL wrk_dealloc( jpi, jpj,      zw2d )
+          CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
+      ELSE
+         IF( ln_diatrc ) THEN
+            zfact = 1.e3 * rfact2r
+            trc2d(:,:  ,jp_pcs0_2d + 4)  = sinking (:,:,ik100)  * zfact * tmask(:,:,1)
+            trc2d(:,:  ,jp_pcs0_2d + 5)  = sinking2(:,:,ik100)  * zfact * tmask(:,:,1)
+            trc2d(:,:  ,jp_pcs0_2d + 6)  = sinkfer (:,:,ik100)  * zfact * tmask(:,:,1)
+            trc2d(:,:  ,jp_pcs0_2d + 7)  = sinksil (:,:,ik100)  * zfact * tmask(:,:,1)
+            trc2d(:,:  ,jp_pcs0_2d + 8)  = sinkcal (:,:,ik100)  * zfact * tmask(:,:,1)
+            trc3d(:,:,:,jp_pcs0_3d + 11) = sinking (:,:,:)      * zfact * tmask(:,:,:)
+            trc3d(:,:,:,jp_pcs0_3d + 12) = sinking2(:,:,:)      * zfact * tmask(:,:,:)
+            trc3d(:,:,:,jp_pcs0_3d + 13) = sinksil (:,:,:)      * zfact * tmask(:,:,:)
+            trc3d(:,:,:,jp_pcs0_3d + 14) = sinkcal (:,:,:)      * zfact * tmask(:,:,:)
+            trc3d(:,:,:,jp_pcs0_3d + 15) = znum3d  (:,:,:)              * tmask(:,:,:)
+            trc3d(:,:,:,jp_pcs0_3d + 16) = wsbio3  (:,:,:)              * tmask(:,:,:)
+            trc3d(:,:,:,jp_pcs0_3d + 17) = wsbio4  (:,:,:)              * tmask(:,:,:)
          ENDIF
-# if ! defined key_iomput
-         trc2d(:,:  ,jp_pcs0_2d + 4)  = sinking (:,:,ik1)    * zrfact2 * tmask(:,:,1)
-         trc2d(:,:  ,jp_pcs0_2d + 5)  = sinking2(:,:,ik1)    * zrfact2 * tmask(:,:,1)
-         trc2d(:,:  ,jp_pcs0_2d + 6)  = sinkfer (:,:,ik1)    * zrfact2 * tmask(:,:,1)
-         trc2d(:,:  ,jp_pcs0_2d + 7)  = sinksil (:,:,ik1)    * zrfact2 * tmask(:,:,1)
-         trc2d(:,:  ,jp_pcs0_2d + 8)  = sinkcal (:,:,ik1)    * zrfact2 * tmask(:,:,1)
-         trc3d(:,:,:,jp_pcs0_3d + 11) = sinking (:,:,:)      * zrfact2 * tmask(:,:,:)
-         trc3d(:,:,:,jp_pcs0_3d + 12) = sinking2(:,:,:)      * zrfact2 * tmask(:,:,:)
-         trc3d(:,:,:,jp_pcs0_3d + 13) = sinksil (:,:,:)      * zrfact2 * tmask(:,:,:)
-         trc3d(:,:,:,jp_pcs0_3d + 14) = sinkcal (:,:,:)      * zrfact2 * tmask(:,:,:)
-         trc3d(:,:,:,jp_pcs0_3d + 15) = znum3d  (:,:,:)                * tmask(:,:,:)
-         trc3d(:,:,:,jp_pcs0_3d + 16) = wsbio3  (:,:,:)                * tmask(:,:,:)
-         trc3d(:,:,:,jp_pcs0_3d + 17) = wsbio4  (:,:,:)                * tmask(:,:,:)
-# endif
-        !
       ENDIF
+
       !
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
@@ -663,4 +753,12 @@
       END DO
       !
+      ik100 = 10        !  last level where depth less than 100 m
+      DO jk = jpkm1, 1, -1
+         IF( gdept_1d(jk) > 100. )  iksed = jk - 1
+      END DO
+      IF (lwp) WRITE(numout,*)
+      IF (lwp) WRITE(numout,*) ' Level corresponding to 100m depth ',  ik100 + 1
+      IF (lwp) WRITE(numout,*)
+      !
       t_oce_co2_exp = 0._wp
       !
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 4996)
@@ -120,5 +120,4 @@
          ENDDO
       ENDIF
-
       CALL p4z_lys( kt )             ! Compute CaCO3 saturation
       CALL p4z_flx( kt )             ! Compute surface fluxes
@@ -151,5 +150,5 @@
       END IF
       !
-      CALL p4z_chk_mass( kt ) ! Mass conservation checking
+      IF( lk_iomput .OR. ln_check_mass )  CALL p4z_chk_mass( kt ) ! Mass conservation checking
 
       IF ( lwm .AND. kt == nittrc000 ) CALL FLUSH    ( numonp )     ! flush output namelist PISCES
@@ -328,4 +327,10 @@
          ENDIF
          !
+         IF( iom_varid( numrtr, 'tcflxcum', ldstop = .FALSE. ) > 0 ) THEN  ! cumulative total flux of carbon
+            CALL iom_get( numrtr, 'tcflxcum' , t_oce_co2_flx_cum  )
+         ELSE
+            t_oce_co2_flx_cum = 0._wp
+         ENDIF
+         !
       ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
          IF( kt == nitrst ) THEN
@@ -337,4 +342,5 @@
          CALL iom_rstput( kt, nitrst, numrtw, 'Silicalim', xksi(:,:) )
          CALL iom_rstput( kt, nitrst, numrtw, 'Silicamax', xksimax(:,:) )
+         CALL iom_rstput( kt, nitrst, numrtw, 'tcflxcum', t_oce_co2_flx_cum )
       ENDIF
       !
@@ -398,5 +404,6 @@
       !!---------------------------------------------------------------------
       !
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index      
+      INTEGER , INTENT( in ) ::   kt      ! ocean time-step index      
+      REAL(wp)               ::  zfact       
       !!
       !!---------------------------------------------------------------------
@@ -409,5 +416,6 @@
       ENDIF
 
-      IF( ln_check_mass .AND. kt == nitend ) THEN      !   Compute the budget of NO3, ALK, Si, Fer
+      IF( iom_use( "pno3tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
+         !   Compute the budget of NO3, ALK, Si, Fer
          no3budget = glob_sum( (   trn(:,:,:,jpno3) + trn(:,:,:,jpnh4)  &
             &                    + trn(:,:,:,jpphy) + trn(:,:,:,jpdia)  &
@@ -417,13 +425,28 @@
             &                    + trn(:,:,:,jpgoc)                     &
 #endif
-            &                    + trn(:,:,:,jpdoc)                     ) * cvol(:,:,:)  ) 
-         ! 
+            &                    + trn(:,:,:,jpdoc)                     ) * cvol(:,:,:)  )
+         !
+         no3budget = no3budget / areatot
+         CALL iom_put( "pno3tot", no3budget )
+      ENDIF
+      !
+      IF( iom_use( "psiltot" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
          silbudget = glob_sum( (   trn(:,:,:,jpsil) + trn(:,:,:,jpgsi)  &
             &                    + trn(:,:,:,jpdsi)                     ) * cvol(:,:,:)  )
-         ! 
+         !
+         silbudget = silbudget / areatot
+         CALL iom_put( "psiltot", silbudget )
+      ENDIF
+      !
+      IF( iom_use( "palktot" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
          alkbudget = glob_sum( (   trn(:,:,:,jpno3) * rno3              &
             &                    + trn(:,:,:,jptal)                     &
             &                    + trn(:,:,:,jpcal) * 2.                ) * cvol(:,:,:)  )
-         ! 
+         !
+         alkbudget = alkbudget / areatot
+         CALL iom_put( "palktot", alkbudget )
+      ENDIF
+      !
+      IF( iom_use( "pfertot" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
          ferbudget = glob_sum( (   trn(:,:,:,jpfer) + trn(:,:,:,jpnfe)  &
             &                    + trn(:,:,:,jpdfe)                     &
@@ -434,25 +457,24 @@
             &                    + trn(:,:,:,jpzoo) * ferat3            &
             &                    + trn(:,:,:,jpmes) * ferat3            ) * cvol(:,:,:)  )
-
-         !
+         !
+         ferbudget = ferbudget / areatot
+         CALL iom_put( "pfertot", ferbudget )
+      ENDIF
+      !
+      IF( ln_check_mass .AND. kt == nitend ) THEN   ! Compute the budget of NO3, ALK, Si, Fer
+         zfact = rfact2r * 12. / 1.e15 * ryyss    ! conversion molC/kt --> PgC/year
          t_atm_co2_flx  = t_atm_co2_flx / glob_sum( e1e2t(:,:) )
-         t_oce_co2_flx  = t_oce_co2_flx         * 12. / 1.e15 * (-1 )
-         tpp            = tpp           * 1000. * 12. / 1.E15
-         t_oce_co2_exp  = t_oce_co2_exp * 1000. * 12. / 1.E15
-         !
-         no3budget = no3budget / areatot
-         silbudget = silbudget / areatot
-         alkbudget = alkbudget / areatot
-         ferbudget = ferbudget / areatot
-         !
-         IF(lwp) THEN
-            WRITE(numco2,9000) ndastp, t_atm_co2_flx, t_oce_co2_flx, tpp, t_oce_co2_exp
-            WRITE(numnut,9500) ndastp, alkbudget, no3budget, silbudget, ferbudget
-         ENDIF
-         !
-      ENDIF
-       !
+         t_oce_co2_flx  = t_oce_co2_flx         * zfact * (-1 )
+         tpp            = tpp           * 1000. * zfact
+         t_oce_co2_exp  = t_oce_co2_exp * 1000. * zfact
+         IF( lwp ) WRITE(numco2,9000) ndastp, t_atm_co2_flx, t_oce_co2_flx, tpp, t_oce_co2_exp
+         IF( lwp ) WRITE(numnut,9500) ndastp, alkbudget        * 1.e+06, &
+             &                                no3budget * rno3 * 1.e+06, &
+             &                                silbudget        * 1.e+06, &
+             &                                ferbudget        * 1.e+09
+      ENDIF
+      !
  9000  FORMAT(i8,f10.5,e18.10,f10.5,f10.5)
- 9500  FORMAT(i8,4e18.10)     
+ 9500  FORMAT(i8,4e18.10)
        !
    END SUBROUTINE p4z_chk_mass
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90	(revision 4996)
@@ -36,4 +36,7 @@
    REAL(wp) ::   rfact2, rfact2r   !: ???
    REAL(wp) ::   xstep             !: Time step duration for biology
+   REAL(wp) ::   ryyss             !: number of seconds per year 
+   REAL(wp) ::   r1_ryyss          !: inverse number of seconds per year 
+
 
    !!*  Biological parameters 
@@ -53,4 +56,5 @@
    REAL(wp) ::  t_oce_co2_exp      !: total carbon export
    REAL(wp) ::  t_oce_co2_flx      !: Total ocean carbon flux
+   REAL(wp) ::  t_oce_co2_flx_cum  !: Cumulative Total ocean carbon flux
    REAL(wp) ::  t_atm_co2_flx      !: global mean of atmospheric pco2
 
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 4996)
@@ -100,4 +100,7 @@
       IF( lk_mpp    )   CALL mpp_sum( ierr )
       IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'pisces_alloc: unable to allocate PISCES arrays' )
+      !
+      ryyss    = nyear_len(1) * rday    ! number of seconds per year
+      r1_ryyss = 1. / ryyss
       !
 
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcwri_pisces.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcwri_pisces.F90	(revision 4995)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcwri_pisces.F90	(revision 4996)
@@ -30,7 +30,8 @@
       !! ** Purpose :   output passive tracers fields 
       !!---------------------------------------------------------------------
-      CHARACTER (len=20)   :: cltra
-      REAL(wp)             :: zrfact
-      INTEGER              :: jn
+      CHARACTER (len=20)           :: cltra
+      REAL(wp)                     :: zfact
+      INTEGER                      :: ji, jj, jk, jn
+      REAL(wp), DIMENSION(jpi,jpj) :: zdic, zo2min, zdepo2min
       !!---------------------------------------------------------------------
  
@@ -40,23 +41,43 @@
       DO jn = jp_pcs0, jp_pcs1
          cltra = TRIM( ctrcnm(jn) )                  ! short title for tracer
-         IF( lk_vvl ) THEN
-            CALL iom_put( cltra, trn(:,:,:,jn) * fse3t_n(:,:,:) )
-         ELSE
-            CALL iom_put( cltra, trn(:,:,:,jn) )
-         ENDIF
-         CALL iom_put( cltra, trn(:,:,:,jn) * zrfact )
+         CALL iom_put( cltra, trn(:,:,:,jn) )
       END DO
 #else
       DO jn = jp_pcs0, jp_pcs1
-         zrfact = 1.0e+6 
-         IF( jn == jpno3 .OR. jn == jpnh4 ) zrfact = rno3 * 1.0e+6 
-         IF( jn == jppo4  )                 zrfact = po4r * 1.0e+6
+         zfact = 1.0e+6 
+         IF( jn == jpno3 .OR. jn == jpnh4 ) zfact = rno3 * 1.0e+6 
+         IF( jn == jppo4  )                 zfact = po4r * 1.0e+6
          cltra = TRIM( ctrcnm(jn) )                  ! short title for tracer
-         IF( lk_vvl ) THEN
-            CALL iom_put( cltra, trn(:,:,:,jn) * fse3t_n(:,:,:) * zrfact )
-         ELSE
-            CALL iom_put( cltra, trn(:,:,:,jn) * zrfact )
-         ENDIF
+         IF( iom_use( cltra ) )  CALL iom_put( cltra, trn(:,:,:,jn) * zfact )
       END DO
+
+      IF( iom_use( "INTDIC" ) ) THEN                     !   DIC content in kg/m2
+         zdic(:,:) = 0.
+         DO jk = 1, jpkm1
+            zdic(:,:) = zdic(:,:) + trn(:,:,jk,jpdic) * fse3t(:,:,jk) * tmask(:,:,jk) * 12.
+         ENDDO
+         CALL iom_put( 'INTDIC', zdic )     
+      ENDIF
+      !
+      IF( iom_use( "O2MIN" ) .OR. iom_use ( "ZO2MIN" ) ) THEN  ! Oxygen minimum concentration and depth 
+         zo2min   (:,:) = trn(:,:,1,jpoxy) * tmask(:,:,1)
+         zdepo2min(:,:) = fsdepw(:,:,1)    * tmask(:,:,1)
+         DO jk = 2, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF( tmask(ji,jj,jk) == 1 ) then
+                     IF( trn(ji,jj,jk,jpoxy) < zo2min(ji,jj) ) then
+                        zo2min   (ji,jj) = trn(ji,jj,jk,jpoxy)
+                        zdepo2min(ji,jj) = fsdepw(ji,jj,jk)
+                     ENDIF
+                  ENDIF
+               END DO
+            END DO
+         END DO
+         !
+         CALL iom_put('O2MIN' , zo2min     )                              ! oxygen minimum concentration
+         CALL iom_put('ZO2MIN', zdepo2min  )                              ! depth of oxygen minimum concentration
+          !
+      ENDIF
 #endif
       !
