Index: trunk/NEMO/OFF_SRC/IOM/iom.F90
===================================================================
--- trunk/NEMO/OFF_SRC/IOM/iom.F90	(revision 1455)
+++ trunk/NEMO/OFF_SRC/IOM/iom.F90	(revision 1457)
@@ -33,4 +33,10 @@
    PUBLIC   !   must be public to be able to access iom_def through iom
    
+#if defined key_iomput
+   LOGICAL, PUBLIC, PARAMETER ::   lk_iomput = .TRUE.       !: iom_put flag
+#else
+   LOGICAL, PUBLIC, PARAMETER ::   lk_iomput = .FALSE.       !: iom_put flag
+#endif
+
    PUBLIC iom_init, iom_open, iom_close, iom_setkt, iom_varid, iom_get, iom_gettime, iom_rstput, iom_put
 
Index: trunk/NEMO/OFF_SRC/opa.F90
===================================================================
--- trunk/NEMO/OFF_SRC/opa.F90	(revision 1455)
+++ trunk/NEMO/OFF_SRC/opa.F90	(revision 1457)
@@ -225,8 +225,6 @@
       
       CALL dta_dyn( nit000 )                 ! Initialization for the dynamics
-     
+      
       CALL trc_ini                           ! Passive tracers
-
-      CALL day_init                          ! Calendar
       !                                     ! Ocean physics
       CALL tra_qsr_init                         ! Solar radiation penetration
Index: trunk/NEMO/OFF_SRC/step.F90
===================================================================
--- trunk/NEMO/OFF_SRC/step.F90	(revision 1455)
+++ trunk/NEMO/OFF_SRC/step.F90	(revision 1457)
@@ -22,4 +22,5 @@
    USE trcrst          ! restart for passive tracers
    USE stpctl          ! time stepping control            (stp_ctl routine)
+   USE iom
 
    IMPLICIT NONE
@@ -70,16 +71,13 @@
       !! * Arguments
       INTEGER, INTENT( in ) ::   kstp   ! ocean time-step index
-
-      !! * local declarations
-      INTEGER ::   indic    ! error indicator if < 0
       !! ---------------------------------------------------------------------
-
-      indic = 1                    ! reset to no error condition
 
       CALL day( kstp )             ! Calendar
 
-      CALL dta_dyn( kstp )          ! Interpolation of the dynamical fields
+      IF( lk_iomput ) CALL iom_setkt( kstp )       ! say to iom that we are at time step kstp
 
-      CALL trc_stp( kstp, indic)           ! time-stepping
+      CALL dta_dyn( kstp )         ! Interpolation of the dynamical fields
+
+      CALL trc_stp( kstp )         ! time-stepping
 
 
Index: trunk/NEMO/OPA_SRC/IOM/iom.F90
===================================================================
--- trunk/NEMO/OPA_SRC/IOM/iom.F90	(revision 1455)
+++ trunk/NEMO/OPA_SRC/IOM/iom.F90	(revision 1457)
@@ -33,4 +33,9 @@
    PUBLIC   !   must be public to be able to access iom_def through iom
    
+#if defined key_iomput
+   LOGICAL, PUBLIC, PARAMETER ::   lk_iomput = .TRUE.       !: iom_put flag
+#else
+   LOGICAL, PUBLIC, PARAMETER ::   lk_iomput = .FALSE.       !: iom_put flag
+#endif
    PUBLIC iom_init, iom_open, iom_close, iom_setkt, iom_varid, iom_get, iom_gettime, iom_rstput, iom_put
 
@@ -892,5 +897,5 @@
 
    !!----------------------------------------------------------------------
-   !!                   INTERFACE iom_rstput
+   !!                   INTERFACE iom_put
    !!----------------------------------------------------------------------
    SUBROUTINE iom_p2d( cdname, pfield2d )
Index: trunk/NEMO/OPA_SRC/step.F90
===================================================================
--- trunk/NEMO/OPA_SRC/step.F90	(revision 1455)
+++ trunk/NEMO/OPA_SRC/step.F90	(revision 1457)
@@ -281,5 +281,5 @@
       ! N.B. ua, va, ta, sa arrays are used as workspace in this section
       !-----------------------------------------------------------------------
-                             CALL trc_stp( kstp, indic )            ! time-stepping
+                             CALL trc_stp( kstp )            ! time-stepping
 #endif
 
Index: trunk/NEMO/TOP_SRC/LOBSTER/trcbio.F90
===================================================================
--- trunk/NEMO/TOP_SRC/LOBSTER/trcbio.F90	(revision 1455)
+++ trunk/NEMO/TOP_SRC/LOBSTER/trcbio.F90	(revision 1457)
@@ -22,4 +22,5 @@
    USE trdmld_trc
    USE trdmld_trc_oce
+   USE iom
    
    IMPLICIT NONE
@@ -29,5 +30,5 @@
 
    !!* Substitution
-#  include "domzgr_substitute.h90"
+#  include "top_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
@@ -71,5 +72,11 @@
       REAL(wp) ::   znh4a, zdeta, zdoma, zzoobod, zboddet, zdomaju
 #if defined key_trc_diaadd
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t
+      REAL(wp) ::   ze3t
+#endif
+#if defined key_trc_diaadd && defined key_iomput
+      REAL(wp), DIMENSION(jpi,jpj,17)     :: zw2d
+# if defined key_trc_dia3d
+      REAL(wp), DIMENSION(jpi,jpj,jpk,3) :: zw3d
+# endif
 #endif
       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   ztrbio
@@ -83,17 +90,15 @@
       ENDIF
 
-#if defined key_trc_diaadd
-      ! convert fluxes in per day
-      ze3t(:,:,:) = 0.e0
-      DO jk = 1, jpkbm1
-         ze3t(:,:,jk) = fse3t(:,:,jk) * 86400.
-      END DO 
-#endif
-
       fbod(:,:) = 0.e0
-#if defined key_trc_diaadd
+#if defined key_trc_diaadd && ! defined key_iomput
       DO jl = jp_lob0_2d, jp_lob1_2d
          trc2d(:,:,jl) = 0.e0
       END DO 
+#endif
+#if defined key_trc_diaadd && defined key_iomput
+      zw2d(:,:,:) = 0.e0
+# if defined key_trc_dia3d
+      zw3d(:,:,:,:) = 0.e0
+# endif
 #endif
 
@@ -107,6 +112,5 @@
          !                                   ! -------------------------- !
          DO jj = 2, jpjm1
-            DO ji = 2, jpim1           !!gm  use of fs_2 fs_jpm1  required here
-
+            DO ji = fs_2, fs_jpim1 
                ! trophic variables( det, zoo, phy, no3, nh4, dom)
                ! ------------------------------------------------
@@ -126,5 +130,4 @@
                zlno3 = zno3 * EXP( -psinut * znh4 ) / ( akno3 + zno3 )
                zlnh4 = znh4 / (znh4+aknh4) 
-
 
                ! sinks and sources
@@ -189,4 +192,12 @@
                zdeta =   zphydet + zzoodet - zdetzoo - zdetnh4 - zdetdom + zboddet
                zdoma =   zphydom + zzoodom + zdetdom - zdomnh4 - zdomaju
+
+               ! tracer flux at totox-point added to the general trend
+               tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) + zdeta
+               tra(ji,jj,jk,jpzoo) = tra(ji,jj,jk,jpzoo) + zzooa
+               tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) + zphya
+               tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) + zno3a
+               tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + znh4a
+               tra(ji,jj,jk,jpdom) = tra(ji,jj,jk,jpdom) + zdoma
 
 #if defined key_trc_diabio
@@ -206,30 +217,4 @@
                trbio(ji,jj,jk,jp_lob0_trd + 14) = zdetnh4
                trbio(ji,jj,jk,jp_lob0_trd + 15) = zdetdom
-#endif
-#if defined key_trc_diaadd
-               trc2d(ji,jj,jp_lob0_2d    ) = trc2d(ji,jj, jp_lob0_2d    ) + zno3phy * ze3t(ji,jj,jk) 
-               trc2d(ji,jj,jp_lob0_2d + 1) = trc2d(ji,jj, jp_lob0_2d + 1) + znh4phy * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d + 2) = trc2d(ji,jj, jp_lob0_2d + 2) + zphydom * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d + 3) = trc2d(ji,jj, jp_lob0_2d + 3) + zphynh4 * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d + 4) = trc2d(ji,jj, jp_lob0_2d + 4) + zphyzoo * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d + 5) = trc2d(ji,jj, jp_lob0_2d + 5) + zphydet * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d + 6) = trc2d(ji,jj, jp_lob0_2d + 6) + zdetzoo * ze3t(ji,jj,jk)
-! trend number 8 is in trcsed.F            
-               trc2d(ji,jj,jp_lob0_2d +  8) = trc2d(ji,jj,jp_lob0_2d +  8) + zzoodet * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d +  9) = trc2d(ji,jj,jp_lob0_2d +  9) + zzoobod * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d + 10) = trc2d(ji,jj,jp_lob0_2d + 10) + zzoonh4 * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d + 11) = trc2d(ji,jj,jp_lob0_2d + 11) + zzoodom * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d + 12) = trc2d(ji,jj,jp_lob0_2d + 12) + znh4no3 * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d + 13) = trc2d(ji,jj,jp_lob0_2d + 13) + zdomnh4 * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d + 14) = trc2d(ji,jj,jp_lob0_2d + 14) + zdetnh4 * ze3t(ji,jj,jk)             
-               trc2d(ji,jj,jp_lob0_2d + 15) = trc2d(ji,jj,jp_lob0_2d + 15) + (  zno3phy + znh4phy - zphynh4   &
-                  &                                 - zphydom - zphyzoo - zphydet ) * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d + 16) = trc2d(ji,jj,jp_lob0_2d + 16) + (  zphyzoo + zdetzoo - zzoodet   &
-                  &                                 - zzoobod - zzoonh4 - zzoodom ) * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d + 17) = trc2d(ji,jj,jp_lob0_2d + 17) + zdetdom * ze3t(ji,jj,jk)
-! trend number 19 is in trcexp.F
-               trc3d(ji,jj,jk,jp_lob0_3d    ) = zno3phy * 86400     
-               trc3d(ji,jj,jk,jp_lob0_3d + 1) = znh4phy * 86400     
-               trc3d(ji,jj,jk,jp_lob0_3d + 2) = znh4no3 * 86400     
 #endif
                IF( l_trdtrc ) THEN
@@ -253,27 +238,69 @@
                 ENDIF
 
-
-               ! tracer flux at totox-point added to the general trend
-               tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) + zdeta
-               tra(ji,jj,jk,jpzoo) = tra(ji,jj,jk,jpzoo) + zzooa
-               tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) + zphya
-               tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) + zno3a
-               tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + znh4a
-               tra(ji,jj,jk,jpdom) = tra(ji,jj,jk,jpdom) + zdoma
-
+#if defined key_trc_diaadd
+               ! convert fluxes in per day
+               ze3t = fse3t(ji,jj,jk) * 86400.
+#if ! defined key_iomput
+               trc2d(ji,jj,jp_lob0_2d    ) = trc2d(ji,jj, jp_lob0_2d    ) + zno3phy * ze3t 
+               trc2d(ji,jj,jp_lob0_2d + 1) = trc2d(ji,jj, jp_lob0_2d + 1) + znh4phy * ze3t
+               trc2d(ji,jj,jp_lob0_2d + 2) = trc2d(ji,jj, jp_lob0_2d + 2) + zphydom * ze3t
+               trc2d(ji,jj,jp_lob0_2d + 3) = trc2d(ji,jj, jp_lob0_2d + 3) + zphynh4 * ze3t
+               trc2d(ji,jj,jp_lob0_2d + 4) = trc2d(ji,jj, jp_lob0_2d + 4) + zphyzoo * ze3t
+               trc2d(ji,jj,jp_lob0_2d + 5) = trc2d(ji,jj, jp_lob0_2d + 5) + zphydet * ze3t
+               trc2d(ji,jj,jp_lob0_2d + 6) = trc2d(ji,jj, jp_lob0_2d + 6) + zdetzoo * ze3t
+               ! trend number 8 is in trcsed.F            
+               trc2d(ji,jj,jp_lob0_2d +  8) = trc2d(ji,jj,jp_lob0_2d +  8) + zzoodet * ze3t
+               trc2d(ji,jj,jp_lob0_2d +  9) = trc2d(ji,jj,jp_lob0_2d +  9) + zzoobod * ze3t
+               trc2d(ji,jj,jp_lob0_2d + 10) = trc2d(ji,jj,jp_lob0_2d + 10) + zzoonh4 * ze3t
+               trc2d(ji,jj,jp_lob0_2d + 11) = trc2d(ji,jj,jp_lob0_2d + 11) + zzoodom * ze3t
+               trc2d(ji,jj,jp_lob0_2d + 12) = trc2d(ji,jj,jp_lob0_2d + 12) + znh4no3 * ze3t
+               trc2d(ji,jj,jp_lob0_2d + 13) = trc2d(ji,jj,jp_lob0_2d + 13) + zdomnh4 * ze3t
+               trc2d(ji,jj,jp_lob0_2d + 14) = trc2d(ji,jj,jp_lob0_2d + 14) + zdetnh4 * ze3t             
+               trc2d(ji,jj,jp_lob0_2d + 15) = trc2d(ji,jj,jp_lob0_2d + 15) + (  zno3phy + znh4phy - zphynh4   &
+                  &                                 - zphydom - zphyzoo - zphydet ) * ze3t
+               trc2d(ji,jj,jp_lob0_2d + 16) = trc2d(ji,jj,jp_lob0_2d + 16) + (  zphyzoo + zdetzoo - zzoodet   &
+                  &                                 - zzoobod - zzoonh4 - zzoodom ) * ze3t
+               trc2d(ji,jj,jp_lob0_2d + 17) = trc2d(ji,jj,jp_lob0_2d + 17) + zdetdom * ze3t
+               ! trend number 19 is in trcexp.F
+#else
+               zw2d(ji,jj,1)  = zw2d(ji,jj,1)  + zno3phy * ze3t 
+               zw2d(ji,jj,2)  = zw2d(ji,jj,2)  + znh4phy * ze3t
+               zw2d(ji,jj,3)  = zw2d(ji,jj,3)  + zphydom * ze3t
+               zw2d(ji,jj,4)  = zw2d(ji,jj,4)  + zphynh4 * ze3t
+               zw2d(ji,jj,5)  = zw2d(ji,jj,5)  + zphyzoo * ze3t
+               zw2d(ji,jj,6)  = zw2d(ji,jj,6)  + zphydet * ze3t
+               zw2d(ji,jj,7)  = zw2d(ji,jj,7)  + zdetzoo * ze3t
+               zw2d(ji,jj,8)  = zw2d(ji,jj,8)  + zzoodet * ze3t
+               zw2d(ji,jj,9)  = zw2d(ji,jj,9)  + zzoobod * ze3t
+               zw2d(ji,jj,10) = zw2d(ji,jj,10) + zzoonh4 * ze3t
+               zw2d(ji,jj,11) = zw2d(ji,jj,11) + zzoodom * ze3t
+               zw2d(ji,jj,12) = zw2d(ji,jj,12) + znh4no3 * ze3t
+               zw2d(ji,jj,13) = zw2d(ji,jj,13) + zdomnh4 * ze3t
+               zw2d(ji,jj,14) = zw2d(ji,jj,14) + zdetnh4 * ze3t             
+               zw2d(ji,jj,15) = zw2d(ji,jj,15) + ( zno3phy + znh4phy - zphynh4 - zphydom - zphyzoo - zphydet ) * ze3t
+               zw2d(ji,jj,16) = zw2d(ji,jj,16) + ( zphyzoo + zdetzoo - zzoodet - zzoobod - zzoonh4 - zzoodom ) * ze3t
+               zw2d(ji,jj,17) = zw2d(ji,jj,17) + zdetdom * ze3t
+#endif
+#if defined key_trc_dia3d 
+# if ! defined key_iomput
+               trc3d(ji,jj,jk,jp_lob0_3d    ) = zno3phy * 86400     
+               trc3d(ji,jj,jk,jp_lob0_3d + 1) = znh4phy * 86400     
+               trc3d(ji,jj,jk,jp_lob0_3d + 2) = znh4no3 * 86400   
+# else
+               zw3d(ji,jj,jk,1) = zno3phy * 86400     
+               zw3d(ji,jj,jk,2) = znh4phy * 86400     
+               zw3d(ji,jj,jk,3) = znh4no3 * 86400   
+# endif
+#endif  
+#endif
             END DO
          END DO
       END DO
 
-      !
-
-!!gm do loop until jpkm1 only!
       !                                      ! -------------------------- !
-      DO jk = jpkb, jpk                      !  Upper ocean (bio-layers)  !
+      DO jk = jpkb, jpkm1                    !  Upper ocean (bio-layers)  !
          !                                   ! -------------------------- !
-
          DO jj = 2, jpjm1
-            DO ji = 2,jpim1             !!gm use of fs_2 & fs_jpim1 required
- 
+            DO ji = fs_2, fs_jpim1 
                ! remineralisation of all quantities towards nitrate 
 
@@ -335,4 +362,12 @@
                zdoma = zphydom + zzoodom + zdetdom - zdomnh4 - zdomaju
 
+               ! tracer flux at totox-point added to the general trend
+               tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) + zdeta
+               tra(ji,jj,jk,jpzoo) = tra(ji,jj,jk,jpzoo) + zzooa
+               tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) + zphya
+               tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) + zno3a
+               tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + znh4a
+               tra(ji,jj,jk,jpdom) = tra(ji,jj,jk,jpdom) + zdoma
+               !
 #if defined key_trc_diabio
                trbio(ji,jj,jk,jp_lob0_trd     ) = zno3phy
@@ -351,31 +386,4 @@
                trbio(ji,jj,jk,jp_lob0_trd + 14) = zdetnh4
                trbio(ji,jj,jk,jp_lob0_trd + 15) = zdetdom
-#endif
-#if defined key_trc_diaadd
-               trc2d(ji,jj, jp_lob0_2d    ) = trc2d(ji,jj, jp_lob0_2d    ) + zno3phy * ze3t(ji,jj,jk)          
-               trc2d(ji,jj, jp_lob0_2d + 1) = trc2d(ji,jj, jp_lob0_2d + 1) + znh4phy * ze3t(ji,jj,jk)
-               trc2d(ji,jj, jp_lob0_2d + 2) = trc2d(ji,jj, jp_lob0_2d + 2) + zphydom * ze3t(ji,jj,jk)
-               trc2d(ji,jj, jp_lob0_2d + 3) = trc2d(ji,jj, jp_lob0_2d + 3) + zphynh4 * ze3t(ji,jj,jk)
-               trc2d(ji,jj, jp_lob0_2d + 4) = trc2d(ji,jj, jp_lob0_2d + 4) + zphyzoo * ze3t(ji,jj,jk)
-               trc2d(ji,jj, jp_lob0_2d + 5) = trc2d(ji,jj, jp_lob0_2d + 5) + zphydet * ze3t(ji,jj,jk)
-               trc2d(ji,jj, jp_lob0_2d + 6) = trc2d(ji,jj, jp_lob0_2d + 6) + zdetzoo * ze3t(ji,jj,jk)
-               ! trend number 8 is in trcsed.F            
-               trc2d(ji,jj,jp_lob0_2d +  8) = trc2d(ji,jj,jp_lob0_2d +  8) + zzoodet * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d +  9) = trc2d(ji,jj,jp_lob0_2d +  9) + zzoobod * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d + 10) = trc2d(ji,jj,jp_lob0_2d + 10) + zzoonh4 * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d + 11) = trc2d(ji,jj,jp_lob0_2d + 11) + zzoodom * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d + 12) = trc2d(ji,jj,jp_lob0_2d + 12) + znh4no3 * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d + 13) = trc2d(ji,jj,jp_lob0_2d + 13) + zdomnh4 * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d + 14) = trc2d(ji,jj,jp_lob0_2d + 14) + zdetnh4 * ze3t(ji,jj,jk)
-             
-               trc2d(ji,jj,jp_lob0_2d + 15) = trc2d(ji,jj,jp_lob0_2d + 15) + (  zno3phy + znh4phy - zphynh4   &
-                  &                                 - zphydom - zphyzoo - zphydet  ) * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d + 16) = trc2d(ji,jj,jp_lob0_2d + 16) + (  zphyzoo + zdetzoo - zzoodet   &
-                  &                                 - zzoobod - zzoonh4 - zzoodom  ) * ze3t(ji,jj,jk)
-               trc2d(ji,jj,jp_lob0_2d + 17) = trc2d(ji,jj,jp_lob0_2d + 17) + zdetdom * ze3t(ji,jj,jk)
-
-               trc3d(ji,jj,jk,jp_lob0_3d    ) =  zno3phy * 86400     
-               trc3d(ji,jj,jk,jp_lob0_3d + 1) =  znh4phy * 86400     
-               trc3d(ji,jj,jk,jp_lob0_3d + 2) =  znh4no3 * 86400     
 #endif
                IF( l_trdtrc ) THEN
@@ -398,14 +406,15 @@
                   !  trend number 17 in trcexp
                 ENDIF
-
-
-               ! tracer flux at totox-point added to the general trend
-               tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) + zdeta
-               tra(ji,jj,jk,jpzoo) = tra(ji,jj,jk,jpzoo) + zzooa
-               tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) + zphya
-               tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) + zno3a
-               tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + znh4a
-               tra(ji,jj,jk,jpdom) = tra(ji,jj,jk,jpdom) + zdoma
-               !
+#if defined key_trc_diaadd && defined key_trc_dia3d
+# if ! defined key_iomput
+               trc3d(ji,jj,jk,jp_lob0_3d    ) =  zno3phy * 86400     
+               trc3d(ji,jj,jk,jp_lob0_3d + 1) =  znh4phy * 86400     
+               trc3d(ji,jj,jk,jp_lob0_3d + 2) =  znh4no3 * 86400     
+# else
+               zw3d(ji,jj,jk,1) = zno3phy * 86400     
+               zw3d(ji,jj,jk,2) = znh4phy * 86400     
+               zw3d(ji,jj,jk,3) = znh4no3 * 86400   
+# endif
+#endif
             END DO
          END DO
@@ -413,11 +422,49 @@
 
 #if defined key_trc_diaadd
-      ! Lateral boundary conditions on trc2d and trc3d
+      ! Lateral boundary conditions 
+# if ! defined key_iomput
       DO jl = jp_lob0_2d, jp_lob1_2d
           CALL lbc_lnk( trc2d(:,:,jl),'T', 1. )
       END DO 
+# else
+      DO jl = 1, 17 
+          CALL lbc_lnk( zw2d(:,:,jl),'T', 1. )
+      END DO
+      ! 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( "TDETDOM", zw2d(:,:,17) )
+# endif
+#endif
+
+#if defined key_trc_diaadd && defined key_trc_dia3d
+      ! Lateral boundary conditions 
+# if ! defined key_iomput
       DO jl = jp_lob0_3d, jp_lob1_3d
           CALL lbc_lnk( trc3d(:,:,1,jl),'T', 1. )
       END DO 
+# else
+      DO jl = 1, 3
+          CALL lbc_lnk( zw3d(:,:,:,jl),'T', 1. )
+      END DO
+      ! save diagnostics
+      CALL iom_put( "FNO3PHY", zw3d(:,:,:,1) )
+      CALL iom_put( "FNH4PHY", zw3d(:,:,:,2) )
+      CALL iom_put( "FNH4NO3", zw3d(:,:,:,3) )
+# endif 
 #endif
 
Index: trunk/NEMO/TOP_SRC/LOBSTER/trcexp.F90
===================================================================
--- trunk/NEMO/TOP_SRC/LOBSTER/trcexp.F90	(revision 1455)
+++ trunk/NEMO/TOP_SRC/LOBSTER/trcexp.F90	(revision 1457)
@@ -23,4 +23,5 @@
    USE trdmld_trc
    USE trdmld_trc_oce
+   USE iom
 
    IMPLICIT NONE
@@ -30,5 +31,5 @@
 
    !!* Substitution
-#  include "domzgr_substitute.h90"
+#  include "top_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
@@ -54,8 +55,6 @@
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index      
       !!
-      INTEGER  ::   ji, jj, jk, jl
-      REAL(wp) ::   zgeolpoc, zfact
-      INTEGER , DIMENSION(jpi,jpj) ::   ikbot
-      REAL(wp), DIMENSION(jpi,jpj) ::   zwork
+      INTEGER  ::   ji, jj, jk, jl, ikbot
+      REAL(wp) ::   zgeolpoc, zfact, zwork, ze3t
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrbio
       CHARACTER (len=25) :: charout
@@ -81,49 +80,31 @@
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3)   &
-     &                             + (1./fse3t(ji,jj,jk)) * dmin3(ji,jj,jk) * fbod(ji,jj)
+            DO ji = fs_2, fs_jpim1
+               ze3t = 1. / fse3t(ji,jj,jk)
+               tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) + ze3t * dmin3(ji,jj,jk) * fbod(ji,jj)
             END DO
          END DO
       END DO
 
-!     Find the last level of the water column
-!     Compute fluxes due to sinking particles (slow)
+      ! Find the last level of the water column
+      ! Compute fluxes due to sinking particles (slow)
    
-      ikbot(:,:) = jpk
-      zwork(:,:) = 0.e0
 
-!!gm ikbot already exist in opa...
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               IF( tmask(ji,jj,jk) == 1 .AND.  tmask(ji,jj,jk+1) == 0 ) THEN
-                  ikbot(ji,jj) = jk
-                  zwork(ji,jj) = vsed * trn(ji,jj,jk,jpdet)
-               ENDIF
-            END DO
+      zgeolpoc = 0.e0         !     Initialization
+      ! Release of nutrients from the "simple" sediment
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1
+            ikbot = mbathy(ji,jj) - 1
+            tra(ji,jj,ikbot,jpno3) = tra(ji,jj,ikbot,jpno3) + sedlam * sedpocn(ji,jj) / fse3t(ji,jj,ikbot) 
+            ! Deposition of organic matter in the sediment
+            zwork = vsed * trn(ji,jj,ikbot,jpdet)
+            sedpoca(ji,jj) = ( zwork + dminl(ji,jj) * fbod(ji,jj)   &
+               &           - sedlam * sedpocn(ji,jj) - sedlostpoc * sedpocn(ji,jj) ) * rdt
+            zgeolpoc = zgeolpoc + sedlostpoc * sedpocn(ji,jj) * e1t(ji,jj) * e2t(ji,jj)
          END DO
       END DO
 
-      zgeolpoc = 0.e0         !     Initialization
-
-      ! Release of nutrients from the "simple" sediment
       DO jj = 2, jpjm1
-         DO ji = 2, jpim1
-            tra(ji,jj,ikbot(ji,jj),jpno3) = tra(ji,jj,ikbot(ji,jj),jpno3)   &
-               &                          + sedlam * sedpocn(ji,jj) / fse3t(ji,jj,ikbot(ji,jj))
-
-            !     Deposition of organic matter in the sediment
-            zgeolpoc = zgeolpoc + sedlostpoc * sedpocn(ji,jj) * e1t(ji,jj) * e2t(ji,jj)
-
-!!gm factorisationof rdt just bellow...
-            sedpoca(ji,jj) = zwork(ji,jj) * rdt + dminl(ji,jj) * fbod(ji,jj) * rdt   &
-               &           - sedlam * sedpocn(ji,jj) * rdt - sedlostpoc * sedpocn(ji,jj) * rdt
-
-         END DO
-      END DO
-
-      DO jj = 2,jpjm1
-         DO ji = 2,jpim1
+         DO ji = fs_2, fs_jpim1
             tra(ji,jj,1,jpno3) = tra(ji,jj,1,jpno3) + zgeolpoc * cmask(ji,jj) / areacot / fse3t(ji,jj,1)
          END DO
@@ -133,7 +114,11 @@
  
       ! Oa & Ek: diagnostics depending on jpdia2d !          left as example
-# if defined key_trc_diaadd
+#if defined key_trc_diaadd
+# if ! defined key_iomput
       trc2d(:,:,jp_lob0_2d + 18) = sedpocn(:,:)
+# else
+     CALL iom_put( "SEDPOC" , sedpocn )
 # endif
+#endif
 
       ! Leap-frog scheme (only in explicit case, otherwise the 
Index: trunk/NEMO/TOP_SRC/LOBSTER/trclsm_lobster.F90
===================================================================
--- trunk/NEMO/TOP_SRC/LOBSTER/trclsm_lobster.F90	(revision 1455)
+++ trunk/NEMO/TOP_SRC/LOBSTER/trclsm_lobster.F90	(revision 1457)
@@ -42,5 +42,5 @@
       INTEGER ::   numnatl
       !!
-#if defined key_trc_diaadd
+#if defined key_trc_diaadd && ! defined key_iomput
       INTEGER :: jl, jn
       ! definition of additional diagnostic as a structure
@@ -78,5 +78,5 @@
 
       NAMELIST/namlobopt/ xkg0, xkr0, xkgp, xkrp, xlg, xlr, rpig
-#if defined key_trc_diaadd
+#if defined key_trc_diaadd && ! defined key_iomput
       NAMELIST/namlobdia/nwritedia, lobdia3d, lobdia2d     ! additional diagnostics
 #endif
@@ -280,5 +280,5 @@
       ENDIF
 
-#if defined key_trc_diaadd
+#if defined key_trc_diaadd && ! defined key_iomput
 
       ! Namelist namlobdia
Index: trunk/NEMO/TOP_SRC/LOBSTER/trcsed.F90
===================================================================
--- trunk/NEMO/TOP_SRC/LOBSTER/trcsed.F90	(revision 1455)
+++ trunk/NEMO/TOP_SRC/LOBSTER/trcsed.F90	(revision 1457)
@@ -20,5 +20,5 @@
    USE trdmld_trc
    USE trdmld_trc_oce
-
+   USE iom
    USE prtctl_trc      ! Print control for debbuging
 
@@ -61,4 +61,7 @@
       REAL(wp) ::   ztra
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwork
+#if defined key_trc_diaadd && defined key_iomput
+      REAL(wp), DIMENSION(jpi,jpj) ::  zw2d
+#endif
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrbio
       CHARACTER (len=25) :: charout
@@ -75,7 +78,10 @@
 
       ! for detritus sedimentation only - jpdet
-
       zwork(:,:,1  ) = 0.e0      ! surface value set to zero
       zwork(:,:,jpk) = 0.e0      ! bottom value  set to zero
+
+#if defined key_trc_diaadd && defined key_iomput
+      zw2d(:,:) = 0.
+# endif
 
       IF( l_trdtrc )THEN
@@ -85,11 +91,9 @@
 
       ! tracer flux at w-point: we use -vsed (downward flux)  with simplification : no e1*e2
-
-      DO jk = 2, jpk
+      DO jk = 2, jpkm1
          zwork(:,:,jk) = -vsed * trn(:,:,jk-1,jpdet)
       END DO
 
       ! tracer flux divergence at t-point added to the general trend
-
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -97,10 +101,14 @@
                ztra  = - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
                tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) + ztra
-# if defined key_trc_diabio
+#if defined key_trc_diabio
                trbio(ji,jj,jk,jp_lob0_trd + 7) = ztra
+#endif
+#if defined key_trc_diaadd
+# if ! defined key_iomput
+               trc2d(ji,jj,jp_lob0_2d + 7) = trc2d(ji,jj,jp_lob0_2d + 7) + ztra * fse3t(ji,jj,jk) * 86400.
+# else
+               zw2d(ji,jj) = zw2d(ji,jj) + ztra * fse3t(ji,jj,jk) * 86400.
 # endif
-# if defined key_trc_diaadd
-               trc2d(ji,jj,jp_lob0_2d + 7) = trc2d(ji,jj,jp_lob0_2d + 7) + ztra * fse3t(ji,jj,jk) * 86400.
-# endif
+#endif
             END DO
          END DO
@@ -112,6 +120,11 @@
 #endif
 #if defined key_trc_diaadd
+# if ! defined key_iomput
       jl = jp_lob0_2d + 7
       CALL lbc_lnk( trc2d(:,:,jl), 'T', 1. )      ! Lateral boundary conditions on trc2d
+# else
+      CALL lbc_lnk( zw2d(:,:), 'T', 1. )      ! Lateral boundary conditions on zw2d
+      CALL iom_put( "TDETSED", zw2d )
+# endif
 #endif
       !
Index: trunk/NEMO/TOP_SRC/PISCES/p4zbio.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zbio.F90	(revision 1455)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zbio.F90	(revision 1457)
@@ -26,5 +26,6 @@
    USE p4zrem          ! 
    USE prtctl_trc
-   
+   USE iom
+  
    IMPLICIT NONE
    PRIVATE
@@ -62,4 +63,7 @@
 #if defined key_kriest
       REAL(wp) ::  zcoef1, zcoef2
+#endif
+#if defined key_trc_dia3d && defined key_kriest && defined key_iomput
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d
 #endif
       CHARACTER (len=25) :: charout
@@ -128,7 +132,11 @@
 
 # if defined key_trc_dia3d && defined key_kriest
-      trc3d(:,:,:,jp_pcs0_3d + 10) = tra(:,:,:,jpcal) &
-     &                              * xnegtr(:,:,:) * 1.e3 * rfact2r  * tmask(:,:,:)
+#if ! defined key_iomput
+      trc3d(:,:,:,jp_pcs0_3d + 10) = tra(:,:,:,jpcal) * xnegtr(:,:,:) * 1.e3 * rfact2r  * tmask(:,:,:)
+#else
+      zw3d(:,:,:) = tra(:,:,:,jpcal) * xnegtr(:,:,:) * 1.e3 * rfact2r  * tmask(:,:,:)
+      IF( jnt == nrdttrc ) CALL iom_put( "PBSi", zw3d )
 # endif
+#endif
       !
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
Index: trunk/NEMO/TOP_SRC/PISCES/p4zflx.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zflx.F90	(revision 1455)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zflx.F90	(revision 1457)
@@ -23,4 +23,5 @@
    USE prtctl_trc
    USE p4zche
+   USE iom
 
    USE lib_mpp
@@ -32,9 +33,10 @@
 
    REAL(wp) :: &  ! pre-industrial atmospheric [co2] (ppm) 	
-     atcox  = 0.20946 ,   &  !:
-     atcco2 = 278.           !:
-
-   REAL(wp) :: &
-     tco2flx                 !: Total flux of carbon per year
+      atcox  = 0.20946 ,    &  !:
+      atcco2 = 278.            !:
+ 
+   REAL(wp) :: &          
+      xconv  = 0.01/3600.,  &  ! coefficients for conversion       
+      tco2flx                 !: Total flux of carbon per year
 
    !!* Substitution
@@ -58,8 +60,11 @@
       INTEGER, INTENT(in) :: kt
       INTEGER  ::   ji, jj, jrorr, nspyr
-      REAL(wp) ::   zttc, zws, zkgwan
+      REAL(wp) ::   ztc, ztc2, ztc3, zws, zkgwan
       REAL(wp) ::   zfld, zflu, zfld16, zflu16, zfact
-      REAL(wp) ::   zph, zah2, zbot, zdic, zalk, zschmitto2, zalka, zschmittco2
+      REAL(wp) ::   zph, zah2, zbot, zdic, zalk, zsch_o2, zalka, zsch_co2
       REAL(wp), DIMENSION(jpi,jpj) ::   zkgco2, zkgo2, zh2co3, ztco2flx
+#if defined key_trc_diaadd && defined key_iomput
+      REAL(wp), DIMENSION(jpi,jpj) ::  zcflx, zoflx, zkg, zdelc
+#endif
       CHARACTER (len=25) :: charout
 
@@ -120,34 +125,26 @@
 !CDIR NOVERRCHK
          DO ji = 1, jpi
-
-            zttc = MIN( 35., tn(ji,jj,1) )
-
-      ! Compute the schmidt Number both O2 and CO2
-      ! ------------------------------------------
-
-            zschmittco2 = 2073.1 - 125.62 * zttc + 3.6276 * zttc**2 - 0.043126 * zttc**3
-            zschmitto2  = 1953.4 - 128.0 * zttc + 3.9918 * zttc**2 - 0.050091 * zttc**3
-
+            ztc  = MIN( 35., tn(ji,jj,1) )
+            ztc2 = ztc * ztc
+            ztc3 = ztc * ztc2 
+            ! Compute the schmidt Number both O2 and CO2
+            zsch_co2 = 2073.1 - 125.62 * ztc + 3.6276 * ztc2 - 0.043126 * ztc3
+            zsch_o2  = 1953.4 - 128.0  * ztc + 3.9918 * ztc2 - 0.050091 * ztc3
             !  wind speed 
-            zws  = wndm(ji,jj)
-
-      ! Compute the piston velocity for O2 and CO2
-      ! ------------------------------------------
-
-            zkgwan = (  0.3 * zws * zws    &
-               &   + 2.5 * ( 0.5246 + zttc * ( 0.016256+zttc*0.00049946 ) ) ) &
+            zws  = wndm(ji,jj) * wndm(ji,jj)
+            ! Compute the piston velocity for O2 and CO2
+            zkgwan = 0.3 * zws  + 2.5 * ( 0.5246 + 0.016256 * ztc + 0.00049946  * ztc2 )
 # if defined key_off_degrad
-               &         * facvol(ji,jj,1)      &
-# endif
-               &   / (100. * 3600.)* ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1)  
-
-            ! COMPUTE GAS EXCHANGE FOR CO2
-            zkgco2(ji,jj) = zkgwan * SQRT( 660./ zschmittco2 )
-            zkgo2(ji,jj)  = zkgwan * SQRT( 660./ zschmitto2 )
-
-         END DO
-      END DO
-
-      ztco2flx(:,:) = 0.
+            zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1) * facvol(ji,jj,1)
+#else
+            zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1)
+#endif 
+            ! compute gas exchange for CO2 and O2
+            zkgco2(ji,jj) = zkgwan * SQRT( 660./ zsch_co2 )
+            zkgo2 (ji,jj) = zkgwan * SQRT( 660./ zsch_o2 )
+         END DO
+      END DO
+
+#if ! defined key_iomput
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -167,5 +164,5 @@
             tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + ( zfld16 - zflu16 ) / fse3t(ji,jj,1)
 
-# if defined key_trc_diaadd
+# if defined key_trc_diaadd 
             ! Save diagnostics
             trc2d(ji,jj,jp_pcs0_2d    ) = ( zfld - zflu ) * 1000.  * tmask(ji,jj,1)
@@ -177,4 +174,30 @@
       END DO
       !
+#else
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ! Compute CO2 flux for the sea and air
+            zfld = atcco2 * tmask(ji,jj,1) * chemc(ji,jj,1) * zkgco2(ji,jj)
+            zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj)
+            tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) / fse3t(ji,jj,1)
+
+            ! compute flux of carbon
+            ztco2flx(ji,jj) = ( zfld - zflu ) * rfact &
+               &             * e1t(ji,jj) * e2t(ji,jj) * tmask(ji,jj,1) * 1000.
+
+            ! Compute O2 flux 
+            zfld16 = atcox * chemc(ji,jj,2) *tmask(ji,jj,1) * zkgo2(ji,jj)
+            zflu16 = trn(ji,jj,1,jpoxy) * tmask(ji,jj,1) * zkgo2(ji,jj)
+            tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + ( zfld16 - zflu16 ) / fse3t(ji,jj,1)
+# if defined key_trc_diaadd 
+            ! Save diagnostics
+            zcflx(ji,jj) = ( zfld - zflu ) * 1000.  * tmask(ji,jj,1)
+            zoflx(ji,jj) = ( zfld16 - zflu16 ) * 1000. * tmask(ji,jj,1)
+            zkg  (ji,jj) = zkgco2(ji,jj) * tmask(ji,jj,1)
+            zdelc(ji,jj) = atcco2 - zh2co3(ji,jj) / ( chemc(ji,jj,1) + rtrn ) * tmask(ji,jj,1)
+# endif
+         END DO
+      END DO
+#endif
 
       ! Total Flux of Carbon
@@ -203,4 +226,10 @@
       ENDIF
 
+# if defined key_trc_diaadd && defined key_iomput
+      CALL iom_put( "Cflx", zcflx )
+      CALL iom_put( "Oflx", zoflx )
+      CALL iom_put( "Kg"  , zkg   )
+      CALL iom_put( "DelC", zdelc )
+#endif
 
    END SUBROUTINE p4z_flx
Index: trunk/NEMO/TOP_SRC/PISCES/p4zlys.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zlys.F90	(revision 1455)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zlys.F90	(revision 1457)
@@ -22,4 +22,5 @@
    USE sms_pisces
    USE prtctl_trc
+   USE iom
 
    IMPLICIT NONE
@@ -61,4 +62,7 @@
       REAL(wp) ::   zomegaca, zexcess, zexcess0
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zco3
+#if defined key_trc_dia3d && defined key_iomput
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d
+#endif
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
@@ -151,13 +155,16 @@
 
 # if defined key_trc_diaadd &&  defined key_trc_dia3d
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               trc3d(ji,jj,jk,jp_pcs0_3d    ) = hi  (ji,jj,jk)          * tmask(ji,jj,jk)
-               trc3d(ji,jj,jk,jp_pcs0_3d + 1) = zco3(ji,jj,jk)          * tmask(ji,jj,jk)
-               trc3d(ji,jj,jk,jp_pcs0_3d + 2) = aksp(ji,jj,jk) / calcon * tmask(ji,jj,jk)
-            ENDDO
-         ENDDO
-      ENDDO
+#  if ! defined key_iomput
+      trc3d(:,:,:,jp_pcs0_3d    ) = hi  (:,:,:)          * tmask(:,:,:)
+      trc3d(:,:,:,jp_pcs0_3d + 1) = zco3(:,:,:)          * tmask(:,:,:)
+      trc3d(:,:,:,jp_pcs0_3d + 2) = aksp(:,:,:) / calcon * tmask(:,:,:)
+#  else
+      zw3d(:,:,:) = hi  (:,:,:)          * tmask(:,:,:)
+      CALL iom_put( "PH", zw3d )
+      zw3d(:,:,:) = zco3(:,:,:)          * tmask(:,:,:)
+      CALL iom_put( "CO3", zw3d )
+      zw3d(:,:,:) = aksp(:,:,:) / calcon * tmask(:,:,:)
+      CALL iom_put( "CO3sat", zw3d )
+#  endif
 # endif
       !
Index: trunk/NEMO/TOP_SRC/PISCES/p4zopt.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zopt.F90	(revision 1455)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zopt.F90	(revision 1457)
@@ -18,4 +18,5 @@
    USE trc_oce        ! ocean-tracer share variables
    USE sms_pisces     ! Source Minus Sink of PISCES
+   USE iom
 
    IMPLICIT NONE
@@ -60,4 +61,8 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zekg, zekr, zekb
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze1 , ze2 , ze3, ze0
+#if defined key_trc_diaadd && defined key_iomput
+     REAL(wp), DIMENSION(jpi,jpj)      ::   zw2d
+     REAL(wp), DIMENSION(jpi,jpj,jpk)  ::   zw3d
+#endif
       !!---------------------------------------------------------------------
 
@@ -226,7 +231,17 @@
       END DO
 
-# if defined key_trc_diaadd
-      trc2d(:,:,  jp_pcs0_2d + 10) = heup (:,:) * tmask(:,:,1)      ! save for outputs
+#if defined key_trc_diaadd
+# if ! defined key_iomput
+      ! save for outputs
+      trc2d(:,:,  jp_pcs0_2d + 10) = heup(:,:  ) * tmask(:,:,1)  
+      trc3d(:,:,:,jp_pcs0_3d + 3)  = etot(:,:,:) * tmask(:,:,:)
+# else
+      ! write diagnostics 
+      zw2d(:,:  ) =  heup(:,:  ) * tmask(:,:,1)
+      zw3d(:,:,:) =  etot(:,:,:) * tmask(:,:,:)
+      IF( jnt == nrdttrc ) CALL iom_put( "heup", zw2d )               
+      IF( jnt == nrdttrc ) CALL iom_put( "PAR" , zw3d )
 # endif
+#endif
       !
    END SUBROUTINE p4z_opt
Index: trunk/NEMO/TOP_SRC/PISCES/p4zprod.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zprod.F90	(revision 1455)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zprod.F90	(revision 1457)
@@ -20,4 +20,5 @@
    USE p4zint
    USE p4zlim
+   USE iom
 
    USE lib_mpp
@@ -78,4 +79,7 @@
 #if defined key_trc_diaadd && defined key_trc_dia3d
       REAL(wp) ::   zrfact2
+#if  defined key_iomput
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d
+#endif
 #endif
       REAL(wp), DIMENSION(jpi,jpj)     ::   zmixnano   , zmixdiat, zstrn
@@ -107,6 +111,5 @@
 
 
-!     Computation of the optimal production
-!     -------------------------------------
+      ! Computation of the optimal production
 
 # if defined key_off_degrad
@@ -117,5 +120,4 @@
 
       ! compute the day length depending on latitude and the day
-      !--------------------------------------------------------
       IF(lwp) write(numout,*)
       IF(lwp) write(numout,*) 'p4zday : - Julian day ', nday_year
@@ -147,7 +149,6 @@
             DO ji = 1, jpi
 
-!      Computation of the P-I slope for nanos and diatoms
-!      --------------------------------------------------
-                IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+               ! Computation of the P-I slope for nanos and diatoms
+               IF( etot(ji,jj,jk) > 1.E-3 ) THEN
                    ztn    = MAX( 0., tn(ji,jj,jk) - 15. )
                    zadap  = 0.+ 1.* ztn / ( 2.+ ztn )
@@ -167,7 +168,5 @@
                      &          / ( prmax(ji,jj,jk) * rjjss * xlimdia(ji,jj,jk) + rtrn )
 
-!     Computation of production function
-!     ----------------------------------
-
+                   ! Computation of production function
                    zprbio(ji,jj,jk) = prmax(ji,jj,jk) * &
                      &                (  1.- EXP( -zpislopen * enano(ji,jj,jk) )  )
@@ -185,9 +184,9 @@
 
                 IF( etot(ji,jj,jk) > 1.E-3 ) THEN
-!    Si/C of diatoms
-!    ------------------------
-!    Si/C increases with iron stress and silicate availability
-!    Si/C is arbitrariliy increased for very high Si concentrations
-!    to mimic the very high ratios observed in the Southern Ocean (silpot2)
+                   !    Si/C of diatoms
+                   !    ------------------------
+                   !    Si/C increases with iron stress and silicate availability
+                   !    Si/C is arbitrariliy increased for very high Si concentrations
+                   !    to mimic the very high ratios observed in the Southern Ocean (silpot2)
 
                   zlim1  = trn(ji,jj,jk,jpsil) / ( trn(ji,jj,jk,jpsil) + xksi1 )
@@ -209,8 +208,6 @@
       END DO
 
-!    Computation of the limitation term due to
-!    A mixed layer deeper than the euphotic depth
-!    --------------------------------------------
-
+      !  Computation of the limitation term due to
+      !  A mixed layer deeper than the euphotic depth
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -221,12 +218,10 @@
          END DO
       END DO
-                                                                                
+ 
+      !  Mixed-layer effect on production                                                                               
       DO jk = 1, jpkm1
          DO jj = 1, jpj
             DO ji = 1, jpi
                IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
-
-!     Mixed-layer effect on production
-!     --------------------------------
                   zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * zmixnano(ji,jj)
                   zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * zmixdiat(ji,jj)
@@ -365,21 +360,31 @@
       zrfact2 = 1.e3 * rfact2r
       !   Supplementary diagnostics
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi 
-               trc3d(ji,jj,jk,jp_pcs0_3d + 3)  = etot(ji,jj,jk)               * tmask(ji,jj,jk)
-               trc3d(ji,jj,jk,jp_pcs0_3d + 4)  = zprorca (ji,jj,jk) * zrfact2 * tmask(ji,jj,jk)
-               trc3d(ji,jj,jk,jp_pcs0_3d + 5)  = zprorcad(ji,jj,jk) * zrfact2 * tmask(ji,jj,jk)
-               trc3d(ji,jj,jk,jp_pcs0_3d + 6)  = zpronew (ji,jj,jk) * zrfact2 * tmask(ji,jj,jk)
-               trc3d(ji,jj,jk,jp_pcs0_3d + 7)  = zpronewd(ji,jj,jk) * zrfact2 * tmask(ji,jj,jk)
-               trc3d(ji,jj,jk,jp_pcs0_3d + 8)  = zprorcad(ji,jj,jk) * zrfact2 * tmask(ji,jj,jk) &
-     &                                         * zysopt(ji,jj,jk)
-               trc3d(ji,jj,jk,jp_pcs0_3d + 9)  = zprofed (ji,jj,jk) * zrfact2 * tmask(ji,jj,jk)
+#  if ! defined key_iomput
+      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(:,:,:)
 #if ! defined key_kriest
-               trc3d(ji,jj,jk,jp_pcs0_3d + 10) = zprofen (ji,jj,jk) * zrfact2 * tmask(ji,jj,jk)
+      trc3d(:,:,:,jp_pcs0_3d + 10) = zprofen (:,:,:) * zrfact2 * tmask(:,:,:)
 #endif
-            ENDDO
-         ENDDO
-      ENDDO
+
+# else
+      zw3d(:,:,:) = zprorca (:,:,:) * zrfact2 * tmask(:,:,:)
+      IF( jnt == nrdttrc ) CALL iom_put( "PPPHY" , zw3d )
+      zw3d(:,:,:) = zprorcad(:,:,:) * zrfact2 * tmask(:,:,:)
+      IF( jnt == nrdttrc ) CALL iom_put( "PPPHY2", zw3d )
+      zw3d(:,:,:) = zpronew (:,:,:) * zrfact2 * tmask(:,:,:)
+      IF( jnt == nrdttrc ) CALL iom_put( "PPZOO" , zw3d )
+      zw3d(:,:,:) = zpronewd(:,:,:) * zrfact2 * tmask(:,:,:)
+      IF( jnt == nrdttrc ) CALL iom_put( "PPZOO2", zw3d )
+      zw3d(:,:,:) = zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) * zysopt(:,:,:)
+      IF( jnt == nrdttrc ) CALL iom_put( "PBSi"  , zw3d )
+      zw3d(:,:,:) = zprofed (:,:,:) * zrfact2 * tmask(:,:,:)
+      IF( jnt == nrdttrc ) CALL iom_put( "PFeD"  , zw3d )
+      zw3d(:,:,:) = zprofen (:,:,:) * zrfact2 * tmask(:,:,:)
+      IF( jnt == nrdttrc ) CALL iom_put( "PFeN"  , zw3d )
+# endif
 #endif
 
Index: trunk/NEMO/TOP_SRC/PISCES/p4zsed.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zsed.F90	(revision 1455)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zsed.F90	(revision 1457)
@@ -89,4 +89,7 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   zsidep
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   znitrpot, zirondep
+#if defined key_trc_diaadd && defined key_iomput
+     REAL(wp), DIMENSION(jpi,jpj)     ::   zw2d
+#endif
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
@@ -94,5 +97,4 @@
 
       IF( ( kt * jnt ) == nittrc000  )   CALL p4z_sed_init      ! Initialization (first time-step only)
-
       IF( (jnt == 1) .and. (bdustfer) )  CALL p4z_sbc( kt )
 
@@ -136,14 +138,10 @@
 
       DO jk = 1, jpkm1
-         trn(:,:,jk,jpfer) = trn(:,:,jk,jpfer)   &
-         &       + zirondep(:,:,jk) + ironsed(:,:,jk) * rfact2
+         trn(:,:,jk,jpfer) = trn(:,:,jk,jpfer) + zirondep(:,:,jk) + ironsed(:,:,jk) * rfact2
       END DO
 
 
 #if ! defined key_sed
-
       ! Initialisation of variables used to compute Sinking Speed
-      ! ---------------------------------------------------------
-
       zsumsedsi  = 0.e0
       zsumsedpo4 = 0.e0
@@ -154,10 +152,8 @@
       ! The factor for calcite comes from the alkalinity effect
       ! -------------------------------------------------------------
-
       DO jj = 1, jpj
          DO ji = 1, jpi
             ikt = MAX( mbathy(ji,jj)-1, 1 )
             zfact = e1t(ji,jj) * e2t(ji,jj) / rjjss * tmask_i(ji,jj)
-
 # if defined key_kriest
             zsumsedsi  = zsumsedsi  + zfact *  trn(ji,jj,ikt,jpdsi) * wscal (ji,jj,ikt)
@@ -168,7 +164,5 @@
                &       + trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt) )
 # endif
-
             zsumsedcal = zsumsedcal + zfact *  trn(ji,jj,ikt,jpcal) * wscal (ji,jj,ikt) * 2.e0
-
          END DO
       END DO
@@ -197,5 +191,4 @@
      &             * wsbio4(ji,jj,ikt)
 # endif
-
             trn(ji,jj,ikt,jpdsi) = trn(ji,jj,ikt,jpdsi) - zconctmp
 
@@ -226,5 +219,4 @@
             ikt = MAX( mbathy(ji,jj) - 1, 1 )
             zfact = zstep / fse3t(ji,jj,ikt)
-
 # if ! defined key_kriest
             zconctmp  = trn(ji,jj,ikt,jpgoc)
@@ -241,5 +233,4 @@
 
 # else
-
             zconctmp  = trn(ji,jj,ikt,jpnum)
             zconctmp2 = trn(ji,jj,ikt,jppoc)
@@ -327,13 +318,16 @@
       END DO
 
-# if defined key_trc_diaadd
-      DO jj = 1,jpj
-         DO ji = 1,jpi
-            trc2d(ji,jj,jp_pcs0_2d + 11) = zirondep(ji,jj,1) * 1.e+3 * rfact2r &
-      &                                   * fse3t(ji,jj,1) * tmask(ji,jj,1)
-            trc2d(ji,jj,jp_pcs0_2d + 12) = znitrpot(ji,jj,1) * 1.e-7 &
-      &                                   * fse3t(ji,jj,1) * 1.e+3 / rfact2 * tmask(ji,jj,1)
-         END DO
-      END DO
+#if defined key_trc_diaadd
+#  if  ! defined key_iomput
+      trc2d(:,:,jp_pcs0_2d + 11) = zirondep(:,:,1)         * 1.e+3 * rfact2r * fse3t(:,:,1) * tmask(:,:,1)
+      trc2d(:,:,jp_pcs0_2d + 12) = znitrpot(:,:,1) * 1.e-7 * 1.e+3 / rfact2  * fse3t(:,:,1) * tmask(:,:,1)
+# else
+      ! write diagnostics
+      zw2d(:,:) = zirondep(:,:,1)         * 1.e+3 * rfact2r * fse3t(:,:,1) * tmask(:,:,1)
+      IF( jnt == nrdttrc ) CALL iom_put( "Fedep", zw2d )
+      zw2d(:,:) = znitrpot(:,:,1) * 1.e-7 * 1.e+3 / rfact2  * fse3t(:,:,1) * tmask(:,:,1)
+      IF( jnt == nrdttrc ) CALL iom_put( "Nfix", zw2d  ) 
+# endif
+
 # endif
       !
Index: trunk/NEMO/TOP_SRC/PISCES/p4zsink.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zsink.F90	(revision 1455)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zsink.F90	(revision 1457)
@@ -14,5 +14,5 @@
    USE sms_pisces
    USE prtctl_trc
-
+   USE iom
 
    IMPLICIT NONE
@@ -34,4 +34,7 @@
    REAL(wp) ::   &
      xstep , xstep2            !: Time step duration for biology
+
+   INTEGER  :: &
+      iksed  = 10              !
 
 #if  defined key_kriest
@@ -89,5 +92,4 @@
       INTEGER, INTENT(in) :: kt, jnt
       INTEGER  :: ji, jj, jk
-      INTEGER  :: iksed
       REAL(wp) :: zagg1, zagg2, zagg3, zagg4, zagg5, zaggsi, zaggsh
       REAL(wp) :: zagg , zaggdoc, znumdoc
@@ -95,6 +97,11 @@
       REAL(wp) :: zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5
       REAL(wp) :: zval1, zval2, zval3, zval4
-#if defined key_trc_dia3d
-      REAL(wp) ::   zrfact2
+#if defined key_trc_diaadd
+      REAL(wp) :: zrfact2
+      INTEGER  :: iksed1
+#if defined key_iomput
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zw2d
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d
+#endif
 #endif
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   znum3d
@@ -113,5 +120,4 @@
 
        znum3d(:,:,:) = 0.e0
-       iksed = 10
        zval1 = 1. + xkr_zeta
        zval2 = 1. + xkr_zeta + xkr_eta
@@ -281,27 +287,43 @@
 #if defined key_trc_diaadd
       zrfact2 = 1.e3 * rfact2r
-      DO jj = 1, jpj 
-         DO ji = 1, jpi
-            trc2d(ji,jj, jp_pcs0_2d + 4) = sinking (ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
-            trc2d(ji,jj, jp_pcs0_2d + 5) = sinking2(ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
-            trc2d(ji,jj, jp_pcs0_2d + 6) = sinkfer (ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
-            trc2d(ji,jj, jp_pcs0_2d + 7) = sinksil (ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
-            trc2d(ji,jj, jp_pcs0_2d + 8) = sinkcal (ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
-         ENDDO
-      ENDDO
-#  if defined key_trc_dia3d
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               trc3d(ji,jj,jk,jp_pcs0_3d + 11) = sinking (ji,jj,jk) * zrfact2 * tmask(ji,jj,jk)
-               trc3d(ji,jj,jk,jp_pcs0_3d + 12) = sinking2(ji,jj,jk) * zrfact2 * tmask(ji,jj,jk)
-               trc3d(ji,jj,jk,jp_pcs0_3d + 13) = sinksil (ji,jj,jk) * zrfact2 * tmask(ji,jj,jk)
-               trc3d(ji,jj,jk,jp_pcs0_3d + 14) = sinkcal (ji,jj,jk) * zrfact2 * tmask(ji,jj,jk)
-               trc3d(ji,jj,jk,jp_pcs0_3d + 15) = znum3d  (ji,jj,jk)           * tmask(ji,jj,jk)
-               trc3d(ji,jj,jk,jp_pcs0_3d + 16) = wsbio3  (ji,jj,jk)           * tmask(ji,jj,jk)
-               trc3d(ji,jj,jk,jp_pcs0_3d + 17) = wsbio4  (ji,jj,jk)           * tmask(ji,jj,jk)
-            ENDDO
-         ENDDO
-      ENDDO
+      iksed1 = iksed + 1
+#  if ! defined key_iomput
+      trc2d(:,:  ,jp_pcs0_2d + 4)  = sinking (:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:  ,jp_pcs0_2d + 5)  = sinking2(:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:  ,jp_pcs0_2d + 6)  = sinkfer (:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:  ,jp_pcs0_2d + 7)  = sinksil (:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:  ,jp_pcs0_2d + 8)  = sinkcal (:,:,iksed1) * 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(:,:,:)
+#else
+      zw2d(:,:  )  = sinking (:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      IF( jnt == nrdttrc ) CALL iom_put( "PMO", zw2d )
+      zw2d(:,:  )  = sinking2(:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      IF( jnt == nrdttrc ) CALL iom_put( "PMO2", zw2d )
+      zw2d(:,:  )  = sinkfer (:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      IF( jnt == nrdttrc ) CALL iom_put( "ExpFe1", zw2d )
+      zw2d(:,:  )  = sinksil (:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      IF( jnt == nrdttrc ) CALL iom_put( "ExpSi", zw2d )
+      zw2d(:,:  )  = sinkcal (:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      IF( jnt == nrdttrc ) CALL iom_put( "ExpCaCO3", zw2d )
+      zw3d(:,:,:)  = sinking (:,:,:)      * zrfact2 * tmask(:,:,:)
+      IF( jnt == nrdttrc ) CALL iom_put( "POCFlx", zw3d )
+      zw3d(:,:,:)  = sinking2(:,:,:)      * zrfact2 * tmask(:,:,:)
+      IF( jnt == nrdttrc ) CALL iom_put( "GOCFlx", zw3d )
+      zw3d(:,:,:)  = sinksil (:,:,:)      * zrfact2 * tmask(:,:,:)
+      IF( jnt == nrdttrc ) CALL iom_put( "SiFlx", zw3d )
+      zw3d(:,:,:)  = sinkcal (:,:,:)      * zrfact2 * tmask(:,:,:)
+      IF( jnt == nrdttrc ) CALL iom_put( "CaCO3Flx", zw3d )
+      zw3d(:,:,:)  = znum3d  (:,:,:)                * tmask(:,:,:)
+      IF( jnt == nrdttrc ) CALL iom_put( "xnum", zw3d )
+      zw3d(:,:,:)  = wsbio3  (:,:,:)                * tmask(:,:,:)
+      IF( jnt == nrdttrc ) CALL iom_put( "W1", zw3d )
+      zw3d(:,:,:)  = wsbio4  (:,:,:)                * tmask(:,:,:)
+      IF( jnt == nrdttrc ) CALL iom_put( "W2", zw3d )
 #  endif
 
@@ -463,5 +485,4 @@
       INTEGER, INTENT(in) :: kt, jnt
       INTEGER  ::   ji, jj, jk
-      INTEGER  ::   iksed
       REAL(wp) ::   zagg1, zagg2, zagg3, zagg4
       REAL(wp) ::   zagg , zaggfe, zaggdoc, zaggdoc2
@@ -469,4 +490,8 @@
 #if defined key_trc_dia3d
       REAL(wp) ::   zrfact2
+      INTEGER  ::   iksed1
+#endif
+#if defined key_iomput
+      REAL(wp), DIMENSION(jpi,jpj) ::   zw2d
 #endif
       CHARACTER (len=25) :: charout
@@ -481,7 +506,4 @@
 !    by data and from the coagulation theory
 !    -----------------------------------------------------------
-
-      iksed = 10
-
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -546,57 +568,38 @@
             DO ji = 1, jpi
                zfact = xstep * xdiss(ji,jj,jk)
-
-!    Part I : Coagulation dependent on turbulence
-!    ----------------------------------------------
-
+               !  Part I : Coagulation dependent on turbulence
 # if defined key_off_degrad
                zagg1 = 940.* zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc) * facvol(ji,jj,jk)
+               zagg2 = 1.054e4 * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc) * facvol(ji,jj,jk)
 # else
                zagg1 = 940.* zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
-# endif
-
-# if defined key_off_degrad
-               zagg2 = 1.054e4 * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc) * facvol(ji,jj,jk)
-# else
                zagg2 = 1.054e4 * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
 # endif
 
-!    Aggregation of small into large particles
-!    Part II : Differential settling
-!    ----------------------------------------------
-
+               ! Part II : Differential settling
+
+               !  Aggregation of small into large particles
 # if defined key_off_degrad
                zagg3 = 0.66 * xstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc) * facvol(ji,jj,jk)
+               zagg4 = 0.e0 * xstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc) * facvol(ji,jj,jk)
 # else
                zagg3 = 0.66 * xstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
-# endif
-
-# if defined key_off_degrad
-               zagg4 = 0.e0 * xstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc) * facvol(ji,jj,jk)
-# else
                zagg4 = 0.e0 * xstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
 # endif
-
                zagg   = zagg1 + zagg2 + zagg3 + zagg4
                zaggfe = zagg * trn(ji,jj,jk,jpsfe) / ( trn(ji,jj,jk,jppoc) + rtrn )
 
-!     Aggregation of DOC to small particles
-!     --------------------------------------
-
+               ! Aggregation of DOC to small particles
+#if defined key_off_degrad
                zaggdoc = ( 80.* trn(ji,jj,jk,jpdoc) + 698. * trn(ji,jj,jk,jppoc) )       &
-# if defined key_off_degrad
-                  &      * facvol(ji,jj,jk)                           &
-# endif
-                  &      * zfact * trn(ji,jj,jk,jpdoc)
-
+                  &      * facvol(ji,jj,jk)  * zfact * trn(ji,jj,jk,jpdoc)
                zaggdoc2 = 1.05e4 * zfact * trn(ji,jj,jk,jpgoc)   &
-# if defined key_off_degrad
-                  &        * facvol(ji,jj,jk)                            &
-# endif      
-                  &        * trn(ji,jj,jk,jpdoc)
-!
-!  Update the trends
-!  -----------------
-!
+                  &      * facvol(ji,jj,jk) * trn(ji,jj,jk,jpdoc)
+#else
+               zaggdoc = ( 80.* trn(ji,jj,jk,jpdoc) + 698. * trn(ji,jj,jk,jppoc) )    &
+                  &      *  zfact * trn(ji,jj,jk,jpdoc)
+               zaggdoc2 = 1.05e4 * zfact * trn(ji,jj,jk,jpgoc) * trn(ji,jj,jk,jpdoc)
+#endif
+               !  Update the trends
                tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zagg + zaggdoc
                tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zagg + zaggdoc2
@@ -604,22 +607,35 @@
                tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zaggfe
                tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc2
-
+               !
             END DO
          END DO
       END DO
 
-# if defined key_trc_diaadd
+#if defined key_trc_diaadd
       zrfact2 = 1.e3 * rfact2r
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            trc2d(ji,jj, jp_pcs0_2d + 4) = sinking (ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
-            trc2d(ji,jj, jp_pcs0_2d + 5) = sinking2(ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
-            trc2d(ji,jj, jp_pcs0_2d + 6) = sinkfer (ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
-            trc2d(ji,jj, jp_pcs0_2d + 7) = sinkfer2(ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
-            trc2d(ji,jj, jp_pcs0_2d + 8) = sinksil (ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
-            trc2d(ji,jj, jp_pcs0_2d + 9) = sinkcal (ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
-         ENDDO
-      ENDDO
-# endif
+      iksed1 = iksed + 1
+#  if ! defined key_iomput
+      trc2d(:,:,jp_pcs0_2d + 4) = sinking (:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:,jp_pcs0_2d + 5) = sinking2(:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:,jp_pcs0_2d + 6) = sinkfer (:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:,jp_pcs0_2d + 7) = sinkfer2(:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:,jp_pcs0_2d + 8) = sinksil (:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      trc2d(:,:,jp_pcs0_2d + 9) = sinkcal (:,:,iksed1) * zrfact2 * tmask(:,:,1)
+#  else
+      ! write diagnostics 
+      zw2d(:,:) = sinking (:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      IF( jnt == nrdttrc ) CALL iom_put( "PMO", zw2d )
+      zw2d(:,:) = sinking2(:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      IF( jnt == nrdttrc ) CALL iom_put( "PMO2", zw2d )
+      zw2d(:,:) = sinkfer (:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      IF( jnt == nrdttrc ) CALL iom_put( "ExpFe1", zw2d )
+      zw2d(:,:) = sinkfer2(:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      IF( jnt == nrdttrc ) CALL iom_put( "ExpFe2", zw2d )
+      zw2d(:,:) = sinksil (:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      IF( jnt == nrdttrc ) CALL iom_put( "ExpSi", zw2d )
+      zw2d(:,:) = sinkcal (:,:,iksed1) * zrfact2 * tmask(:,:,1)
+      IF( jnt == nrdttrc ) CALL iom_put( "ExpCaCO3", zw2d )
+#  endif
+#endif
       !
        IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
Index: trunk/NEMO/TOP_SRC/PISCES/trclsm_pisces.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/trclsm_pisces.F90	(revision 1455)
+++ trunk/NEMO/TOP_SRC/PISCES/trclsm_pisces.F90	(revision 1457)
@@ -47,5 +47,5 @@
       CHARACTER (len=32) ::  clname
       !!
-#if defined key_trc_diaadd
+#if defined key_trc_diaadd && ! defined key_iomput
       INTEGER ::  jl, jn
       ! definition of additional diagnostic as a structure
@@ -64,5 +64,5 @@
       NAMELIST/nampiskrp/ xkr_eta, xkr_zeta, xkr_mass_min, xkr_mass_max
 #endif
-#if defined key_trc_diaadd
+#if defined key_trc_diaadd && ! defined key_iomput
       NAMELIST/nampisdia/ nwritedia, pisdia3d, pisdia2d     ! additional diagnostics
 #endif
@@ -122,5 +122,5 @@
 #endif
       !
-#if defined key_trc_diaadd
+#if defined key_trc_diaadd && ! defined key_iomput
 
       ! Namelist namlobdia
Index: trunk/NEMO/TOP_SRC/TRP/trcstp.F90
===================================================================
--- trunk/NEMO/TOP_SRC/TRP/trcstp.F90	(revision 1455)
+++ 	(revision )
@@ -1,89 +1,0 @@
-MODULE trcstp
-   !!======================================================================
-   !!                       ***  MODULE trcstp  ***
-   !! Time-stepping    : time loop of opa for passive tracer
-   !!======================================================================
-#if defined key_top
-   !!----------------------------------------------------------------------
-   !!   trc_stp      : passive tracer system time-stepping
-   !!----------------------------------------------------------------------
-   !! * Modules used
-   USE oce_trc          ! ocean dynamics and active tracers variables
-   USE trp_trc
-   USE trctrp           ! passive tracers transport
-   USE trcsms           ! passive tracers sources and sinks
-   USE prtctl_trc       ! Print control for debbuging
-   USE trcdia
-   USE trcrst
-   USE trdmld_trc_oce
-   USE trdmld_trc
-
-   IMPLICIT NONE
-   PRIVATE
-
-   !! * Routine accessibility
-   PUBLIC trc_stp           ! called by step
-   !!----------------------------------------------------------------------
-   !!   TOP 1.0 , LOCEAN-IPSL (2005) 
-   !! $Id$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
-   !!----------------------------------------------------------------------
-
-CONTAINS
-
-   SUBROUTINE trc_stp( kt, kindic )
-      !!-------------------------------------------------------------------
-      !!                     ***  ROUTINE trc_stp  ***
-      !!                      
-      !! ** Purpose : Time loop of opa for passive tracer
-      !! 
-      !! ** Method  : 
-      !!              Compute the passive tracers trends 
-      !!              Update the passive tracers
-      !!
-      !! History :
-      !!   9.0  !  04-03  (C. Ethe)  Original
-      !!-------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in ) ::  kt  ! ocean time-step index
-      INTEGER, INTENT( in ) ::  kindic
-      CHARACTER (len=25)    ::  charout
-
-      ! this ROUTINE is called only every ndttrc time step
-      IF( MOD( kt , ndttrc ) /= 0 ) RETURN
-
-      IF(ln_ctl) THEN
-         WRITE(charout,FMT="('kt =', I4,'  d/m/y =',I2,I2,I4)") kt, nday, nmonth, nyear
-         CALL prt_ctl_trc_info(charout)
-      ENDIF
-
-      IF( kt == nittrc000 .AND. lk_trdmld_trc  )  &
-         &                   CALL trd_mld_trc_init        ! trends: Mixed-layer
-
-                             CALL trc_rst_opn( kt )       ! Open tracer restart file 
-
-                             CALL trc_sms( kt )           ! tracers: sink and source
-
-                             CALL trc_trp( kt )           ! transport of passive tracers
-
-      IF( lrst_trc )         CALL trc_rst_wri( kt )       ! write tracer restart file
-
-                             CALL trc_dia( kt, kindic )   ! diagnostics
-
-      IF( lk_trdmld_trc  )   CALL trd_mld_trc( kt )     ! trends: Mixed-layer
-
-   END SUBROUTINE trc_stp
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default key                                     NO passive tracers
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE trc_stp( kt )        ! Empty routine
-      INTEGER, INTENT(in) :: kt
-      WRITE(*,*) 'trc_stp: You should not have seen this print! error?', kt
-   END SUBROUTINE trc_stp
-#endif
-
-   !!======================================================================
-END MODULE trcstp
Index: trunk/NEMO/TOP_SRC/trc.F90
===================================================================
--- trunk/NEMO/TOP_SRC/trc.F90	(revision 1455)
+++ trunk/NEMO/TOP_SRC/trc.F90	(revision 1457)
@@ -64,5 +64,5 @@
    INTEGER , PUBLIC ::   nwritetrc   !: time step frequency for concentration outputs (namelist)
    
-# if defined key_trc_diaadd
+# if defined key_trc_diaadd && ! defined key_iomput
    !! additional 2D/3D outputs namelist
    !! --------------------------------------------------
Index: trunk/NEMO/TOP_SRC/trcdia.F90
===================================================================
--- trunk/NEMO/TOP_SRC/trcdia.F90	(revision 1455)
+++ trunk/NEMO/TOP_SRC/trcdia.F90	(revision 1457)
@@ -11,5 +11,5 @@
    !!                  !  2008-05 (C. Ethe re-organization)
    !!----------------------------------------------------------------------
-#if defined key_top
+#if defined key_top && ! defined key_iomput
    !!----------------------------------------------------------------------
    !!   'key_top'                                                TOP models
@@ -29,5 +29,4 @@
    USE lib_mpp
    USE ioipsl
-   USE iom
 
    IMPLICIT NONE
@@ -70,5 +69,5 @@
 CONTAINS
 
-   SUBROUTINE trc_dia( kt, kindic )  
+   SUBROUTINE trc_dia( kt )  
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE trc_dia  ***
@@ -76,5 +75,6 @@
       !! ** Purpose :   output passive tracers fields 
       !!---------------------------------------------------------------------
-      INTEGER, INTENT( in ) :: kt, kindic
+      INTEGER, INTENT( in ) :: kt
+      INTEGER               :: kindic
       !!---------------------------------------------------------------------
       
@@ -117,6 +117,4 @@
       ! Initialisation
       ! --------------
-
-      CALL iom_setkt( kt + ndttrc - 1 )
 
       ! local variable for debugging
@@ -212,5 +210,4 @@
          cltra = ctrcnm(jn)      ! short title for tracer
          IF( lutsav(jn) ) CALL histwrite( nit5, cltra, it, trn(:,:,:,jn), ndimt50, ndext50 )
-         CALL iom_put( ctrcnm(jn), trn(:,:,:,jn) )
       END DO
 
@@ -219,5 +216,4 @@
       IF( kt == nitend .OR. kindic < 0 )   CALL histclo( nit5 )
       !
-      CALL iom_setkt( kt )
 
    END SUBROUTINE trcdit_wr
@@ -257,5 +253,4 @@
       ! -----------------
       
-      CALL iom_setkt( kt + ndttrc - 1 )
 
       ! local variable for debugging
@@ -445,5 +440,4 @@
             END DO
          END IF
-         CALL iom_put( ctrcnm(jn), trn(:,:,:,jn) )
       END DO
 
@@ -456,5 +450,4 @@
       ENDIF
       !
-      CALL iom_setkt( kt )
 
    END SUBROUTINE trcdid_wr
@@ -500,6 +493,4 @@
       ! Initialisation
       ! --------------
-   
-      CALL iom_setkt( kt + ndttrc - 1 )
       
       ! local variable for debugging
@@ -597,5 +588,4 @@
          cltra = ctrc3d(jl)   ! short title for 3D diagnostic
          CALL histwrite( nitd, cltra, it, trc3d(:,:,:,jl), ndimt50 ,ndext50)
-         CALL iom_put( cltra, trc3d(:,:,:,jl) )
       END DO
 
@@ -604,5 +594,4 @@
          cltra = ctrc2d(jl)   ! short title for 2D diagnostic
          CALL histwrite(nitd, cltra, it, trc2d(:,:,jl), ndimt51  ,ndext51)
-         CALL iom_put( cltra, trc2d(:,:,jl) )
       END DO
 
@@ -611,5 +600,4 @@
       IF( kt == nitend .OR. kindic < 0 )   CALL histclo(nitd)
       !
-      CALL iom_setkt( kt )
 
    END SUBROUTINE trcdii_wr
@@ -657,5 +645,4 @@
       ! --------------
 
-      CALL iom_setkt( kt + ndttrc - 1 )
       
       ! local variable for debugging
@@ -735,5 +722,4 @@
          cltra = ctrbio(jl)  ! short title for biological diagnostic
          CALL histwrite(nitb, cltra, it, trbio(:,:,:,jl), ndimt50,ndext50)
-         CALL iom_put( cltra, trbio(:,:,:,jl) )
       END DO
 
@@ -742,5 +728,4 @@
       IF( kt == nitend .OR. kindic < 0 )   CALL histclo( nitb )
       !
-      CALL iom_setkt( kt )
 
    END SUBROUTINE trcdib_wr
@@ -759,5 +744,6 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE trc_dia                      ! Empty routine   
+   SUBROUTINE trc_dia( kt )                      ! Empty routine   
+      INTEGER, INTENT(in) :: kt
    END SUBROUTINE trc_dia   
 
Index: trunk/NEMO/TOP_SRC/trcini.F90
===================================================================
--- trunk/NEMO/TOP_SRC/trcini.F90	(revision 1455)
+++ trunk/NEMO/TOP_SRC/trcini.F90	(revision 1457)
@@ -27,5 +27,8 @@
    USE trcini_c14b     ! C14 bomb initialisation
    USE trcini_my_trc   ! MY_TRC   initialisation
-   USE trcdta    
+   USE trcdta   
+#if defined key_off_tra 
+   USE daymod
+#endif
    USE zpshde_trc      ! partial step: hor. derivative 
    USE in_out_manager  ! I/O manager
@@ -107,4 +110,7 @@
 
       IF( .NOT. lrsttr ) THEN 
+#if defined key_off_tra
+         CALL day_init      ! calendar
+#endif
 # if defined key_dtatrc
          ! Initialization of tracer from a file that may also be used for damping
@@ -117,4 +123,7 @@
       ELSE
          CALL trc_rst_read      ! restart from a file
+#if defined key_off_tra
+         CALL day_init          ! calendar
+#endif
       ENDIF
 
Index: trunk/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- trunk/NEMO/TOP_SRC/trcrst.F90	(revision 1455)
+++ trunk/NEMO/TOP_SRC/trcrst.F90	(revision 1457)
@@ -154,5 +154,4 @@
                   zco3    = ( zcaralk - trn(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1
                   zbicarb = ( 2. * trn(ji,jj,jk,jpdic) - zcaralk )
-!!                  write(numout,*) 'plante :',ji,jj,jk,ztmas,ztmas1,ak23(ji,jj,jk),zbicarb ,zco3
                   hi(ji,jj,jk) = ( ak23(ji,jj,jk) * zbicarb / zco3 ) * ztmas + 1.e-9 * ztmas1
                END DO
Index: trunk/NEMO/TOP_SRC/trcstp.F90
===================================================================
--- trunk/NEMO/TOP_SRC/trcstp.F90	(revision 1457)
+++ trunk/NEMO/TOP_SRC/trcstp.F90	(revision 1457)
@@ -0,0 +1,90 @@
+MODULE trcstp
+   !!======================================================================
+   !!                       ***  MODULE trcstp  ***
+   !! Time-stepping    : time loop of opa for passive tracer
+   !!======================================================================
+#if defined key_top
+   !!----------------------------------------------------------------------
+   !!   trc_stp      : passive tracer system time-stepping
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE oce_trc          ! ocean dynamics and active tracers variables
+   USE trp_trc
+   USE trctrp           ! passive tracers transport
+   USE trcsms           ! passive tracers sources and sinks
+   USE prtctl_trc       ! Print control for debbuging
+   USE trcdia
+   USE trcwri
+   USE trcrst
+   USE trdmld_trc_oce
+   USE trdmld_trc
+   USE iom
+   USE in_out_manager
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC trc_stp           ! called by step
+   !!----------------------------------------------------------------------
+   !!   TOP 1.0 , LOCEAN-IPSL (2005) 
+   !! $Id: trcstp.F90 1285 2009-02-03 13:38:51Z cetlod $ 
+   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE trc_stp( kt )
+      !!-------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_stp  ***
+      !!                      
+      !! ** Purpose : Time loop of opa for passive tracer
+      !! 
+      !! ** Method  : 
+      !!              Compute the passive tracers trends 
+      !!              Update the passive tracers
+      !!
+      !! History :
+      !!   9.0  !  04-03  (C. Ethe)  Original
+      !!-------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT( in ) ::  kt  ! ocean time-step index
+      CHARACTER (len=25)    ::  charout
+
+      ! this ROUTINE is called only every ndttrc time step
+      IF( MOD( kt , ndttrc ) /= 0 ) RETURN
+
+      IF(ln_ctl) THEN
+         WRITE(charout,FMT="('kt =', I4,'  d/m/y =',I2,I2,I4)") kt, nday, nmonth, nyear
+         CALL prt_ctl_trc_info(charout)
+      ENDIF
+
+      IF( kt == nittrc000 .AND. lk_trdmld_trc  )  &
+         &                   CALL trd_mld_trc_init        ! trends: Mixed-layer
+                             CALL trc_rst_opn( kt )       ! Open tracer restart file 
+                             CALL trc_sms( kt )           ! tracers: sink and source
+                             CALL trc_trp( kt )           ! transport of passive tracers
+      IF( kt == nittrc000 )  CALL iom_close( numrtr )     ! close input  passive tracers restart file
+      IF( lrst_trc )         CALL trc_rst_wri( kt )       ! write tracer restart file
+      IF( lk_iomput ) THEN
+                             CALL trc_wri( kt )           ! output of passive tracers
+      ELSE
+                             CALL trc_dia( kt )   ! diagnostics
+      ENDIF
+      IF( lk_trdmld_trc  )   CALL trd_mld_trc( kt )     ! trends: Mixed-layer
+
+   END SUBROUTINE trc_stp
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default key                                     NO passive tracers
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_stp( kt )        ! Empty routine
+      INTEGER, INTENT(in) :: kt
+      WRITE(*,*) 'trc_stp: You should not have seen this print! error?', kt
+   END SUBROUTINE trc_stp
+#endif
+
+   !!======================================================================
+END MODULE trcstp
Index: trunk/NEMO/TOP_SRC/trcwri.F90
===================================================================
--- trunk/NEMO/TOP_SRC/trcwri.F90	(revision 1457)
+++ trunk/NEMO/TOP_SRC/trcwri.F90	(revision 1457)
@@ -0,0 +1,74 @@
+MODULE trcwri
+   !!======================================================================
+   !!                       *** MODULE trcwri ***
+   !!    TOP :   Output of passive tracers 
+   !!======================================================================
+   !!             1.0  !  
+   !!                  !  2009-05 (C. Ethe )
+   !!----------------------------------------------------------------------
+#if defined key_top &&  defined key_iomput
+   !!----------------------------------------------------------------------
+   !!   'key_top' && 'key_iomput'                              TOP models
+   !!----------------------------------------------------------------------
+   !! trc_wri     :  outputs of concentration fields
+   !!----------------------------------------------------------------------
+   USE trc
+   USE iom
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC trc_wri      
+
+   !! * Substitutions
+#  include "top_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 1.0 , LOCEAN-IPSL (2005) 
+   !! $Id: trcdia.F90 1450 2009-05-15 14:12:12Z cetlod $ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE trc_wri( kt )  
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_wri  ***
+      !!
+      !! ** Purpose :   output passive tracers fields 
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT( in ) :: kt       ! ocean time-step
+      INTEGER               :: jn
+      CHARACTER (len=20)    :: cltra
+      !!---------------------------------------------------------------------
+ 
+      ! Initialisation
+      ! --------------
+
+      CALL iom_setkt( kt + ndttrc - 1 ) ! set the passive tracer time step
+
+      ! write the tracer concentrations in the file
+      ! ---------------------------------------
+      DO jn = 1, jptra
+         cltra = ctrcnm(jn)                   ! short title for tracer
+         CALL iom_put( cltra, trn(:,:,:,jn) )
+      END DO
+      !
+      CALL iom_setkt( kt )       ! set the model time step
+
+      !
+   END SUBROUTINE trc_wri
+
+#else
+   !!----------------------------------------------------------------------
+   !!  Dummy module :                                     No passive tracer
+   !!----------------------------------------------------------------------
+   PUBLIC trc_wri
+CONTAINS
+   SUBROUTINE trc_wri( kt )                     ! Empty routine   
+   INTEGER, INTENT(in) :: kt
+   END SUBROUTINE trc_wri
+
+#endif
+
+   !!======================================================================
+END MODULE trcwri
