Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zbio.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zbio.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zbio.F90	(revision 2528)
@@ -41,7 +41,7 @@
 #  include "top_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -84,11 +84,11 @@
       CALL p4z_sink ( kt, jnt )     ! vertical flux of particulate organic matter
       CALL p4z_opt  ( kt, jnt )     ! Optic: PAR in the water column
-      CALL p4z_lim  ( kt, jnt )     ! co-limitations by the various nutrients
+      CALL p4z_lim  ( kt      )     ! co-limitations by the various nutrients
       CALL p4z_prod ( kt, jnt )     ! phytoplankton growth rate over the global ocean. 
       !                             ! (for each element : C, Si, Fe, Chl )
-      CALL p4z_rem  ( kt, jnt )     ! remineralization terms of organic matter+scavenging of Fe
-      CALL p4z_mort ( kt, jnt )     ! phytoplankton mortality
+      CALL p4z_rem  ( kt      )     ! remineralization terms of organic matter+scavenging of Fe
+      CALL p4z_mort ( kt      )     ! phytoplankton mortality
       !                             ! zooplankton sources/sinks routines 
-      CALL p4z_micro( kt, jnt )           ! microzooplankton
+      CALL p4z_micro( kt      )           ! microzooplankton
       CALL p4z_meso ( kt, jnt )           ! mesozooplankton
 
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zche.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zche.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zche.F90	(revision 2528)
@@ -149,7 +149,7 @@
 #include "top_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -181,8 +181,8 @@
 
             !                             ! SET ABSOLUTE TEMPERATURE
-            ztkel = tn(ji,jj,1) + 273.16
+            ztkel = tsn(ji,jj,1,jp_tem) + 273.16
             zqtt  = ztkel * 0.01
             zqtt2 = zqtt * zqtt
-            zsal  = sn(ji,jj,1) + (1.- tmask(ji,jj,1) ) * 35.
+            zsal  = tsn(ji,jj,1,jp_sal) + (1.- tmask(ji,jj,1) ) * 35.
             zlqtt = LOG( zqtt )
 
@@ -214,7 +214,7 @@
 
                ! SET ABSOLUTE TEMPERATURE
-               ztkel   = tn(ji,jj,jk) + 273.16
+               ztkel   = tsn(ji,jj,jk,jp_tem) + 273.16
                zqtt    = ztkel * 0.01
-               zsal    = sn(ji,jj,jk) + ( 1.-tmask(ji,jj,jk) ) * 35.
+               zsal    = tsn(ji,jj,jk,jp_sal) + ( 1.-tmask(ji,jj,jk) ) * 35.
                zsqrt  = SQRT( zsal )
                zsal15  = zsqrt * zsal
@@ -224,5 +224,5 @@
                zis2   = zis * zis
                zisqrt = SQRT( zis )
-               ztc     = tn(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20.
+               ztc     = tsn(ji,jj,jk,jp_tem) + ( 1.- tmask(ji,jj,jk) ) * 20.
 
                ! CHLORINITY (WOOSTER ET AL., 1969)
@@ -249,6 +249,6 @@
                   &    + ( cb8 + cb9 * zsqrt + cb10 * zsal ) * zlogt + cb11 * zsqrt * ztkel             &
                   &    + LOG(  ( 1.+ zst / zcks + zft / zckf ) / ( 1.+ zst / zcks )  )
-!!gm zsal**2 to be replaced by a *...
-               zck1    = c10 * ztr + c11 + c12 * zlogt + c13 * zsal + c14 * zsal**2
+
+               zck1    = c10 * ztr + c11 + c12 * zlogt + c13 * zsal + c14 * zsal * zsal
                zck2    = c20 * ztr + c21 + c22 * zsal   + c23 * zsal**2
 
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zflx.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zflx.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zflx.F90	(revision 2528)
@@ -28,4 +28,5 @@
 #endif
    USE lib_mpp
+   USE lib_fortran
 
    IMPLICIT NONE
@@ -33,28 +34,21 @@
 
    PUBLIC   p4z_flx  
-
-   REAL(wp) :: &  ! pre-industrial atmospheric [co2] (ppm) 	
-      atcox  = 0.20946 ,    &  !:
-      atcco2 = 278.            !:
-
-   REAL(wp) :: &
-      xconv  = 0.01/3600      !: coefficients for conversion 
-
-   INTEGER  ::  nspyr         !: number of timestep per year
-
-#if defined key_cpl_carbon_cycle
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::  &
-      oce_co2            !: ocean carbon flux
-   REAL(wp) :: &
-      t_atm_co2_flx,  &  !: Total atmospheric carbon flux per year
-      t_oce_co2_flx      !: Total ocean carbon flux per year
-#endif
+   PUBLIC   p4z_flx_init  
+
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::  oce_co2            !: ocean carbon flux 
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::  satmco2            !: atmospheric pco2
+   REAL(wp)                             ::  t_oce_co2_flx      !: Total ocean carbon flux 
+   REAL(wp)                             ::  t_atm_co2_flx      !: global mean of atmospheric pco2
+   REAL(wp)                             ::  area               !: ocean surface
+   REAL(wp)                             ::  atcco2 = 278.      !: pre-industrial atmospheric [co2] (ppm) 	
+   REAL(wp)                             ::  atcox  = 0.20946   !:
+   REAL(wp)                             ::  xconv  = 0.01/3600 !: coefficients for conversion 
 
    !!* Substitution
 #  include "top_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -75,17 +69,18 @@
       REAL(wp) ::   zph, zah2, zbot, zdic, zalk, zsch_o2, zalka, zsch_co2
       REAL(wp), DIMENSION(jpi,jpj) ::   zkgco2, zkgo2, zh2co3
-#if defined key_trc_diaadd && defined key_iomput
-      REAL(wp), DIMENSION(jpi,jpj) ::  zcflx, zoflx, zkg, zdpco2, zdpo2
+#if defined key_diatrc && defined key_iomput
+      REAL(wp), DIMENSION(jpi,jpj) ::  zoflx, zkg, zdpco2, zdpo2
 #endif
       CHARACTER (len=25) :: charout
 
       !!---------------------------------------------------------------------
-
-
-      IF( kt == nittrc000  )   CALL p4z_flx_init      ! Initialization (first time-step only)
 
       ! SURFACE CHEMISTRY (PCO2 AND [H+] IN
       !     SURFACE LAYER); THE RESULT OF THIS CALCULATION
       !     IS USED TO COMPUTE AIR-SEA FLUX OF CO2
+
+#if defined key_cpl_carbon_cycle
+      satmco2(:,:) = atm_co2(:,:)
+#endif
 
       DO jrorr = 1, 10
@@ -128,5 +123,5 @@
 !CDIR NOVERRCHK
          DO ji = 1, jpi
-            ztc  = MIN( 35., tn(ji,jj,1) )
+            ztc  = MIN( 35., tsn(ji,jj,1,jp_tem) )
             ztc2 = ztc * ztc
             ztc3 = ztc * ztc2 
@@ -138,5 +133,5 @@
             ! 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
+# if defined key_degrad
             zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1) * facvol(ji,jj,1)
 #else
@@ -152,14 +147,9 @@
          DO ji = 1, jpi
             ! Compute CO2 flux for the sea and air
-#if ! defined key_cpl_carbon_cycle
-            zfld = atcco2 * tmask(ji,jj,1) * chemc(ji,jj,1) * zkgco2(ji,jj)
+            zfld = satmco2(ji,jj) * tmask(ji,jj,1) * chemc(ji,jj,1) * zkgco2(ji,jj)
             zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj)
-#else
-            zfld = atm_co2(ji,jj) * tmask(ji,jj,1) * chemc(ji,jj,1) * zkgco2(ji,jj)
-            zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj)
-            ! compute flux of carbon
             oce_co2(ji,jj) = ( zfld - zflu ) * rfact &
                &             * e1t(ji,jj) * e2t(ji,jj) * tmask(ji,jj,1) * 1000.
-#endif
+            ! compute the trend
             tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) / fse3t(ji,jj,1)
 
@@ -169,20 +159,18 @@
             tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + ( zfld16 - zflu16 ) / fse3t(ji,jj,1)
 
-#if defined key_trc_diaadd 
+#if defined key_diatrc 
             ! Save diagnostics
 #  if ! defined key_iomput
-            trc2d(ji,jj,jp_pcs0_2d    ) = ( zfld - zflu )     * 1000. * tmask(ji,jj,1)
+            zfact = 1. / ( e1t(ji,jj) * e2t(ji,jj) ) / rfact
+            trc2d(ji,jj,jp_pcs0_2d    ) = oce_co2(ji,jj) * zfact
             trc2d(ji,jj,jp_pcs0_2d + 1) = ( zfld16 - zflu16 ) * 1000. * tmask(ji,jj,1)
             trc2d(ji,jj,jp_pcs0_2d + 2) = zkgco2(ji,jj) * tmask(ji,jj,1)
-            trc2d(ji,jj,jp_pcs0_2d + 3) = ( atcco2 - zh2co3(ji,jj) / ( chemc(ji,jj,1) + rtrn ) ) &
+            trc2d(ji,jj,jp_pcs0_2d + 3) = ( satmco2(ji,jj) - zh2co3(ji,jj) / ( chemc(ji,jj,1) + rtrn ) ) &
                &                            * tmask(ji,jj,1)
 #  else
-            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)
-            zdpco2(ji,jj) = ( atcco2 - zh2co3(ji,jj)      / ( chemc(ji,jj,1) + rtrn ) ) &
-              &             * tmask(ji,jj,1)
-            zdpo2 (ji,jj) = ( atcox  - trn(ji,jj,1,jpoxy) / ( chemc(ji,jj,2) + rtrn ) ) &
-              &             * tmask(ji,jj,1)
+            zdpco2(ji,jj) = ( satmco2(ji,jj) - zh2co3(ji,jj) / ( chemc(ji,jj,1) + rtrn ) ) * tmask(ji,jj,1)
+            zdpo2 (ji,jj) = ( atcox  - trn(ji,jj,1,jpoxy) / ( chemc(ji,jj,2) + rtrn ) ) * tmask(ji,jj,1)
 #  endif
 #endif
@@ -190,35 +178,21 @@
       END DO
 
-#if defined key_cpl_carbon_cycle
-      ! Total Flux of Carbon
-      DO jj = 1, jpj 
-        DO ji = 1, jpi
-           t_atm_co2_flx = t_atm_co2_flx + atm_co2(ji,jj) * tmask_i(ji,jj)
-           t_oce_co2_flx = t_oce_co2_flx + oce_co2(ji,jj) * tmask_i(ji,jj)
-        END DO
-      END DO
-
-      IF( MOD( kt, nspyr ) == 0 ) THEN
-        IF( lk_mpp ) THEN
-          CALL mpp_sum( t_atm_co2_flx )   ! sum over the global domain
-          CALL mpp_sum( t_oce_co2_flx )   ! sum over the global domain
-        ENDIF
-        ! Conversion in GtC/yr ; negative for outgoing from ocean
-        t_oce_co2_flx = (-1.) * t_oce_co2_flx  * 12. / 1.e15
-        !
-        WRITE(numout,*) ' Atmospheric pCO2    :'
-        WRITE(numout,*) '-------------------- : ',kt,'  ',t_atm_co2_flx
-        WRITE(numout,*) '(ppm)'
-        WRITE(numout,*) 'Total Flux of Carbon out of the ocean :'
-        WRITE(numout,*) '-------------------- : ',t_oce_co2_flx
-        WRITE(numout,*) '(GtC/yr)'
-        t_atm_co2_flx = 0.
-        t_oce_co2_flx = 0.
-# if defined key_iomput
-        CALL iom_put( "tatpco2" , t_atm_co2_flx  )
-        CALL iom_put( "tco2flx" , t_oce_co2_flx  )
-#endif
+      t_oce_co2_flx = t_oce_co2_flx + glob_sum( oce_co2(:,:) )                     ! Cumulative Total Flux of Carbon
+      IF( kt == nitend ) THEN
+         t_atm_co2_flx = glob_sum( satmco2(:,:) * e1t(:,:) * e2t(:,:) )            ! Total atmospheric pCO2
+         !
+         t_oce_co2_flx = (-1.) * t_oce_co2_flx  * 12. / 1.e15                      ! Conversion in PgC ; negative for out of the ocean
+         t_atm_co2_flx = t_atm_co2_flx  / area                                     ! global mean of atmospheric pCO2
+         !
+         IF( lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) ' Global mean of atmospheric pCO2 (ppm) at it= ', kt, ' date= ', ndastp
+            WRITE(numout,*) '------------------------------------------------------- :  ',t_atm_co2_flx
+            WRITE(numout,*)
+            WRITE(numout,*) ' Cumulative total Flux of Carbon out of the ocean (PgC) :'
+            WRITE(numout,*) '-------------------------------------------------------  ',t_oce_co2_flx
+         ENDIF
+         !
       ENDIF
-#endif
 
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
@@ -228,6 +202,6 @@
       ENDIF
 
-# if defined key_trc_diaadd && defined key_iomput
-      CALL iom_put( "Cflx" , zcflx  )
+# if defined key_diatrc && defined key_iomput
+      CALL iom_put( "Cflx" , oce_co2(:,:) / ( e1t(:,:) * e2t(:,:) ) / rfact  )
       CALL iom_put( "Oflx" , zoflx  )
       CALL iom_put( "Kg"   , zkg    )
@@ -246,5 +220,5 @@
       !!
       !! ** Method  :   Read the nampisext namelist and check the parameters
-      !!      called at the first timestep (nittrc000)
+      !!      called at the first timestep (nit000)
       !! ** input   :   Namelist nampisext
       !!
@@ -263,13 +237,13 @@
       ENDIF
 
-      ! number of time step per year  
-      nspyr = INT( nyear_len(1) * rday / rdt )
-
-#if defined key_cpl_carbon_cycle
+      ! interior global domain surface
+      area = glob_sum( e1t(:,:) * e2t(:,:) )  
+
       ! Initialization of Flux of Carbon
-      oce_co2(:,:) = 0.
-      t_atm_co2_flx = 0.
-      t_oce_co2_flx = 0.
-#endif
+      oce_co2(:,:)  = 0._wp
+      t_atm_co2_flx = 0._wp
+      ! Initialisation of atmospheric pco2
+      satmco2(:,:)  = atcco2
+      t_oce_co2_flx = 0._wp
 
    END SUBROUTINE p4z_flx_init
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zint.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zint.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zint.F90	(revision 2528)
@@ -32,7 +32,7 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -55,6 +55,6 @@
       ! -------------------------------------------
 
-      tgfunc (:,:,:) = EXP( 0.063913 * tn(:,:,:) )
-      tgfunc2(:,:,:) = EXP( 0.07608  * tn(:,:,:) )
+      tgfunc (:,:,:) = EXP( 0.063913 * tsn(:,:,:,jp_tem) )
+      tgfunc2(:,:,:) = EXP( 0.07608  * tsn(:,:,:,jp_tem) )
 
       ! Computation of the silicon dependant half saturation
@@ -69,5 +69,5 @@
       END DO
 
-      IF( nday_year == 365 ) THEN
+      IF( nday_year == nyear_len(1) ) THEN
          xksi    = xksimax
          xksimax = 0.e0
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zlim.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zlim.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zlim.F90	(revision 2528)
@@ -23,4 +23,5 @@
 
    PUBLIC p4z_lim    
+   PUBLIC p4z_lim_init    
 
    !! * Shared module variables
@@ -43,12 +44,12 @@
 #  include "top_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
 CONTAINS
 
-   SUBROUTINE p4z_lim( kt, jnt )
+   SUBROUTINE p4z_lim( kt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_lim  ***
@@ -59,5 +60,5 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
+      INTEGER, INTENT(in)  :: kt
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   zlim1, zlim2, zlim3, zlim4, zno3, zferlim
@@ -67,10 +68,7 @@
 
 
-      IF( ( kt * jnt ) == nittrc000  )   CALL p4z_lim_init      ! Initialization (first time-step only)
-
-
-!  Tuning of the iron concentration to a minimum
-!  level that is set to the detection limit
-!  -------------------------------------
+      !  Tuning of the iron concentration to a minimum
+      !  level that is set to the detection limit
+      !  -------------------------------------
 
       DO jk = 1, jpkm1
@@ -85,8 +83,7 @@
       END DO
 
-!  Computation of a variable Ks for iron on diatoms
-!  taking into account that increasing biomass is
-!  made of generally bigger cells
-!  ------------------------------------------------
+      !  Computation of a variable Ks for iron on diatoms taking into account
+      !  that increasing biomass is made of generally bigger cells
+      !  ------------------------------------------------
 
       DO jk = 1, jpkm1
@@ -107,11 +104,9 @@
       END DO
 
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-    
-!      Michaelis-Menten Limitation term for nutrients
-!      Small flagellates
-!      -----------------------------------------------
+     !  Michaelis-Menten Limitation term for nutrients Small flagellates
+     !      -----------------------------------------------
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
               zdenom = 1. / &
                   & ( conc0 * concnnh4 + concnnh4 * trn(ji,jj,jk,jpno3) + conc0 * trn(ji,jj,jk,jpnh4) )
@@ -132,10 +127,9 @@
       END DO
 
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-
-!   Michaelis-Menten Limitation term for nutrients Diatoms
-!   ----------------------------------------------
+      !   Michaelis-Menten Limitation term for nutrients Diatoms
+      !   ----------------------------------------------
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
               zdenom = 1. / &
                   & ( conc1 * concdnh4 + concdnh4 * trn(ji,jj,jk,jpno3) + conc1 * trn(ji,jj,jk,jpnh4) )
@@ -161,5 +155,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               ztemp = MAX( 0., tn(ji,jj,jk) )
+               ztemp = MAX( 0., tsn(ji,jj,jk,jp_tem) )
                xfracal(ji,jj,jk) = caco3r * xlimphy(ji,jj,jk)   &
                   &                       * MAX( 0.0001, ztemp / ( 2.+ ztemp ) )   &
@@ -181,5 +175,5 @@
       !!
       !! ** Method  :   Read the nampislim namelist and check the parameters
-      !!      called at the first timestep (nittrc000)
+      !!      called at the first timestep (nit000)
       !!
       !! ** input   :   Namelist nampislim
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zlys.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zlys.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zlys.F90	(revision 2528)
@@ -27,5 +27,6 @@
    PRIVATE
 
-   PUBLIC   p4z_lys    ! called in p4zprg.F90
+   PUBLIC   p4z_lys         ! called in trcsms_pisces.F90
+   PUBLIC   p4z_lys_init    ! called in trcsms_pisces.F90
 
    !! * Shared module variables
@@ -42,7 +43,7 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -65,5 +66,5 @@
       REAL(wp) ::   zomegaca, zexcess, zexcess0
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zco3
-#if defined key_trc_dia3d && defined key_iomput
+#if defined key_diatrc && defined key_iomput
       REAL(wp) ::   zrfact2
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zcaldiss
@@ -72,9 +73,7 @@
       !!---------------------------------------------------------------------
 
-      IF( kt == nittrc000  )   CALL p4z_lys_init      ! Initialization (first time-step only)
-
       zco3(:,:,:) = 0.
 
-# if defined key_trc_dia3d && defined key_iomput
+# if defined key_diatrc && defined key_iomput
       zcaldiss(:,:,:) = 0.
 # endif
@@ -146,5 +145,5 @@
                !       (ACCORDING TO THIS FORMULATION ALSO SOME PARTICULATE
                !       CACO3 GETS DISSOLVED EVEN IN THE CASE OF OVERSATURATION)
-# if defined key_off_degrad
+# if defined key_degrad
               zdispot = kdca * zexcess * trn(ji,jj,jk,jpcal) * facvol(ji,jj,jk)
 # else
@@ -160,5 +159,5 @@
               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) +      zremco3
 
-# if defined key_trc_dia3d && defined key_iomput
+# if defined key_diatrc && defined key_iomput
               zcaldiss(ji,jj,jk) = zremco3  ! calcite dissolution
 # endif
@@ -167,5 +166,5 @@
       END DO
 
-# if defined key_trc_diaadd &&  defined key_trc_dia3d
+# if defined key_diatrc
 #  if ! defined key_iomput
       trc3d(:,:,:,jp_pcs0_3d    ) = hi  (:,:,:)          * tmask(:,:,:)
@@ -197,5 +196,5 @@
       !!
       !! ** Method  :   Read the nampiscal namelist and check the parameters
-      !!      called at the first timestep (nittrc000)
+      !!      called at the first timestep (nit000)
       !!
       !! ** input   :   Namelist nampiscal
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zmeso.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zmeso.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zmeso.F90	(revision 2528)
@@ -26,5 +26,6 @@
    PRIVATE
 
-   PUBLIC   p4z_meso         ! called in p4zbio.F90
+   PUBLIC   p4z_meso              ! called in p4zbio.F90
+   PUBLIC   p4z_meso_init         ! called in trcsms_pisces.F90
 
    !! * Shared module variables
@@ -47,12 +48,12 @@
 #  include "top_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
 CONTAINS
 
-   SUBROUTINE p4z_meso( kt,jnt )
+   SUBROUTINE p4z_meso( kt, jnt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_meso  ***
@@ -65,35 +66,18 @@
       INTEGER  :: ji, jj, jk
       REAL(wp) :: zcompadi, zcompaph, zcompapoc, zcompaz
-      REAL(wp) :: zfact, zstep, zcompam, zdenom, zgraze2
+      REAL(wp) :: zfact, zcompam, zdenom, zgraze2, zstep
       REAL(wp) :: zgrarem2, zgrafer2, zgrapoc2, zprcaca, zmortz2
 #if defined key_kriest
       REAL znumpoc
 #endif
-      REAL(wp),DIMENSION(jpi,jpj,jpk) :: zrespz2,ztortz2,zgrazd,zgrazz,zgrazpof
-      REAL(wp),DIMENSION(jpi,jpj,jpk) :: zgrazn,zgrazpoc,zgraznf,zgrazf
-      REAL(wp),DIMENSION(jpi,jpj,jpk) :: zgrazfff,zgrazffe
+      REAL(wp) :: zrespz2,ztortz2,zgrazd,zgrazz,zgrazpof
+      REAL(wp) :: zgrazn,zgrazpoc,zgraznf,zgrazf
+      REAL(wp) :: zgrazfff,zgrazffe
       CHARACTER (len=25) :: charout
-#if defined key_trc_diaadd && defined key_trc_dia3d && defined key_iomput
+#if defined key_diatrc && defined key_iomput
       REAL(wp) :: zrfact2
 #endif
 
       !!---------------------------------------------------------------------
-
-
-      IF( ( kt * jnt ) == nittrc000  )   CALL p4z_meso_init      ! Initialization (first time-step only)
-
-      zrespz2 (:,:,:) = 0.
-      ztortz2 (:,:,:) = 0.
-      zgrazd  (:,:,:) = 0.
-      zgrazz  (:,:,:) = 0.
-      zgrazpof(:,:,:) = 0.
-      zgrazn  (:,:,:) = 0.
-      zgrazpoc(:,:,:) = 0.
-      zgraznf (:,:,:) = 0.
-      zgrazf  (:,:,:) = 0.
-      zgrazfff(:,:,:) = 0.
-      zgrazffe(:,:,:) = 0.
-
-      zstep = rfact2 / rday      ! Time step duration for biology
 
       DO jk = 1, jpkm1
@@ -102,29 +86,22 @@
 
                zcompam = MAX( ( trn(ji,jj,jk,jpmes) - 1.e-9 ), 0.e0 )
-# if defined key_off_degrad
-               zfact   = zstep * tgfunc(ji,jj,jk) * zcompam * facvol(ji,jj,jk)
+# if defined key_degrad
+               zstep   = xstep * facvol(ji,jj,jk)
 # else
+               zstep   = xstep
+# endif
                zfact   = zstep * tgfunc(ji,jj,jk) * zcompam
-# endif
-
-!     Respiration rates of both zooplankton
-!     -------------------------------------
-               zrespz2(ji,jj,jk)  = resrat2 * zfact * ( 1. + 3. * nitrfac(ji,jj,jk) )        &
+
+               !  Respiration rates of both zooplankton
+               !  -------------------------------------
+               zrespz2  = resrat2 * zfact * ( 1. + 3. * nitrfac(ji,jj,jk) )        &
                   &     * trn(ji,jj,jk,jpmes) / ( xkmort + trn(ji,jj,jk,jpmes) )
 
-!     Zooplankton mortality. A square function has been selected with
-!     no real reason except that it seems to be more stable and may
-!     mimic predation.
-!     ---------------------------------------------------------------
-               ztortz2(ji,jj,jk) = mzrat2 * 1.e6 * zfact * trn(ji,jj,jk,jpmes)
+               !  Zooplankton mortality. A square function has been selected with
+               !  no real reason except that it seems to be more stable and may mimic predation
+               !  ---------------------------------------------------------------
+               ztortz2 = mzrat2 * 1.e6 * zfact * trn(ji,jj,jk,jpmes)
                !
-            END DO
-         END DO
-      END DO
-
-
-      DO jk = 1,jpkm1
-         DO jj = 1,jpj
-            DO ji = 1,jpi
+
                zcompadi  = MAX( ( trn(ji,jj,jk,jpdia) - 1.e-8 ), 0.e0 )
                zcompaz   = MAX( ( trn(ji,jj,jk,jpzoo) - 1.e-8 ), 0.e0 )
@@ -132,6 +109,6 @@
                zcompapoc = MAX( ( trn(ji,jj,jk,jppoc) - 1.e-8 ), 0.e0 )
 
-!     Microzooplankton grazing
-!     ------------------------
+               !  Microzooplankton grazing
+               !     ------------------------
                zdenom = 1. / (  xkgraz2 + xprefc   * trn(ji,jj,jk,jpdia)   &
                   &                     + xprefz   * trn(ji,jj,jk,jpzoo)   &
@@ -139,107 +116,65 @@
                   &                     + xprefpoc * trn(ji,jj,jk,jppoc)  )
 
-               zgraze2 = grazrat2 * zstep * Tgfunc2(ji,jj,jk) * zdenom    &
-# if defined key_off_degrad
-                  &     * facvol(ji,jj,jk)          &
+               zgraze2 = grazrat2 * zstep * Tgfunc2(ji,jj,jk) * zdenom * trn(ji,jj,jk,jpmes) 
+
+               zgrazd   = zgraze2  * xprefc   * zcompadi
+               zgrazz   = zgraze2  * xprefz   * zcompaz
+               zgrazn   = zgraze2  * xprefp   * zcompaph
+               zgrazpoc = zgraze2  * xprefpoc * zcompapoc
+
+               zgraznf  = zgrazn   * trn(ji,jj,jk,jpnfe) / (trn(ji,jj,jk,jpphy) + rtrn)
+               zgrazf   = zgrazd   * trn(ji,jj,jk,jpdfe) / (trn(ji,jj,jk,jpdia) + rtrn)
+               zgrazpof = zgrazpoc * trn(ji,jj,jk,jpsfe) / (trn(ji,jj,jk,jppoc) + rtrn)
+               
+               !  Mesozooplankton flux feeding on GOC
+               !  ----------------------------------
+# if ! defined key_kriest
+               zgrazffe = grazflux * zstep * wsbio4(ji,jj,jk)          &
+                  &                 * tgfunc2(ji,jj,jk) * trn(ji,jj,jk,jpgoc) * trn(ji,jj,jk,jpmes)
+               zgrazfff = zgrazffe * trn(ji,jj,jk,jpbfe) / (trn(ji,jj,jk,jpgoc) + rtrn)
+# else
+               !!--------------------------- KRIEST3 -------------------------------------------
+               !!               zgrazffe = 0.5 * 1.3e-2 / 5.5e-7 * 0.3 * zstep * wsbio3(ji,jj,jk)     &
+               !!                  &     * tgfunc(ji,jj,jk) * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpmes)    &
+               !! #  if defined key_degrad
+               !!                  &     * facvol(ji,jj,jk)          &
+               !! #  endif
+               !!                  &     /  (trn(ji,jj,jk,jppoc) * 1.e7 + 0.1)
+               !!--------------------------- KRIEST3 -------------------------------------------
+
+              zgrazffe = grazflux * zstep * wsbio3(ji,jj,jk)     &
+                  &                * tgfunc2(ji,jj,jk) * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpmes)
+              zgrazfff = zgrazffe * trn(ji,jj,jk,jpsfe) / (trn(ji,jj,jk,jppoc) + rtrn)
 # endif
-                  &     * trn(ji,jj,jk,jpmes)
-
-               zgrazd(ji,jj,jk)   = zgraze2 * xprefc   * zcompadi
-               zgrazz(ji,jj,jk)   = zgraze2 * xprefz   * zcompaz
-               zgrazn(ji,jj,jk)   = zgraze2 * xprefp   * zcompaph
-               zgrazpoc(ji,jj,jk) = zgraze2 * xprefpoc * zcompapoc
-
-               zgraznf(ji,jj,jk)  = zgrazn(ji,jj,jk)   * trn(ji,jj,jk,jpnfe) &
-                  &                                     / (trn(ji,jj,jk,jpphy) + rtrn)
-               zgrazf(ji,jj,jk)   = zgrazd(ji,jj,jk)   * trn(ji,jj,jk,jpdfe) &
-                  &                                    / (trn(ji,jj,jk,jpdia) + rtrn)
-               zgrazpof(ji,jj,jk) = zgrazpoc(ji,jj,jk) * trn(ji,jj,jk,jpsfe) &
-                  &                                   / (trn(ji,jj,jk,jppoc) + rtrn)
-            END DO
-         END DO
-      END DO
       
-      
-      DO jk = 1,jpkm1
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-               
-!    Mesozooplankton flux feeding on GOC
-!    ----------------------------------
-# if ! defined key_kriest
-#   if ! defined key_off_degrad
-               zgrazffe(ji,jj,jk) = grazflux * zstep * wsbio4(ji,jj,jk)          &
-                  &                 * tgfunc2(ji,jj,jk) * trn(ji,jj,jk,jpgoc) * trn(ji,jj,jk,jpmes)
-#   else
-               zgrazffe(ji,jj,jk) = grazflux * zstep * wsbio4(ji,jj,jk) * facvol(ji,jj,jk)         &
-                  &                 * tgfunc2(ji,jj,jk) * trn(ji,jj,jk,jpgoc) * trn(ji,jj,jk,jpmes)
-#  endif
-               zgrazfff(ji,jj,jk) = zgrazffe(ji,jj,jk)       &
-                  &                 * trn(ji,jj,jk,jpbfe) / (trn(ji,jj,jk,jpgoc) + rtrn)
-# else
-!!--------------------------- KRIEST3 -------------------------------------------
-!!               zgrazffe(ji,jj,jk) = 0.5 * 1.3e-2 / 5.5e-7 * 0.3 * zstep * wsbio3(ji,jj,jk)     &
-!!                  &     * tgfunc(ji,jj,jk) * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpmes)    &
-#  if defined key_off_degrad
-!!                  &     * facvol(ji,jj,jk)          &
-#  endif
-!!                  &     /  (trn(ji,jj,jk,jppoc) * 1.e7 + 0.1)
-!!--------------------------- KRIEST3 -------------------------------------------
-
-#  if ! defined key_off_degrad
-              zgrazffe(ji,jj,jk) = grazflux * zstep * wsbio3(ji,jj,jk)     &
-                  &                * tgfunc2(ji,jj,jk) * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpmes)
-#  else
-              zgrazffe(ji,jj,jk) = grazflux * zstep * wsbio3(ji,jj,jk) * facvol(ji,jj,jk)    &
-                  &               * tgfunc2(ji,jj,jk) * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpmes)
-#  endif
-
-               zgrazfff(ji,jj,jk) = zgrazffe(ji,jj,jk)      &
-                  &                * trn(ji,jj,jk,jpsfe) / (trn(ji,jj,jk,jppoc) + rtrn)
-# endif
-            END DO
-         END DO
-      END DO
-      
-#if defined key_trc_dia3d
-      ! Total grazing ( grazing by microzoo is already computed in p4zmicro ) 
-      grazing(:,:,:) = grazing(:,:,:) + (  zgrazd  (:,:,:) + zgrazz  (:,:,:) + zgrazn(:,:,:) &
-                     &                   + zgrazpoc(:,:,:) + zgrazffe(:,:,:)  )
-#endif
-
-
-      DO jk = 1,jpkm1
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-
-!    Mesozooplankton efficiency
-!    --------------------------
-               zgrarem2 = ( zgrazd(ji,jj,jk) + zgrazz(ji,jj,jk) + zgrazn(ji,jj,jk) &
-                  &     + zgrazpoc(ji,jj,jk) + zgrazffe(ji,jj,jk) )   &
-                  &     * ( 1. - epsher2 - unass2 )
+#if defined key_diatrc
+              ! Total grazing ( grazing by microzoo is already computed in p4zmicro ) 
+              grazing(ji,jj,jk) = grazing(ji,jj,jk) + (  zgrazd + zgrazz + zgrazn + zgrazpoc + zgrazffe )
+#endif
+
+              !    Mesozooplankton efficiency
+              !    --------------------------
+              zgrarem2 = ( zgrazd + zgrazz + zgrazn + zgrazpoc + zgrazffe ) * ( 1. - epsher2 - unass2 )
 #if ! defined key_kriest
-               zgrafer2 = (zgrazf(ji,jj,jk) + zgraznf(ji,jj,jk) + zgrazz(ji,jj,jk) &
-                  &     * ferat3 + zgrazpof(ji,jj,jk) + zgrazfff (ji,jj,jk))*(1.-epsher2-unass2) &
-                  &     + epsher2 * ( &
-                  &      zgrazd(ji,jj,jk)   * MAX((trn(ji,jj,jk,jpdfe) / (trn(ji,jj,jk,jpdia) + rtrn)-ferat3),0.) &
-                  &     + zgrazn(ji,jj,jk)   * MAX((trn(ji,jj,jk,jpnfe) / (trn(ji,jj,jk,jpphy) + rtrn)-ferat3),0.) &
-                  &    + zgrazpoc(ji,jj,jk) * MAX((trn(ji,jj,jk,jpsfe) / (trn(ji,jj,jk,jppoc) + rtrn)-ferat3),0.) &
-                  &    + zgrazffe(ji,jj,jk) * MAX((trn(ji,jj,jk,jpbfe) / (trn(ji,jj,jk,jpgoc) + rtrn)-ferat3),0.)  )
+              zgrafer2 = ( zgrazf + zgraznf + zgrazz * ferat3 + zgrazpof + zgrazfff ) * ( 1.- epsher2 - unass2 ) & 
+                  &     + epsher2 * ( zgrazd   * MAX((trn(ji,jj,jk,jpdfe) / (trn(ji,jj,jk,jpdia) + rtrn)-ferat3),0.) &
+                  &                 + zgrazn   * MAX((trn(ji,jj,jk,jpnfe) / (trn(ji,jj,jk,jpphy) + rtrn)-ferat3),0.) &
+                  &                 + zgrazpoc * MAX((trn(ji,jj,jk,jpsfe) / (trn(ji,jj,jk,jppoc) + rtrn)-ferat3),0.) &
+                  &                 + zgrazffe * MAX((trn(ji,jj,jk,jpbfe) / (trn(ji,jj,jk,jpgoc) + rtrn)-ferat3),0.)  )
 #else
-               zgrafer2 = (zgrazf(ji,jj,jk) + zgraznf(ji,jj,jk) + zgrazz(ji,jj,jk) &
-                  &    * ferat3 + zgrazpof(ji,jj,jk) + zgrazfff(ji,jj,jk) )*(1.-epsher2-unass2) &
-                  &    + epsher2 * ( &
-                  &    zgrazd(ji,jj,jk)   * MAX((trn(ji,jj,jk,jpdfe) / (trn(ji,jj,jk,jpdia) + rtrn)-ferat3),0.) &
-                  &    + zgrazn(ji,jj,jk)   * MAX((trn(ji,jj,jk,jpnfe) / (trn(ji,jj,jk,jpphy) + rtrn)-ferat3),0.) &
-                  &    + zgrazpoc(ji,jj,jk) * MAX((trn(ji,jj,jk,jpsfe) / (trn(ji,jj,jk,jppoc) + rtrn)-ferat3),0.) &
-                  &    + zgrazffe(ji,jj,jk) * MAX((trn(ji,jj,jk,jpsfe) / (trn(ji,jj,jk,jppoc) + rtrn)-ferat3),0.)  )
-
-#endif
-               zgrapoc2 = (zgrazd(ji,jj,jk) + zgrazz(ji,jj,jk)  + zgrazn(ji,jj,jk) &
-                  &    + zgrazpoc(ji,jj,jk) + zgrazffe(ji,jj,jk)) * unass2
+              zgrafer2 = ( zgrazf + zgraznf + zgrazz * ferat3 + zgrazpof + zgrazfff ) * ( 1. - epsher2 - unass2 ) &
+                  &    + epsher2 * ( zgrazd   * MAX((trn(ji,jj,jk,jpdfe) / (trn(ji,jj,jk,jpdia) + rtrn)-ferat3),0.) &
+                  &                + zgrazn   * MAX((trn(ji,jj,jk,jpnfe) / (trn(ji,jj,jk,jpphy) + rtrn)-ferat3),0.) &
+                  &                + zgrazpoc * MAX((trn(ji,jj,jk,jpsfe) / (trn(ji,jj,jk,jppoc) + rtrn)-ferat3),0.) &
+                  &                + zgrazffe * MAX((trn(ji,jj,jk,jpsfe) / (trn(ji,jj,jk,jppoc) + rtrn)-ferat3),0.)  )
+
+#endif
+               !   Update the arrays TRA which contain the biological sources and sinks
+
+               zgrapoc2 =  zgrazd + zgrazz  + zgrazn + zgrazpoc + zgrazffe
 
                tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zgrarem2 * sigma2
                tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zgrarem2 * sigma2
-               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zgrarem2 * (1.-sigma2)
+               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zgrarem2 * ( 1. - sigma2 )
                tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - o2ut * zgrarem2 * sigma2
                tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zgrafer2
@@ -247,39 +182,23 @@
                
 #if defined key_kriest
-               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zgrapoc2
-               tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zgrapoc2 * xkr_dmeso
+               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zgrapoc2 * unass2
+               tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zgrapoc2 * unass2 * xkr_dmeso
 #else
-               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zgrapoc2
-#endif
-            END DO
-         END DO
-      END DO
-
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !
-               !   Update the arrays TRA which contain the biological sources and sinks
-               !   --------------------------------------------------------------------
-               zmortz2 = ztortz2(ji,jj,jk) + zrespz2(ji,jj,jk)
-               tra(ji,jj,jk,jpmes) = tra(ji,jj,jk,jpmes) - zmortz2  &
-                  &    + epsher2 * ( zgrazd(ji,jj,jk) + zgrazz(ji,jj,jk) + zgrazn(ji,jj,jk) &
-                  &    + zgrazpoc(ji,jj,jk) + zgrazffe(ji,jj,jk) )
-               tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) - zgrazd(ji,jj,jk)
-               tra(ji,jj,jk,jpzoo) = tra(ji,jj,jk,jpzoo) - zgrazz(ji,jj,jk)
-               tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) - zgrazn(ji,jj,jk)
-               tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) - zgrazn(ji,jj,jk) * trn(ji,jj,jk,jpnch)  &
-                  &    / ( trn(ji,jj,jk,jpphy) + rtrn )
-               tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) - zgrazd(ji,jj,jk) * trn(ji,jj,jk,jpdch) &
-                  &    / ( trn(ji,jj,jk,jpdia) + rtrn )
-               tra(ji,jj,jk,jpbsi) = tra(ji,jj,jk,jpbsi) - zgrazd(ji,jj,jk) * trn(ji,jj,jk,jpbsi) &
-                  &    / ( trn(ji,jj,jk,jpdia) + rtrn )
-               tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) +  zgrazd(ji,jj,jk) * trn(ji,jj,jk,jpbsi) &
-                  &    / ( trn(ji,jj,jk,jpdia) + rtrn )
-               tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) -  zgraznf(ji,jj,jk)
-               tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) -  zgrazf(ji,jj,jk)
-
-               zprcaca = xfracal(ji,jj,jk) * unass2 * zgrazn(ji,jj,jk)
-#if defined key_trc_dia3d
+               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zgrapoc2 * unass2
+#endif
+               zmortz2 = ztortz2 + zrespz2
+               tra(ji,jj,jk,jpmes) = tra(ji,jj,jk,jpmes) - zmortz2 + epsher2 * zgrapoc2
+               tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) - zgrazd
+               tra(ji,jj,jk,jpzoo) = tra(ji,jj,jk,jpzoo) - zgrazz
+               tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) - zgrazn
+               tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) - zgrazn * trn(ji,jj,jk,jpnch) / ( trn(ji,jj,jk,jpphy) + rtrn )
+               tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) - zgrazd * trn(ji,jj,jk,jpdch) / ( trn(ji,jj,jk,jpdia) + rtrn )
+               tra(ji,jj,jk,jpbsi) = tra(ji,jj,jk,jpbsi) - zgrazd * trn(ji,jj,jk,jpbsi) / ( trn(ji,jj,jk,jpdia) + rtrn )
+               tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) + zgrazd * trn(ji,jj,jk,jpbsi) / ( trn(ji,jj,jk,jpdia) + rtrn )
+               tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) - zgraznf
+               tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) - zgrazf
+
+               zprcaca = xfracal(ji,jj,jk) * unass2 * zgrazn
+#if defined key_diatrc
                prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca  ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)
 #endif
@@ -290,22 +209,15 @@
 #if defined key_kriest
                znumpoc = trn(ji,jj,jk,jpnum) / ( trn(ji,jj,jk,jppoc) + rtrn )
-               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zmortz2  &
-                  &    - zgrazpoc(ji,jj,jk) - zgrazffe(ji,jj,jk)    
-               tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) - zgrazpoc(ji,jj,jk) * znumpoc &
-                  &    + zmortz2  * xkr_dmeso &
-                  &    - zgrazffe(ji,jj,jk)   * znumpoc * wsbio4(ji,jj,jk) &
-                  &    / ( wsbio3(ji,jj,jk) + rtrn )
+               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zmortz2 - zgrazpoc - zgrazffe
+               tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) - zgrazpoc * znumpoc &
+                  &    + zmortz2  * xkr_dmeso - zgrazffe * znumpoc * wsbio4(ji,jj,jk) / ( wsbio3(ji,jj,jk) + rtrn )
                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ferat3 * zmortz2 &
-               &       + unass2 * ( ferat3 * zgrazz(ji,jj,jk) + zgraznf(ji,jj,jk) &
-               &       + zgrazf(ji,jj,jk) + zgrazpof(ji,jj,jk) + zgrazfff(ji,jj,jk) ) &
-               &       - zgrazfff(ji,jj,jk) - zgrazpof(ji,jj,jk)
+               &       + unass2 * ( ferat3 * zgrazz + zgraznf + zgrazf + zgrazpof + zgrazfff ) - zgrazfff - zgrazpof
 #else
-               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zgrazpoc(ji,jj,jk)
-               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zmortz2 - zgrazffe(ji,jj,jk)
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zgrazpof(ji,jj,jk)
+               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zgrazpoc
+               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zmortz2 - zgrazffe
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zgrazpof
                tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + ferat3 * zmortz2 &
-               &       + unass2 * ( ferat3 * zgrazz(ji,jj,jk) + zgraznf(ji,jj,jk) &
-               &       + zgrazf(ji,jj,jk) + zgrazpof(ji,jj,jk) + zgrazfff(ji,jj,jk) ) &
-               &       - zgrazfff(ji,jj,jk)
+               &       + unass2 * ( ferat3 * zgrazz + zgraznf + zgrazf + zgrazpof + zgrazfff ) - zgrazfff
 #endif
 
@@ -314,5 +226,5 @@
       END DO
       !
-#if defined key_trc_diaadd && defined key_trc_dia3d && defined key_iomput
+#if defined key_diatrc && defined key_iomput
       zrfact2 = 1.e3 * rfact2r
       ! Total grazing of phyto by zoo
@@ -342,5 +254,5 @@
       !!
       !! ** Method  :   Read the nampismes namelist and check the parameters
-      !!      called at the first timestep (nittrc000)
+      !!      called at the first timestep (nit000)
       !!
       !! ** input   :   Namelist nampismes
@@ -373,4 +285,5 @@
       ENDIF
 
+
    END SUBROUTINE p4z_meso_init
 
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zmicro.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zmicro.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zmicro.F90	(revision 2528)
@@ -26,5 +26,6 @@
    PRIVATE
 
-   PUBLIC   p4z_micro    ! called in p4zbio.F90
+   PUBLIC   p4z_micro         ! called in p4zbio.F90
+   PUBLIC   p4z_micro_init    ! called in trcsms_pisces.F90
 
    !! * Shared module variables
@@ -45,12 +46,12 @@
 #  include "top_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
 CONTAINS
 
-   SUBROUTINE p4z_micro( kt,jnt )
+   SUBROUTINE p4z_micro( kt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_micro  ***
@@ -60,29 +61,19 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
+      INTEGER, INTENT(in) ::   kt ! ocean time step
       INTEGER  :: ji, jj, jk
       REAL(wp) :: zcompadi, zcompadi2, zcompaz , zcompaph, zcompapoc
-      REAL(wp) :: zgraze  , zdenom  , zdenom2
-      REAL(wp) :: zfact   , zstep   , zinano , zidiat, zipoc
+      REAL(wp) :: zgraze  , zdenom  , zdenom2, zstep
+      REAL(wp) :: zfact   , zinano , zidiat, zipoc
       REAL(wp) :: zgrarem, zgrafer, zgrapoc, zprcaca, zmortz
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrespz,ztortz
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrazp, zgrazm, zgrazsd
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrazmf, zgrazsf, zgrazpf
+      REAL(wp) :: zrespz, ztortz
+      REAL(wp) :: zgrazp, zgrazm, zgrazsd
+      REAL(wp) :: zgrazmf, zgrazsf, zgrazpf
       CHARACTER (len=25) :: charout
 
       !!---------------------------------------------------------------------
 
-      IF( ( kt * jnt ) == nittrc000  )   CALL p4z_micro_init      ! Initialization (first time-step only)
-
-      zrespz (:,:,:) = 0.
-      ztortz (:,:,:) = 0.
-      zgrazp (:,:,:) = 0.
-      zgrazm (:,:,:) = 0.
-      zgrazsd(:,:,:) = 0.
-      zgrazmf(:,:,:) = 0.
-      zgrazsf(:,:,:) = 0.
-      zgrazpf(:,:,:) = 0.
-
-#if defined key_trc_dia3d
+
+#if defined key_diatrc
       grazing(:,:,:) = 0.  !: Initialisation of  grazing
 #endif
@@ -93,33 +84,22 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-
                zcompaz = MAX( ( trn(ji,jj,jk,jpzoo) - 1.e-9 ), 0.e0 )
-# if defined key_off_degrad
-               zfact   = zstep * tgfunc(ji,jj,jk) * zcompaz *facvol(ji,jj,jk)
+# if defined key_degrad
+               zstep   = xstep * facvol(ji,jj,jk)
 # else
+               zstep   = xstep
+# endif
                zfact   = zstep * tgfunc(ji,jj,jk) * zcompaz
-# endif
-
-!     Respiration rates of both zooplankton
-!     -------------------------------------
-
-               zrespz(ji,jj,jk) = resrat * zfact  * ( 1.+ 3.* nitrfac(ji,jj,jk) )     &
+
+               !  Respiration rates of both zooplankton
+               !  -------------------------------------
+               zrespz = resrat * zfact  * ( 1.+ 3.* nitrfac(ji,jj,jk) )     &
                   &            * trn(ji,jj,jk,jpzoo) / ( xkmort + trn(ji,jj,jk,jpzoo) )
 
-!     Zooplankton mortality. A square function has been selected with
-!     no real reason except that it seems to be more stable and may
-!     mimic predation.
-!     ---------------------------------------------------------------
-               ztortz(ji,jj,jk) = mzrat * 1.e6 * zfact * trn(ji,jj,jk,jpzoo)
-
-            END DO
-         END DO
-      END DO
-
-
- 
-      DO jk = 1,jpkm1
-         DO jj = 1,jpj
-            DO ji = 1,jpi
+               !  Zooplankton mortality. A square function has been selected with
+               !  no real reason except that it seems to be more stable and may mimic predation.
+               !  ---------------------------------------------------------------
+               ztortz = mzrat * 1.e6 * zfact * trn(ji,jj,jk,jpzoo)
+
                zcompadi  = MAX( ( trn(ji,jj,jk,jpdia) - 1.e-8 ), 0.e0 )
                zcompadi2 = MIN( zcompadi, 5.e-7 )
@@ -131,9 +111,5 @@
                zdenom2 = 1./ ( xpref2p * zcompaph + xpref2c * zcompapoc + xpref2d * zcompadi2 + rtrn )
 
-               zgraze = grazrat * zstep * tgfunc(ji,jj,jk)     &
-# if defined key_off_degrad
-                  &      * facvol(ji,jj,jk)         &
-# endif
-                  &      * trn(ji,jj,jk,jpzoo)
+               zgraze = grazrat * zstep * tgfunc(ji,jj,jk) * trn(ji,jj,jk,jpzoo)
 
                zinano = xpref2p * zcompaph  * zdenom2
@@ -143,35 +119,25 @@
                zdenom = 1./ ( xkgraz + zinano * zcompaph + zipoc * zcompapoc + zidiat * zcompadi2 )
 
-               zgrazp(ji,jj,jk)  = zgraze * zinano * zcompaph * zdenom
-               zgrazm(ji,jj,jk)  = zgraze * zipoc  * zcompapoc * zdenom
-               zgrazsd(ji,jj,jk) = zgraze * zidiat * zcompadi2 * zdenom
-
-               zgrazpf (ji,jj,jk) = zgrazp(ji,jj,jk)  * trn(ji,jj,jk,jpnfe) / (trn(ji,jj,jk,jpphy) + rtrn)
-               zgrazmf(ji,jj,jk)  = zgrazm(ji,jj,jk)  * trn(ji,jj,jk,jpsfe) / (trn(ji,jj,jk,jppoc) + rtrn)
-               zgrazsf(ji,jj,jk)  = zgrazsd(ji,jj,jk) * trn(ji,jj,jk,jpdfe) / (trn(ji,jj,jk,jpdia) + rtrn)
-
-            END DO
-         END DO
-      END DO
-      
-#if defined key_trc_dia3d
-      ! Grazing by microzooplankton
-      grazing(:,:,:) = grazing(:,:,:) + zgrazp(:,:,:) + zgrazm(:,:,:) + zgrazsd(:,:,:) 
-#endif
-
-      DO jk = 1,jpkm1
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-!    Various remineralization and excretion terms
-!    --------------------------------------------
-
-               zgrarem = (  zgrazp(ji,jj,jk) + zgrazm(ji,jj,jk)  + zgrazsd(ji,jj,jk)  ) &
-                  &          * ( 1.- epsher - unass )
-               zgrafer = (  zgrazpf(ji,jj,jk) + zgrazsf(ji,jj,jk)  + zgrazmf(ji,jj,jk)  ) &
-                  &        * ( 1.- epsher - unass ) + epsher *  &
-                  &  ( zgrazm(ji,jj,jk)  * MAX((trn(ji,jj,jk,jpsfe) /(trn(ji,jj,jk,jppoc)+ rtrn)-ferat3),0.e0) &
-                  &   + zgrazp(ji,jj,jk)  * MAX((trn(ji,jj,jk,jpnfe)/(trn(ji,jj,jk,jpphy)+ rtrn)-ferat3),0.e0) &
-                  &   + zgrazsd(ji,jj,jk) * MAX((trn(ji,jj,jk,jpdfe)/(trn(ji,jj,jk,jpdia)+ rtrn)-ferat3),0.e0 )  )
-               zgrapoc = (  zgrazp(ji,jj,jk) + zgrazm(ji,jj,jk) + zgrazsd(ji,jj,jk)  ) * unass
+               zgrazp  = zgraze * zinano * zcompaph * zdenom
+               zgrazm  = zgraze * zipoc  * zcompapoc * zdenom
+               zgrazsd = zgraze * zidiat * zcompadi2 * zdenom
+
+               zgrazpf = zgrazp  * trn(ji,jj,jk,jpnfe) / (trn(ji,jj,jk,jpphy) + rtrn)
+               zgrazmf = zgrazm  * trn(ji,jj,jk,jpsfe) / (trn(ji,jj,jk,jppoc) + rtrn)
+               zgrazsf = zgrazsd * trn(ji,jj,jk,jpdfe) / (trn(ji,jj,jk,jpdia) + rtrn)
+#if defined key_diatrc
+               ! Grazing by microzooplankton
+               grazing(ji,jj,jk) = grazing(ji,jj,jk) + zgrazp + zgrazm + zgrazsd 
+#endif
+
+               !    Various remineralization and excretion terms
+               !    --------------------------------------------
+               zgrarem = ( zgrazp + zgrazm + zgrazsd ) * ( 1.- epsher - unass )
+               zgrafer = ( zgrazpf + zgrazsf + zgrazmf ) * ( 1.- epsher - unass ) &
+                  &      + epsher * ( zgrazm  * MAX((trn(ji,jj,jk,jpsfe) / (trn(ji,jj,jk,jppoc)+ rtrn)-ferat3),0.e0) & 
+                  &                 + zgrazp  * MAX((trn(ji,jj,jk,jpnfe) / (trn(ji,jj,jk,jpphy)+ rtrn)-ferat3),0.e0) &
+                  &                 + zgrazsd * MAX((trn(ji,jj,jk,jpdfe) / (trn(ji,jj,jk,jpdia)+ rtrn)-ferat3),0.e0 )  )
+
+               zgrapoc = (  zgrazp + zgrazm + zgrazsd ) 
 
                !  Update of the TRA arrays
@@ -183,42 +149,28 @@
                tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - o2ut * zgrarem * sigma1
                tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zgrafer
-               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zgrapoc
+               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zgrapoc * unass
                tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zgrarem * sigma1
 #if defined key_kriest
-               tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zgrapoc * xkr_ddiat
-#endif
-            END DO
-         END DO
-      END DO
-
-!
-!   Update the arrays TRA which contain the biological sources and sinks
-!   --------------------------------------------------------------------
-
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-
-               zmortz = ztortz(ji,jj,jk) + zrespz(ji,jj,jk)
-               tra(ji,jj,jk,jpzoo) = tra(ji,jj,jk,jpzoo) - zmortz  &
-                 &     + epsher * ( zgrazp(ji,jj,jk) + zgrazm(ji,jj,jk) + zgrazsd(ji,jj,jk))
-               tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) - zgrazp(ji,jj,jk)
-               tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) - zgrazsd(ji,jj,jk)
-               tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) - zgrazp(ji,jj,jk)  &
-                 &     * trn(ji,jj,jk,jpnch)/(trn(ji,jj,jk,jpphy)+rtrn)
-               tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) - zgrazsd(ji,jj,jk) &
-                 &     * trn(ji,jj,jk,jpdch)/(trn(ji,jj,jk,jpdia)+rtrn)
-               tra(ji,jj,jk,jpbsi) = tra(ji,jj,jk,jpbsi) - zgrazsd(ji,jj,jk) &
-                 &     * trn(ji,jj,jk,jpbsi)/(trn(ji,jj,jk,jpdia)+rtrn)
-               tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) + zgrazsd(ji,jj,jk) &
-                 &     * trn(ji,jj,jk,jpbsi)/(trn(ji,jj,jk,jpdia)+rtrn)
-               tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) - zgrazpf(ji,jj,jk)
-               tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) - zgrazsf(ji,jj,jk)
-               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zmortz - zgrazm(ji,jj,jk)
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ferat3 * zmortz   &
-                 &     + unass * ( zgrazpf(ji,jj,jk) + zgrazsf (ji,jj,jk)) &
-                 &     - (1.-unass) * zgrazmf(ji,jj,jk)
-               zprcaca = xfracal(ji,jj,jk) * unass * zgrazp(ji,jj,jk)
-#if defined key_trc_dia3d
+               tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zgrapoc * unass * xkr_ddiat
+#endif
+
+               !
+               !   Update the arrays TRA which contain the biological sources and sinks
+               !   --------------------------------------------------------------------
+
+               zmortz = ztortz + zrespz
+               tra(ji,jj,jk,jpzoo) = tra(ji,jj,jk,jpzoo) - zmortz + epsher * zgrapoc 
+               tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) - zgrazp
+               tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) - zgrazsd
+               tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) - zgrazp  * trn(ji,jj,jk,jpnch)/(trn(ji,jj,jk,jpphy)+rtrn)
+               tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) - zgrazsd * trn(ji,jj,jk,jpdch)/(trn(ji,jj,jk,jpdia)+rtrn)
+               tra(ji,jj,jk,jpbsi) = tra(ji,jj,jk,jpbsi) - zgrazsd * trn(ji,jj,jk,jpbsi)/(trn(ji,jj,jk,jpdia)+rtrn)
+               tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) + zgrazsd * trn(ji,jj,jk,jpbsi)/(trn(ji,jj,jk,jpdia)+rtrn)
+               tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) - zgrazpf
+               tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) - zgrazsf
+               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zmortz - zgrazm
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ferat3 * zmortz + unass * ( zgrazpf + zgrazsf ) - (1.-unass) * zgrazmf
+               zprcaca = xfracal(ji,jj,jk) * unass * zgrazp
+#if defined key_diatrc
                prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca  ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)
 #endif
@@ -228,5 +180,5 @@
                tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) + zprcaca
 #if defined key_kriest
-               tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + ( zmortz - zgrazm(ji,jj,jk) ) * xkr_ddiat
+               tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + ( zmortz - zgrazm ) * xkr_ddiat
 #endif
             END DO
@@ -251,5 +203,5 @@
       !!
       !! ** Method  :   Read the nampiszoo namelist and check the parameters
-      !!      called at the first timestep (nittrc000)
+      !!      called at the first timestep (nit000)
       !!
       !! ** input   :   Namelist nampiszoo
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zmort.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zmort.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zmort.F90	(revision 2528)
@@ -25,4 +25,5 @@
 
    PUBLIC   p4z_mort    
+   PUBLIC   p4z_mort_init    
 
 
@@ -35,20 +36,16 @@
      mpratm = 0.01_wp           !:
 
-   !! * Module variables
-   REAL(wp) :: zstep
-
-
 
    !!* Substitution
 #  include "top_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
 CONTAINS
 
-   SUBROUTINE p4z_mort( kt, jnt )
+   SUBROUTINE p4z_mort( kt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_mort  ***
@@ -59,10 +56,6 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
-      !!---------------------------------------------------------------------
-
-      IF( ( kt * jnt ) == nittrc000  )   CALL p4z_mort_init      ! Initialization (first time-step only)
-
-      zstep = rfact2 / rday      ! Time step duration for biology
+      INTEGER, INTENT(in) ::   kt ! ocean time step
+      !!---------------------------------------------------------------------
 
       CALL p4z_nano            ! nanophytoplankton
@@ -83,11 +76,11 @@
       INTEGER  :: ji, jj, jk
       REAL(wp) :: zcompaph
-      REAL(wp) :: zfactfe,zfactch,zprcaca,zfracal
-      REAL(wp) :: ztortp,zrespp,zmortp
+      REAL(wp) :: zfactfe, zfactch, zprcaca, zfracal
+      REAL(wp) :: ztortp , zrespp , zmortp , zstep
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
 
 
-#if defined key_trc_dia3d
+#if defined key_diatrc
      prodcal(:,:,:) = 0.  !: Initialisation of calcite production variable
 #endif
@@ -99,23 +92,17 @@
                zcompaph = MAX( ( trn(ji,jj,jk,jpphy) - 1e-8 ), 0.e0 )
 
-!     Squared mortality of Phyto similar to a sedimentation term during
-!     blooms (Doney et al. 1996)
-!     -----------------------------------------------------------------
-               zrespp = wchl * 1.e6 * zstep * xdiss(ji,jj,jk)   &
-# if defined key_off_degrad
-                  &        * facvol(ji,jj,jk)     &
+# if defined key_degrad
+               zstep =  xstep * facvol(ji,jj,jk)  
+# else
+               zstep =  xstep  
 # endif
-                  &        * zcompaph * trn(ji,jj,jk,jpphy)
-
-!     Phytoplankton mortality. This mortality loss is slightly
-!     increased when nutrients are limiting phytoplankton growth
-!     as observed for instance in case of iron limitation.
-!     ----------------------------------------------------------
-               ztortp = mprat * zstep * trn(ji,jj,jk,jpphy)          &
-# if defined key_off_degrad
-                  &          * facvol(ji,jj,jk)     &
-# endif
-                  &   / ( xkmort + trn(ji,jj,jk,jpphy) ) * zcompaph
-
+               !     Squared mortality of Phyto similar to a sedimentation term during
+               !     blooms (Doney et al. 1996)
+               zrespp = wchl * 1.e6 * zstep * xdiss(ji,jj,jk) * zcompaph * trn(ji,jj,jk,jpphy) 
+
+               !     Phytoplankton mortality. This mortality loss is slightly
+               !     increased when nutrients are limiting phytoplankton growth
+               !     as observed for instance in case of iron limitation.
+               ztortp = mprat * xstep * trn(ji,jj,jk,jpphy) / ( xkmort + trn(ji,jj,jk,jpphy) ) * zcompaph
 
                zmortp = zrespp + ztortp
@@ -130,5 +117,5 @@
                tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) - zmortp * zfactfe
                zprcaca = xfracal(ji,jj,jk) * zmortp
-#if defined key_trc_dia3d
+#if defined key_diatrc
                prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca  ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)
 #endif
@@ -169,5 +156,5 @@
       INTEGER  ::  ji, jj, jk
       REAL(wp) ::  zfactfe,zfactsi,zfactch, zcompadi
-      REAL(wp) ::  zrespp2, ztortp2, zmortp2
+      REAL(wp) ::  zrespp2, ztortp2, zmortp2, zstep
       CHARACTER (len=25) :: charout
  
@@ -175,8 +162,8 @@
 
 
-!    Aggregation term for diatoms is increased in case of nutrient
-!    stress as observed in reality. The stressed cells become more
-!    sticky and coagulate to sink quickly out of the euphotic zone
-!     ------------------------------------------------------------
+      !    Aggregation term for diatoms is increased in case of nutrient
+      !    stress as observed in reality. The stressed cells become more
+      !    sticky and coagulate to sink quickly out of the euphotic zone
+      !     ------------------------------------------------------------
 
       DO jk = 1, jpkm1
@@ -186,28 +173,27 @@
                zcompadi = MAX( ( trn(ji,jj,jk,jpdia) - 1e-8), 0. )
 
-!    Aggregation term for diatoms is increased in case of nutrient
-!    stress as observed in reality. The stressed cells become more
-!    sticky and coagulate to sink quickly out of the euphotic zone
-!     ------------------------------------------------------------
-
+               !    Aggregation term for diatoms is increased in case of nutrient
+               !    stress as observed in reality. The stressed cells become more
+               !    sticky and coagulate to sink quickly out of the euphotic zone
+               !     ------------------------------------------------------------
+
+# if defined key_degrad
+               zstep =  xstep * facvol(ji,jj,jk)  
+# else
+               zstep =  xstep  
+# endif
+               !  Phytoplankton respiration 
+               !     ------------------------
                zrespp2  = 1.e6 * zstep * (  wchl + wchld * ( 1.- xlimdia(ji,jj,jk) )  )    &
-# if defined key_off_degrad
-                  &       * facvol(ji,jj,jk)       &
-# endif
                   &       * xdiss(ji,jj,jk) * zcompadi * trn(ji,jj,jk,jpdia)
-                                                                               
-
-!     Phytoplankton mortality. 
-!     ------------------------
-               ztortp2  = mprat2 * zstep * trn(ji,jj,jk,jpdia)     &
-# if defined key_off_degrad
-                  &        * facvol(ji,jj,jk)       &
-# endif
-                  &      / ( xkmort + trn(ji,jj,jk,jpdia) ) * zcompadi
-
-                zmortp2 = zrespp2 + ztortp2
-
-!   Update the arrays tra which contains the biological sources and sinks
-!   ---------------------------------------------------------------------
+
+               !     Phytoplankton mortality. 
+               !     ------------------------
+               ztortp2  = mprat2 * zstep * trn(ji,jj,jk,jpdia)  / ( xkmort + trn(ji,jj,jk,jpdia) ) * zcompadi 
+
+               zmortp2 = zrespp2 + ztortp2
+
+               !   Update the arrays tra which contains the biological sources and sinks
+               !   ---------------------------------------------------------------------
                zfactch = trn(ji,jj,jk,jpdch) / ( trn(ji,jj,jk,jpdia) + rtrn )
                zfactfe = trn(ji,jj,jk,jpdfe) / ( trn(ji,jj,jk,jpdia) + rtrn )
@@ -249,5 +235,5 @@
       !!
       !! ** Method  :   Read the nampismort namelist and check the parameters
-      !!      called at the first timestep (nittrc000)
+      !!      called at the first timestep
       !!
       !! ** input   :   Namelist nampismort
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zopt.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zopt.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zopt.F90	(revision 2528)
@@ -16,5 +16,4 @@
    USE trc            ! tracer variables
    USE oce_trc        ! tracer-ocean share variables
-   USE trc_oce        ! ocean-tracer share variables
    USE sms_pisces     ! Source Minus Sink of PISCES
    USE iom
@@ -23,12 +22,12 @@
    PRIVATE
 
-   PUBLIC   p4z_opt   ! called in p4zbio.F90 module
+   PUBLIC   p4z_opt        ! called in p4zbio.F90 module
+   PUBLIC   p4z_opt_init   ! called in trcsms_pisces.F90 module
 
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   etot, enano, ediat   !: PAR for phyto, nano and diat 
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   emoy                 !: averaged PAR in the mixed layer
 
-   INTEGER  ::   nksrp   ! levels below which the light cannot penetrate ( depth larger than 391 m)
-   REAL(wp) ::   &
-      parlux = 0.43 / 3.e0
+   INTEGER  ::  nksrp   ! levels below which the light cannot penetrate ( depth larger than 391 m)
+   REAL(wp) ::  parlux = 0.43 / 3.e0
 
    REAL(wp), DIMENSION(3,61), PUBLIC ::   xkrgb  !: tabulated attenuation coefficients for RGB absorption
@@ -37,12 +36,12 @@
 #  include "top_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
 CONTAINS
 
-   SUBROUTINE p4z_opt(kt, jnt)
+   SUBROUTINE p4z_opt( kt, jnt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_opt  ***
@@ -54,5 +53,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
-      INTEGER  ::   ji, jj, jk, jc
+      INTEGER  ::   ji, jj, jk
       INTEGER  ::   irgb
       REAL(wp) ::   zchl, zxsi0r
@@ -64,23 +63,6 @@
 
 
-      !                                        !* tabulated attenuation coef. 
-      IF( kt * jnt == nittrc000 ) THEN
-         !                                ! level of light extinction
-         nksrp = trc_oce_ext_lev( rn_si2, 0.33e2 )
-         IF(lwp) THEN
-           WRITE(numout,*)
-           WRITE(numout,*) ' level max of computation of qsr = ', nksrp, ' ref depth = ', gdepw_0(nksrp+1), ' m'
-         ENDIF
-!!         CALL trc_oce_rgb( xkrgb )     ! tabulated attenuation coefficients
-         CALL trc_oce_rgb_read( xkrgb )     ! tabulated attenuation coefficients
-         etot (:,:,:) = 0.e0
-         enano(:,:,:) = 0.e0
-         ediat(:,:,:) = 0.e0
-         IF( ln_qsr_bio ) etot3(:,:,:) = 0.e0
-      ENDIF
-
-
-!     Initialisation of variables used to compute PAR
-!     -----------------------------------------------
+      !     Initialisation of variables used to compute PAR
+      !     -----------------------------------------------
       ze1 (:,:,jpk) = 0.e0
       ze2 (:,:,jpk) = 0.e0
@@ -227,5 +209,5 @@
       END DO
 
-#if defined key_trc_diaadd
+#if defined key_diatrc
 # if ! defined key_iomput
       ! save for outputs
@@ -243,4 +225,24 @@
    END SUBROUTINE p4z_opt
 
+   SUBROUTINE p4z_opt_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE p4z_opt_init  ***
+      !!
+      !! ** Purpose :   Initialization of tabulated attenuation coef
+      !!
+      !!
+      !!----------------------------------------------------------------------
+
+      CALL trc_oce_rgb( xkrgb )                  ! tabulated attenuation coefficients
+!!      CALL trc_oce_rgb_read( xkrgb )               ! tabulated attenuation coefficients
+      nksrp = trc_oce_ext_lev( r_si2, 0.33e2 )     ! max level of light extinction (Blue Chl=0.01)
+      IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksrp, ' ref depth = ', gdepw_0(nksrp+1), ' m'
+      !
+                         etot (:,:,:) = 0.e0
+                         enano(:,:,:) = 0.e0
+                         ediat(:,:,:) = 0.e0
+      IF( ln_qsr_bio )   etot3(:,:,:) = 0.e0
+      ! 
+   END SUBROUTINE p4z_opt_init
 #else
    !!----------------------------------------------------------------------
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zprod.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zprod.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zprod.F90	(revision 2528)
@@ -23,9 +23,11 @@
 
    USE lib_mpp
+   USE lib_fortran
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   p4z_prod    ! called in p4zbio.F90
+   PUBLIC   p4z_prod         ! called in p4zbio.F90
+   PUBLIC   p4z_prod_init    ! called in trcsms_pisces.F90
 
    !! * Shared module variables
@@ -41,21 +43,18 @@
      grosip    = 0.151_wp
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk)  ::        &
-     &                   prmax
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk)  ::  prmax 
    
    REAL(wp) ::   &
+      rday1                      ,  &  !: 0.6 / rday
       texcret                    ,  &  !: 1 - excret 
       texcret2                   ,  &  !: 1 - excret2        
-      rpis180                    ,  &  !: rpi / 180
       tpp                              !: Total primary production
-
-   INTEGER  ::  nspyr                  !: number of timesteps per year
 
    !!* Substitution
 #  include "top_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -78,6 +77,6 @@
       REAL(wp) ::   zmxltst, zmxlday, zlim1
       REAL(wp) ::   zpislopen  , zpislope2n
-      REAL(wp) ::   zrum, zcodel, zargu, zvol
-#if defined key_trc_diaadd && defined key_trc_dia3d
+      REAL(wp) ::   zrum, zcodel, zargu, zval, zvol
+#if defined key_diatrc
       REAL(wp) ::   zrfact2
 #endif
@@ -90,8 +89,4 @@
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
-
-
-      IF( ( kt * jnt ) == nittrc000  )   CALL p4z_prod_init      ! Initialization (first time-step only)
-
 
       zprorca (:,:,:) = 0.0
@@ -109,21 +104,13 @@
       ! Computation of the optimal production
 
-# if defined key_off_degrad
-      prmax(:,:,:) = 0.6 / rday * tgfunc(:,:,:) * facvol(:,:,:)
+# if defined key_degrad
+      prmax(:,:,:) = rday1 * tgfunc(:,:,:) * facvol(:,:,:)
 # else
-      prmax(:,:,:) = 0.6 / rday * tgfunc(:,:,:)
+      prmax(:,:,:) = rday1 * tgfunc(:,:,:)
 # endif
 
       ! compute the day length depending on latitude and the day
-      IF(lwp) write(numout,*)
-      IF(lwp) write(numout,*) 'p4zday : - Julian day ', nday_year
-      IF(lwp) write(numout,*) '~~~~~~'
-
-      IF( nleapy == 1 .AND. MOD( nyear, 4 ) == 0 ) THEN
-         zrum = FLOAT( nday_year - 80 ) / 366.
-      ELSE
-         zrum = FLOAT( nday_year - 80 ) / 365.
-      ENDIF
-      zcodel = ASIN(  SIN( zrum * rpi * 2. ) * SIN( rpis180 * 23.5 )  )
+      zrum = FLOAT( nday_year - 80 ) / REAL(nyear_len(1), wp)
+      zcodel = ASIN(  SIN( zrum * rpi * 2. ) * SIN( rad * 23.5 )  )
 
       ! day length in hours
@@ -131,7 +118,9 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rpis180 )
+            zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
             zargu = MAX( -1., MIN(  1., zargu ) )
-            zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rpis180 / 15. )
+            zval  = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
+            IF( zval < 1.e0 )   zval = 24.
+            zstrn(ji,jj) = 24. / zval
          END DO
       END DO
@@ -147,5 +136,5 @@
                ! 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. )
+                   ztn    = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )
                    zadap  = 0.+ 1.* ztn / ( 2.+ ztn )
                    zadap2 = 0.e0
@@ -227,7 +216,4 @@
       END DO
 
-
-      WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24.
-      zstrn(:,:) = 24. / zstrn(:,:)
 
 !CDIR NOVERRCHK
@@ -331,27 +317,18 @@
 
      ! Total primary production per year
-     DO jk = 1, jpkm1
-        DO jj = 1, jpj
-          DO ji = 1, jpi
-             zvol = cvol(ji,jj,jk)
-#if defined key_off_degrad
-             zvol = zvol * facvol(ji,jj,jk)
+
+#if defined key_degrad
+     tpp = tpp + glob_sum( ( zprorca(:,:,:) + zprorcad(:,:,:) ) * cvol(:,:,:) * facvol(:,:,:) )
+#else
+     tpp = tpp + glob_sum( ( zprorca(:,:,:) + zprorcad(:,:,:) ) * cvol(:,:,:) )
 #endif
-             tpp  = tpp + ( zprorca(ji,jj,jk) + zprorcad(ji,jj,jk) ) &
-                          * zvol * tmask(ji,jj,jk) * tmask_i(ji,jj)
-          END DO
-        END DO
-      END DO
-
-
-      IF( MOD( kt, nspyr ) == 0 .AND. jnt == nrdttrc ) THEN
-        IF( lk_mpp ) CALL mpp_sum( tpp )
-        WRITE(numout,*) 'Total PP :'
+
+     IF( kt == nitend .AND. jnt == nrdttrc ) THEN
+        WRITE(numout,*) 'Total PP (Gtc) :'
         WRITE(numout,*) '-------------------- : ',tpp * 12. / 1.E12
-        WRITE(numout,*) '(GtC/yr)'
-        tpp = 0.
+        WRITE(numout,*) 
       ENDIF
 
-#if defined key_trc_diaadd && defined key_trc_dia3d && ! defined key_iomput
+#if defined key_diatrc && ! defined key_iomput
       !   Supplementary diagnostics
       zrfact2 = 1.e3 * rfact2r
@@ -367,5 +344,5 @@
 #endif
 
-#if defined key_trc_diaadd && defined key_trc_dia3d && defined key_iomput
+#if defined key_diatrc && defined key_iomput
       zrfact2 = 1.e3 * rfact2r
       IF ( jnt == nrdttrc ) then
@@ -396,5 +373,5 @@
       !!
       !! ** Method  :   Read the nampisprod namelist and check the parameters
-      !!      called at the first timestep (nittrc000)
+      !!      called at the first timestep (nit000)
       !!
       !! ** input   :   Namelist nampisprod
@@ -423,8 +400,5 @@
       ENDIF
 
-      ! number of timesteps per year
-      nspyr  = INT( nyear_len(1) * rday / rdt )
-
-      rpis180   = rpi / 180.
+      rday1     = 0.6 / rday 
       texcret   = 1.0 - excret
       texcret2  = 1.0 - excret2
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zrem.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zrem.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zrem.F90	(revision 2528)
@@ -27,5 +27,6 @@
    PRIVATE
 
-   PUBLIC   p4z_rem    ! called in p4zbio.F90
+   PUBLIC   p4z_rem         ! called in p4zbio.F90
+   PUBLIC   p4z_rem_init    ! called in trcsms_pisces.F90
 
    !! * Shared module variables
@@ -41,18 +42,16 @@
      &                   denitr                     !: denitrification array
 
-   REAL(wp) ::   &
-     xstep            !: Time step duration for biology
 
    !!* Substitution
 #  include "top_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
 CONTAINS
 
-   SUBROUTINE p4z_rem(kt, jnt)
+   SUBROUTINE p4z_rem( kt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_rem  ***
@@ -62,5 +61,5 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
+      INTEGER, INTENT(in) ::   kt ! ocean time step
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   zremip, zremik , zlam1b
@@ -72,5 +71,5 @@
       REAL(wp) ::   zofer2, zdenom, zdenom2
 #endif
-      REAL(wp) ::   zlamfac, zonitr
+      REAL(wp) ::   zlamfac, zonitr, zstep
       REAL(wp), DIMENSION(jpi,jpj)     ::   ztempbac
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdepbac, zfesatur, zolimi
@@ -78,12 +77,4 @@
 
       !!---------------------------------------------------------------------
-
-
-      IF( ( kt * jnt ) == nittrc000  )  THEN
-         CALL p4z_rem_init                ! Initialization (first time-step only)
-         xstep = rfact2 / rday            ! Time step duration for the biology
-         nitrfac(:,:,:) = 0.0
-         denitr (:,:,:) = 0.0  
-      ENDIF
 
 
@@ -94,7 +85,7 @@
        ztempbac(:,:)   = 0.0
 
-!      Computation of the mean phytoplankton concentration as
-!      a crude estimate of the bacterial biomass
-!      --------------------------------------------------
+      !  Computation of the mean phytoplankton concentration as
+      !  a crude estimate of the bacterial biomass
+      !   --------------------------------------------------
 
       DO jk = 1, jpkm1
@@ -114,39 +105,33 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-
-!    DENITRIFICATION FACTOR COMPUTED FROM O2 LEVELS
-!    ----------------------------------------------
-
+               ! denitrification factor computed from O2 levels
                nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - trn(ji,jj,jk,jpoxy) )    &
                   &                                / ( oxymin + trn(ji,jj,jk,jpoxy) )  )
-            END DO
-         END DO
-      END DO
-
-      nitrfac(:,:,:) = MIN( 1., nitrfac(:,:,:) )
-
-
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-
-!     DOC ammonification. Depends on depth, phytoplankton biomass
-!     and a limitation term which is supposed to be a parameterization
-!     of the bacterial activity. 
-!     ----------------------------------------------------------------
-               zremik = xremik * xstep / 1.e-6 * xlimbac(ji,jj,jk)         &
-# if defined key_off_degrad
-                  &            * facvol(ji,jj,jk)              &
+               nitrfac(ji,jj,jk) = MIN( 1., nitrfac(ji,jj,jk) )
+            END DO
+         END DO
+      END DO
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+# if defined key_degrad
+               zstep = xstep * facvol(ji,jj,jk)
+# else
+               zstep = xstep
 # endif
-                  &            * zdepbac(ji,jj,jk)
+               ! DOC ammonification. Depends on depth, phytoplankton biomass
+               !     and a limitation term which is supposed to be a parameterization
+               !     of the bacterial activity. 
+               zremik = xremik * zstep / 1.e-6 * xlimbac(ji,jj,jk) * zdepbac(ji,jj,jk) 
                zremik = MAX( zremik, 5.5e-4 * xstep )
 
-!     Ammonification in oxic waters with oxygen consumption
-!     -----------------------------------------------------
+               !     Ammonification in oxic waters with oxygen consumption
+               !     -----------------------------------------------------
                zolimi(ji,jj,jk) = MIN(  ( trn(ji,jj,jk,jpoxy) - rtrn ) / o2ut,  &
                   &                    zremik * ( 1.- nitrfac(ji,jj,jk) ) * trn(ji,jj,jk,jpdoc)  ) 
 
-!     Ammonification in suboxic waters with denitrification
-!     -------------------------------------------------------
+               !     Ammonification in suboxic waters with denitrification
+               !     -------------------------------------------------------
                denitr(ji,jj,jk) = MIN(  ( trn(ji,jj,jk,jpno3) - rtrn ) / rdenit,   &
                   &                     zremik * nitrfac(ji,jj,jk) * trn(ji,jj,jk,jpdoc)  )
@@ -167,22 +152,21 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-
-!    NH4 nitrification to NO3. Ceased for oxygen concentrations
-!    below 2 umol/L. Inhibited at strong light 
-!    ----------------------------------------------------------
-               zonitr  = nitrif * xstep * trn(ji,jj,jk,jpnh4) / ( 1.+ emoy(ji,jj,jk) )     &
-# if defined key_off_degrad
-                  &      * facvol(ji,jj,jk)              &
+# if defined key_degrad
+               zstep = xstep * facvol(ji,jj,jk)
+# else
+               zstep = xstep
 # endif
-                  &      * ( 1.- nitrfac(ji,jj,jk) )
-
-!
-!   Update of the tracers trends
-!   ----------------------------
-
-              tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) - zonitr
-              tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) + zonitr
-              tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - o2nit * zonitr
-              tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - rno3  * zonitr
+               !    NH4 nitrification to NO3. Ceased for oxygen concentrations
+               !    below 2 umol/L. Inhibited at strong light 
+               !    ----------------------------------------------------------
+               zonitr  = nitrif * zstep * trn(ji,jj,jk,jpnh4) / ( 1.+ emoy(ji,jj,jk) ) * ( 1.- nitrfac(ji,jj,jk) ) 
+
+               !   Update of the tracers trends
+               !   ----------------------------
+
+               tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) - zonitr
+               tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) + zonitr
+               tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - o2nit * zonitr
+               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - rno3  * zonitr
 
             END DO
@@ -200,11 +184,11 @@
             DO ji = 1, jpi
 
-!    Bacterial uptake of iron. No iron is available in DOC. So
-!    Bacteries are obliged to take up iron from the water. Some
-!    studies (especially at Papa) have shown this uptake to be
-!    significant
-!    ----------------------------------------------------------
+               !    Bacterial uptake of iron. No iron is available in DOC. So
+               !    Bacteries are obliged to take up iron from the water. Some
+               !    studies (especially at Papa) have shown this uptake to be significant
+               !    ----------------------------------------------------------
                zbactfer = 15.e-6 * rfact2 * 4.* 0.4 * prmax(ji,jj,jk)           &
-                  &               * ( xlimphy(ji,jj,jk) * zdepbac(ji,jj,jk))**2           &
+                  &               * ( xlimphy(ji,jj,jk) * zdepbac(ji,jj,jk))           &
+                  &               * ( xlimphy(ji,jj,jk) * zdepbac(ji,jj,jk))           &
                   &                  / ( xkgraz2 + zdepbac(ji,jj,jk) )                    &
                   &                  * ( 0.5 + SIGN( 0.5, trn(ji,jj,jk,jpfer) -2.e-11 )  )
@@ -216,5 +200,4 @@
                tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zbactfer
 #endif
-
             END DO
          END DO
@@ -230,18 +213,18 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-
-!    POC disaggregation by turbulence and bacterial activity. 
-!    -------------------------------------------------------------
-               zremip = xremip * xstep * tgfunc(ji,jj,jk)   &
-# if defined key_off_degrad
-                  &            * facvol(ji,jj,jk)              &
+# if defined key_degrad
+               zstep = xstep * facvol(ji,jj,jk)
+# else
+               zstep = xstep
 # endif
-                  &            * ( 1.- 0.5 * nitrfac(ji,jj,jk) )
-
-!    POC disaggregation rate is reduced in anoxic zone as shown by
-!    sediment traps data. In oxic area, the exponent of the martin s
-!    law is around -0.87. In anoxic zone, it is around -0.35. This
-!    means a disaggregation constant about 0.5 the value in oxic zones
-!    -----------------------------------------------------------------
+               !    POC disaggregation by turbulence and bacterial activity. 
+               !    -------------------------------------------------------------
+               zremip = xremip * zstep * tgfunc(ji,jj,jk) * ( 1.- 0.5 * nitrfac(ji,jj,jk) ) 
+
+               !    POC disaggregation rate is reduced in anoxic zone as shown by
+               !    sediment traps data. In oxic area, the exponent of the martin s
+               !    law is around -0.87. In anoxic zone, it is around -0.35. This
+               !    means a disaggregation constant about 0.5 the value in oxic zones
+               !    -----------------------------------------------------------------
                zorem  = zremip * trn(ji,jj,jk,jppoc)
                zofer  = zremip * trn(ji,jj,jk,jpsfe)
@@ -253,6 +236,6 @@
 #endif
 
-!  Update the appropriate tracers trends
-!  -------------------------------------
+               !  Update the appropriate tracers trends
+               !  -------------------------------------
 
                tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zorem
@@ -282,21 +265,20 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-
-!     Remineralization rate of BSi depedant on T and saturation
-!     ---------------------------------------------------------
+# if defined key_degrad
+               zstep = xstep * facvol(ji,jj,jk)
+# else
+               zstep = xstep
+# endif
+               !     Remineralization rate of BSi depedant on T and saturation
+               !     ---------------------------------------------------------
                zsatur  = ( sio3eq(ji,jj,jk) - trn(ji,jj,jk,jpsil) ) / ( sio3eq(ji,jj,jk) + rtrn )
                zsatur  = MAX( rtrn, zsatur )
                zsatur2 = zsatur * ( 1. + tn(ji,jj,jk) / 400.)**4
                znusil  = 0.225  * ( 1. + tn(ji,jj,jk) / 15.) * zsatur + 0.775 * zsatur2**9
-#    if defined key_off_degrad
-               zsiremin = xsirem * xstep * znusil * facvol(ji,jj,jk)
-# else
-               zsiremin = xsirem * xstep * znusil
-#    endif
+               zsiremin = xsirem * zstep * znusil
                zosil = zsiremin * trn(ji,jj,jk,jpdsi)
 
                tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) - zosil
                tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) + zosil
-
                !
             END DO
@@ -317,24 +299,25 @@
 !CDIR NOVERRCHK
             DO ji = 1, jpi
-!
-!      Compute de different ratios for scavenging of iron
-!      --------------------------------------------------
+# if defined key_degrad
+               zstep = xstep * facvol(ji,jj,jk)
+# else
+               zstep = xstep
+# endif
+               !  Compute de different ratios for scavenging of iron
+               !  --------------------------------------------------
 
 #if  defined key_kriest
-                zdenom1 = trn(ji,jj,jk,jppoc) / &
+               zdenom1 = trn(ji,jj,jk,jppoc) / &
            &           ( trn(ji,jj,jk,jppoc) + trn(ji,jj,jk,jpdsi) + trn(ji,jj,jk,jpcal) + rtrn )
 #else
-                zdenom = 1. / ( trn(ji,jj,jk,jppoc) + trn(ji,jj,jk,jpgoc)  &
+               zdenom = 1. / ( trn(ji,jj,jk,jppoc) + trn(ji,jj,jk,jpgoc)  &
            &            + trn(ji,jj,jk,jpdsi) + trn(ji,jj,jk,jpcal) + rtrn )
 
-                zdenom1 = trn(ji,jj,jk,jppoc) * zdenom
-                zdenom2 = trn(ji,jj,jk,jpgoc) * zdenom
-#endif
-
-
-!     scavenging rate of iron. this scavenging rate depends on the
-!     load in particles on which they are adsorbed. The
-!     parameterization has been taken from studies on Th
-!     ------------------------------------------------------------
+               zdenom1 = trn(ji,jj,jk,jppoc) * zdenom
+               zdenom2 = trn(ji,jj,jk,jpgoc) * zdenom
+#endif
+               !  scavenging rate of iron. this scavenging rate depends on the load in particles
+               !  on which they are adsorbed. The  parameterization has been taken from studies on Th
+               !     ------------------------------------------------------------
                zkeq = fekeq(ji,jj,jk)
                zfeequi = ( -( 1. + zfesatur(ji,jj,jk) * zkeq - zkeq * trn(ji,jj,jk,jpfer) )               &
@@ -349,15 +332,10 @@
                   &                      + trn(ji,jj,jk,jpcal) + trn(ji,jj,jk,jpdsi)  ) * 1.e6
 #endif
-
-# if defined key_off_degrad
-               zscave = zfeequi * zlam1b * xstep  * facvol(ji,jj,jk)
-# else
-               zscave = zfeequi * zlam1b * xstep
-# endif
-
-!  Increased scavenging for very high iron concentrations
-!  found near the coasts due to increased lithogenic particles
-!  and let s say it unknown processes (precipitation, ...)
-!  -----------------------------------------------------------
+               zscave = zfeequi * zlam1b * zstep
+
+               !  Increased scavenging for very high iron concentrations
+               !  found near the coasts due to increased lithogenic particles
+               !  and let s say it unknown processes (precipitation, ...)
+               !  -----------------------------------------------------------
                zlamfac = MAX( 0.e0, ( gphit(ji,jj) + 55.) / 30. )
                zlamfac = MIN( 1.  , zlamfac )
@@ -374,9 +352,5 @@
 #endif
 
-# if defined key_off_degrad
-               zaggdfe = zlam1b * xstep * 0.5 * ( trn(ji,jj,jk,jpfer) - zfeequi ) * facvol(ji,jj,jk)
-# else
-               zaggdfe = zlam1b * xstep * 0.5 * ( trn(ji,jj,jk,jpfer) - zfeequi )
-# endif
+               zaggdfe = zlam1b * zstep * 0.5 * ( trn(ji,jj,jk,jpfer) - zfeequi )
 
                tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfe
@@ -400,6 +374,6 @@
        ENDIF
 
-!     Update the arrays TRA which contain the biological sources and sinks
-!     --------------------------------------------------------------------
+       !     Update the arrays TRA which contain the biological sources and sinks
+       !     --------------------------------------------------------------------
 
       DO jk = 1, jpkm1
@@ -429,5 +403,5 @@
       !!
       !! ** Method  :   Read the nampisrem namelist and check the parameters
-      !!      called at the first timestep (nittrc000)
+      !!      called at the first timestep
       !!
       !! ** input   :   Namelist nampisrem
@@ -452,9 +426,8 @@
       ENDIF
 
+      nitrfac(:,:,:) = 0.0
+      denitr (:,:,:) = 0.0  
+
    END SUBROUTINE p4z_rem_init
-
-
-
-
 
 #else
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zsed.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zsed.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zsed.F90	(revision 2528)
@@ -19,4 +19,5 @@
    USE sms_pisces
    USE lib_mpp
+   USE lib_fortran
    USE prtctl_trc
    USE p4zbio
@@ -34,4 +35,5 @@
 
    PUBLIC   p4z_sed   
+   PUBLIC   p4z_sed_init   
 
    !! * Shared module variables
@@ -47,18 +49,19 @@
 
    !! * Module variables
-   INTEGER ::                   &
-     ryyss,                     &  !: number of seconds per year
-     rmtss                         !: number of seconds per month
-
+   REAL(wp) :: ryyss               !: number of seconds per year 
+   REAL(wp) :: ryyss1              !: inverse of ryyss
+   REAL(wp) :: rmtss               !: number of seconds per month
+   REAL(wp) :: rday1               !: inverse of rday
+
+   INTEGER , PARAMETER :: &
+        jpmth = 12, jpyr = 1
    INTEGER ::                   &
       numdust,                  &  !: logical unit for surface fluxes data
       nflx1 , nflx2,            &  !: first and second record used
       nflx11, nflx12      ! ???
-   REAL(wp), DIMENSION(jpi,jpj,2) ::    &  !:
-     dustmo                                !: 2 consecutive set of dust fields 
-   REAL(wp), DIMENSION(jpi,jpj)   ::    &
-     rivinp, cotdep, nitdep, dust
-   REAL(wp), DIMENSION(jpi,jpj,jpk)  ::   &
-     ironsed
+   REAL(wp), DIMENSION(jpi,jpj,jpmth) ::  dustmo    !: set of dust fields
+   REAL(wp), DIMENSION(jpi,jpj)      ::  rivinp, cotdep, nitdep, dust 
+   REAL(wp), DIMENSION(jpi,jpj)      ::  e1e2t
+   REAL(wp), DIMENSION(jpi,jpj,jpk)  ::  ironsed 
    REAL(wp) :: sumdepsi, rivalkinput, rivpo4input, nitdepinput
 
@@ -66,12 +69,12 @@
 #  include "top_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
    !! $Header:$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
 CONTAINS
 
-   SUBROUTINE p4z_sed(kt, jnt)
+   SUBROUTINE p4z_sed( kt, jnt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_sed  ***
@@ -84,30 +87,18 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
-      INTEGER  ::   ji, jj, jk
-      INTEGER  ::   ikt
+      INTEGER  ::   ji, jj, jk, ikt
 #if ! defined key_sed
       REAL(wp) ::   zsumsedsi, zsumsedpo4, zsumsedcal
+      REAL(wp) ::   zrivalk, zrivsil, zrivpo4
 #endif
-      REAL(wp) ::   zconctmp , zdenitot  , znitrpottot
-      REAL(wp) ::   zlim, zconctmp2, zstep, zfact
+      REAL(wp) ::   zdenitot, znitrpottot, zlim, zfact
+      REAL(wp) ::   zwsbio3, zwsbio4, zwscal
       REAL(wp), DIMENSION(jpi,jpj)     ::   zsidep
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zwork, zwork1
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   znitrpot, zirondep
-#if defined key_diaadd || defined key_trc_dia3d 
-      REAL(wp) :: zrfact2
-# if defined key_iomput
-     REAL(wp), DIMENSION(jpi,jpj)    ::    zw2d 
-# endif
-#endif
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
 
-
-      IF( ( kt * jnt ) == nittrc000   )   CALL p4z_sed_init      ! Initialization (first time-step only)
-      IF( (jnt == 1) .and. ( ln_dustfer ) )  CALL p4z_sbc( kt )
-
-      zstep = rfact2 / rday      ! Time step duration for the biology
-
-      zirondep(:,:,:) = 0.e0          ! Initialisation of variables used to compute deposition
-      zsidep  (:,:)   = 0.e0
+      IF( jnt == 1  .AND.  ln_dustfer  )  CALL p4z_sbc( kt )
 
       ! Iron and Si deposition at the surface
@@ -116,5 +107,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            zirondep(ji,jj,1) = ( dustsolub * dust(ji,jj) / ( 55.85 * rmtss ) + 3.e-10 / ryyss )   &
+            zirondep(ji,jj,1) = ( dustsolub * dust(ji,jj) / ( 55.85 * rmtss ) + 3.e-10 * ryyss1 )   &
                &             * rfact2 / fse3t(ji,jj,1)
             zsidep  (ji,jj)   = 8.8 * 0.075 * dust(ji,jj) * rfact2 / ( fse3t(ji,jj,1) * 28.1 * rmtss )
@@ -150,9 +141,4 @@
 
 #if ! defined key_sed
-      ! Initialisation of variables used to compute Sinking Speed
-      zsumsedsi  = 0.e0
-      zsumsedpo4 = 0.e0
-      zsumsedcal = 0.e0
-
       ! Loss of biogenic silicon, Caco3 organic carbon in the sediments. 
       ! First, the total loss is computed.
@@ -161,24 +147,23 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            ikt = MAX( mbathy(ji,jj)-1, 1 )
-            zfact = e1t(ji,jj) * e2t(ji,jj) / rday * tmask_i(ji,jj)
+            ikt = mbkt(ji,jj) 
 # if defined key_kriest
-            zsumsedsi  = zsumsedsi  + zfact *  trn(ji,jj,ikt,jpdsi) * wscal (ji,jj,ikt)
-            zsumsedpo4 = zsumsedpo4 + zfact *  trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt)
+            zwork (ji,jj) = trn(ji,jj,ikt,jpdsi) * wscal (ji,jj,ikt)
+            zwork1(ji,jj) = trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt)
 # else
-            zsumsedsi  = zsumsedsi  + zfact *  trn(ji,jj,ikt,jpdsi) * wsbio4(ji,jj,ikt)
-            zsumsedpo4 = zsumsedpo4 + zfact *( trn(ji,jj,ikt,jpgoc) * wsbio4(ji,jj,ikt)   &
-               &       + trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt) )
+            zwork (ji,jj) = trn(ji,jj,ikt,jpdsi) * wsbio4(ji,jj,ikt)
+            zwork1(ji,jj) = trn(ji,jj,ikt,jpgoc) * wsbio4(ji,jj,ikt) + 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
-
-      IF( lk_mpp ) THEN
-         CALL mpp_sum( zsumsedsi  )   ! sums over the global domain
-         CALL mpp_sum( zsumsedcal )   ! sums over the global domain
-         CALL mpp_sum( zsumsedpo4 )   ! sums over the global domain
-      ENDIF
-
+         END DO
+      END DO
+      zsumsedsi  = glob_sum( zwork (:,:) * e1e2t(:,:) ) * rday1
+      zsumsedpo4 = glob_sum( zwork1(:,:) * e1e2t(:,:) ) * rday1
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ikt = mbkt(ji,jj) 
+            zwork (ji,jj) = trn(ji,jj,ikt,jpcal) * wscal (ji,jj,ikt)
+         END DO
+      END DO
+      zsumsedcal = glob_sum( zwork (:,:) * e1e2t(:,:) ) * 2.0 * rday1
 #endif
 
@@ -191,70 +176,50 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            ikt = MAX( mbathy(ji,jj) - 1, 1 )
-            zconctmp = trn(ji,jj,ikt,jpdsi) * zstep / fse3t(ji,jj,ikt)   &
-# if ! defined key_kriest
-     &             * wscal (ji,jj,ikt)
+            ikt = mbkt(ji,jj)
+            zfact = xstep / fse3t(ji,jj,ikt)
+            zwsbio3 = 1._wp - zfact * wsbio3(ji,jj,ikt)
+            zwsbio4 = 1._wp - zfact * wsbio4(ji,jj,ikt)
+            zwscal  = 1._wp - zfact * wscal (ji,jj,ikt)
+            !
+# if defined key_kriest
+            trn(ji,jj,ikt,jpdsi) = trn(ji,jj,ikt,jpdsi) * zwsbio4
+            trn(ji,jj,ikt,jpnum) = trn(ji,jj,ikt,jpnum) * zwsbio4
+            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) * zwsbio3
+            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) * zwsbio3
 # else
-     &             * wsbio4(ji,jj,ikt)
+            trn(ji,jj,ikt,jpdsi) = trn(ji,jj,ikt,jpdsi) * zwscal 
+            trn(ji,jj,ikt,jpgoc) = trn(ji,jj,ikt,jpgoc) * zwsbio4
+            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) * zwsbio3
+            trn(ji,jj,ikt,jpbfe) = trn(ji,jj,ikt,jpbfe) * zwsbio4
+            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) * zwsbio3
 # endif
-            trn(ji,jj,ikt,jpdsi) = trn(ji,jj,ikt,jpdsi) - zconctmp
+            trn(ji,jj,ikt,jpcal) = trn(ji,jj,ikt,jpcal) * zwscal
+         END DO
+      END DO
 
 #if ! defined key_sed
-            trn(ji,jj,ikt,jpsil) = trn(ji,jj,ikt,jpsil) + zconctmp   &
-            &      * ( 1.- ( sumdepsi + rivalkinput / ryyss / 6. ) / zsumsedsi )
-#endif
-         END DO
-      END DO
-
+      zrivsil =  1._wp - ( sumdepsi + rivalkinput * ryyss1 / 6. ) / zsumsedsi 
+      zrivalk =  1._wp - ( rivalkinput * ryyss1 ) / zsumsedcal 
+      zrivpo4 =  1._wp - ( rivpo4input * ryyss1 ) / zsumsedpo4 
       DO jj = 1, jpj
          DO ji = 1, jpi
-            ikt = MAX( mbathy(ji,jj) - 1, 1 )
-            zconctmp = trn(ji,jj,ikt,jpcal) * wscal(ji,jj,ikt) * zstep / fse3t(ji,jj,ikt)
-            trn(ji,jj,ikt,jpcal) = trn(ji,jj,ikt,jpcal) - zconctmp
-
-#if ! defined key_sed
-            trn(ji,jj,ikt,jptal) = trn(ji,jj,ikt,jptal) + zconctmp   &
-               &   * ( 1.- ( rivalkinput / ryyss ) / zsumsedcal ) * 2.e0
-            trn(ji,jj,ikt,jpdic) = trn(ji,jj,ikt,jpdic) + zconctmp   &
-               &   * ( 1.- ( rivalkinput / ryyss ) / zsumsedcal )
-#endif
-         END DO
-      END DO
-
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ikt = MAX( mbathy(ji,jj) - 1, 1 )
-            zfact = zstep / fse3t(ji,jj,ikt)
-# if ! defined key_kriest
-            zconctmp  = trn(ji,jj,ikt,jpgoc)
-            zconctmp2 = trn(ji,jj,ikt,jppoc)
-            trn(ji,jj,ikt,jpgoc) = trn(ji,jj,ikt,jpgoc) - zconctmp  * wsbio4(ji,jj,ikt) * zfact
-            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) - zconctmp2 * wsbio3(ji,jj,ikt) * zfact
-#if ! defined key_sed
-            trn(ji,jj,ikt,jpdoc) = trn(ji,jj,ikt,jpdoc)    &
-            &      + ( zconctmp  * wsbio4(ji,jj,ikt) + zconctmp2 * wsbio3(ji,jj,ikt) ) * zfact   &
-            &      * ( 1.- rivpo4input / (ryyss * zsumsedpo4 ) )
-#endif
-            trn(ji,jj,ikt,jpbfe) = trn(ji,jj,ikt,jpbfe) - trn(ji,jj,ikt,jpbfe) * wsbio4(ji,jj,ikt) * zfact
-            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) - trn(ji,jj,ikt,jpsfe) * wsbio3(ji,jj,ikt) * zfact
-
+            ikt = mbkt(ji,jj)
+            zfact = xstep / fse3t(ji,jj,ikt)
+            zwsbio3 = zfact * wsbio3(ji,jj,ikt)
+            zwsbio4 = zfact * wsbio4(ji,jj,ikt)
+            zwscal  = zfact * wscal (ji,jj,ikt)
+            trn(ji,jj,ikt,jptal) =  trn(ji,jj,ikt,jptal) + trn(ji,jj,ikt,jpcal) * zwscal  * zrivalk * 2.0
+            trn(ji,jj,ikt,jpdic) =  trn(ji,jj,ikt,jpdic) + trn(ji,jj,ikt,jpcal) * zwscal  * zrivalk
+# if defined key_kriest
+            trn(ji,jj,ikt,jpsil) =  trn(ji,jj,ikt,jpsil) + trn(ji,jj,ikt,jpdsi) * zwsbio4 * zrivsil 
+            trn(ji,jj,ikt,jpdoc) =  trn(ji,jj,ikt,jpdoc) + trn(ji,jj,ikt,jppoc) * zwsbio3 * zrivpo4 
 # else
-            zconctmp  = trn(ji,jj,ikt,jpnum)
-            zconctmp2 = trn(ji,jj,ikt,jppoc)
-            trn(ji,jj,ikt,jpnum) = trn(ji,jj,ikt,jpnum)   &
-            &      - zconctmp  * wsbio4(ji,jj,ikt) * zfact
-            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc)   &
-            &      - zconctmp2 * wsbio3(ji,jj,ikt) * zfact
-#if ! defined key_sed
-            trn(ji,jj,ikt,jpdoc) = trn(ji,jj,ikt,jpdoc)    &
-            &      + ( zconctmp2 * wsbio3(ji,jj,ikt) )   &
-            &      * zfact * ( 1.- rivpo4input / ( ryyss * zsumsedpo4 ) )
-#endif
-            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe)   &
-            &      - trn(ji,jj,ikt,jpsfe) * wsbio3(ji,jj,ikt) * zfact
-
+            trn(ji,jj,ikt,jpsil) =  trn(ji,jj,ikt,jpsil) + trn(ji,jj,ikt,jpdsi) * zwscal  * zrivsil 
+            trn(ji,jj,ikt,jpdoc) =  trn(ji,jj,ikt,jpdoc)   &
+            &                     + ( trn(ji,jj,ikt,jppoc) * zwsbio3 + trn(ji,jj,ikt,jpgoc) * zwsbio4 ) * zrivpo4
 # endif
          END DO
       END DO
+# endif
 
       ! Nitrogen fixation (simple parameterization). The total gain
@@ -263,14 +228,5 @@
       ! -------------------------------------------------------------
 
-      zdenitot = 0.e0
-      DO jk = 1, jpkm1
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-               zdenitot = zdenitot + denitr(ji,jj,jk) * rdenit * cvol(ji,jj,jk) * xnegtr(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-
-      IF( lk_mpp )   CALL mpp_sum( zdenitot )      ! sum over the global domain
+      zdenitot = glob_sum( denitr(:,:,:)  * cvol(:,:,:) * xnegtr(:,:,:) ) * rdenit
 
       ! Potential nitrogen fixation dependant on temperature and iron
@@ -285,6 +241,6 @@
                zlim = ( 1.- xnanono3(ji,jj,jk) - xnanonh4(ji,jj,jk) )
                IF( zlim <= 0.2 )   zlim = 0.01
-               znitrpot(ji,jj,jk) = MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) / rday )   &
-# if defined key_off_degrad
+               znitrpot(ji,jj,jk) = MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) * rday1 )   &
+# if defined key_degrad
                &                  * facvol(ji,jj,jk)   &
 # endif
@@ -295,14 +251,5 @@
       END DO
 
-      znitrpottot = 0.e0
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               znitrpottot = znitrpottot + znitrpot(ji,jj,jk) * cvol(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-
-      IF( lk_mpp )   CALL mpp_sum( znitrpottot )  ! sum over the global domain
+      znitrpottot = glob_sum( znitrpot(:,:,:) * cvol(:,:,:) )
 
       ! Nitrogen change due to nitrogen fixation
@@ -312,10 +259,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-# if ! defined key_c1d && ( defined key_orca_r4 || defined key_orca_r2 || defined key_orca_r05 || defined key_orca_r025 )
-!!             zfact = znitrpot(ji,jj,jk) * zdenitot / znitrpottot
                zfact = znitrpot(ji,jj,jk) * 1.e-7
-# else
-               zfact = znitrpot(ji,jj,jk) * 1.e-7
-# endif
                trn(ji,jj,jk,jpnh4) = trn(ji,jj,jk,jpnh4) + zfact
                trn(ji,jj,jk,jpoxy) = trn(ji,jj,jk,jpoxy) + zfact   * o2nit
@@ -325,18 +267,18 @@
       END DO
 
-#if defined key_trc_diaadd || defined key_trc_dia3d
-      zrfact2 = 1.e+3 * rfact2r
+#if defined key_diatrc
+      zfact = 1.e+3 * rfact2r
 #  if  ! defined key_iomput
-      trc2d(:,:,jp_pcs0_2d + 11) = zirondep(:,:,1)         * zrfact2 * fse3t(:,:,1) * tmask(:,:,1)
-      trc2d(:,:,jp_pcs0_2d + 12) = znitrpot(:,:,1) * 1.e-7 * zrfact2 * fse3t(:,:,1) * tmask(:,:,1)
-# else
-      ! surface downward net flux of iron
-      zw2d(:,:)   =  ( zirondep(:,:,1) + ironsed(:,:,1) * rfact2 ) * zrfact2 * fse3t(:,:,1) * tmask(:,:,1) 
-      IF( jnt == nrdttrc ) CALL iom_put( "Irondep", zw2d )
-      ! nitrogen fixation at surface
-      zw2d(:,:)   =  znitrpot(:,:,1) * 1.e-7 * zrfact2  * fse3t(:,:,1) * tmask(:,:,1)
-      IF( jnt == nrdttrc ) CALL iom_put( "Nfix" , zw2d )
-# endif
-# endif
+      trc2d(:,:,jp_pcs0_2d + 11) = zirondep(:,:,1)         * zfact * fse3t(:,:,1) * tmask(:,:,1)
+      trc2d(:,:,jp_pcs0_2d + 12) = znitrpot(:,:,1) * 1.e-7 * zfact * fse3t(:,:,1) * tmask(:,:,1)
+#  else
+      zwork (:,:)  =  ( zirondep(:,:,1) + ironsed(:,:,1) * rfact2 ) * zfact * fse3t(:,:,1) * tmask(:,:,1) 
+      zwork1(:,:)  =  znitrpot(:,:,1) * 1.e-7                       * zfact * fse3t(:,:,1) * tmask(:,:,1)
+      IF( jnt == nrdttrc ) THEN
+         CALL iom_put( "Irondep", zwork  )  ! surface downward net flux of iron
+         CALL iom_put( "Nfix"   , zwork1 )  ! nitrogen fixation at surface
+      ENDIF
+#  endif
+#endif
       !
        IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
@@ -348,5 +290,5 @@
    END SUBROUTINE p4z_sed
 
-   SUBROUTINE p4z_sbc(kt)
+   SUBROUTINE p4z_sbc( kt )
 
       !!----------------------------------------------------------------------
@@ -365,10 +307,6 @@
 
       !! * Local declarations
-      INTEGER ::   &
-         imois, imois2,       &  ! temporary integers
-         i15  , iman             !    "          "
-      REAL(wp) ::   &
-         zxy                     !    "         "
-
+      INTEGER :: imois, i15, iman 
+      REAL(wp) :: zxy
 
       !!---------------------------------------------------------------------
@@ -381,28 +319,9 @@
       imois = nmonth + i15 - 1
       IF( imois == 0 ) imois = iman
-      imois2 = nmonth
-
-      ! 1. first call kt=nit000
-      ! -----------------------
-
-      IF( kt == nit000 ) THEN
-         ! initializations
-         nflx1  = 0
-         nflx11 = 0
-         ! open the file
-         IF(lwp) THEN
-            WRITE(numout,*) ' '
-            WRITE(numout,*) ' **** Routine p4z_sbc'
-         ENDIF
-         CALL iom_open ( 'dust.orca.nc', numdust )
-      ENDIF
-
-
-     ! Read monthly file
-      ! ----------------
-
+
+      ! Calendar computation
       IF( kt == nit000 .OR. imois /= nflx1 ) THEN
 
-         ! Calendar computation
+         IF( kt == nit000 )  nflx1  = 0
 
          ! nflx1 number of the first file record used in the simulation
@@ -410,36 +329,20 @@
 
          nflx1 = imois
-         nflx2 = nflx1+1
+         nflx2 = nflx1 + 1
          nflx1 = MOD( nflx1, iman )
          nflx2 = MOD( nflx2, iman )
          IF( nflx1 == 0 )   nflx1 = iman
          IF( nflx2 == 0 )   nflx2 = iman
-         IF(lwp) WRITE(numout,*) 'first record file used nflx1 ',nflx1
-         IF(lwp) WRITE(numout,*) 'last  record file used nflx2 ',nflx2
-
-         ! Read monthly fluxes data
-
-         ! humidity
-         CALL iom_get ( numdust, jpdom_data, 'dust', dustmo(:,:,1), nflx1 )
-         CALL iom_get ( numdust, jpdom_data, 'dust', dustmo(:,:,2), nflx2 )
-
-         IF(lwp .AND. nitend-nit000 <= 100 ) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) ' read clio flx ok'
-            WRITE(numout,*)
-               WRITE(numout,*)
-               WRITE(numout,*) 'Clio month: ',nflx1,'  field: dust'
-               CALL prihre( dustmo(:,:,1),jpi,jpj,1,jpi,20,1,jpj,10,1e9,numout )
-         ENDIF
+         IF(lwp) WRITE(numout,*) 
+         IF(lwp) WRITE(numout,*) ' p4z_sbc : first record file used nflx1 ',nflx1
+         IF(lwp) WRITE(numout,*) ' p4z_sbc : last  record file used nflx2 ',nflx2
 
       ENDIF
 
-     ! 3. at every time step interpolation of fluxes
+      ! 3. at every time step interpolation of fluxes
       ! ---------------------------------------------
 
       zxy = FLOAT( nday + 15 - 30 * i15 ) / 30
-      dust(:,:) = ( (1.-zxy) * dustmo(:,:,1) + zxy * dustmo(:,:,2) )
-
-      IF( kt == nitend ) CALL iom_close (numdust)
+      dust(:,:) = ( (1.-zxy) * dustmo(:,:,nflx1) + zxy * dustmo(:,:,nflx2) )
 
    END SUBROUTINE p4z_sbc
@@ -454,5 +357,5 @@
       !!
       !! ** Method  :   Read the files and compute the budget
-      !!      called at the first timestep (nittrc000)
+      !!      called at the first timestep (nit000)
       !!
       !! ** input   :   external netcdf files
@@ -460,6 +363,5 @@
       !!----------------------------------------------------------------------
 
-      INTEGER ::   ji, jj, jk, jm
-      INTEGER , PARAMETER ::   jpmois = 12, jpan = 1
+      INTEGER :: ji, jj, jk, jm
       INTEGER :: numriv, numbath, numdep
 
@@ -469,5 +371,4 @@
       REAL(wp) , DIMENSION (jpi,jpj)     ::   riverdoc, river, ndepo
       REAL(wp) , DIMENSION (jpi,jpj,jpk) ::   cmask
-      REAL(wp) , DIMENSION(jpi,jpj,12)    ::   zdustmo
 
       NAMELIST/nampissed/ ln_dustfer, ln_river, ln_ndepo, ln_sedinput, sedfeinput, dustsolub
@@ -495,10 +396,10 @@
          IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
          CALL iom_open ( 'dust.orca.nc', numdust )
-         DO jm = 1, jpmois
-            CALL iom_get( numdust, jpdom_data, 'dust', zdustmo(:,:,jm), jm )
+         DO jm = 1, jpmth
+            CALL iom_get( numdust, jpdom_data, 'dust', dustmo(:,:,jm), jm )
          END DO
          CALL iom_close( numdust )
       ELSE
-         zdustmo(:,:,:) = 0.e0
+         dustmo(:,:,:) = 0.e0
          dust(:,:) = 0.0
       ENDIF
@@ -510,6 +411,6 @@
          IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
          CALL iom_open ( 'river.orca.nc', numriv )
-         CALL iom_get  ( numriv, jpdom_data, 'riverdic', river   (:,:), jpan )
-         CALL iom_get  ( numriv, jpdom_data, 'riverdoc', riverdoc(:,:), jpan )
+         CALL iom_get  ( numriv, jpdom_data, 'riverdic', river   (:,:), jpyr )
+         CALL iom_get  ( numriv, jpdom_data, 'riverdoc', riverdoc(:,:), jpyr )
          CALL iom_close( numriv )
       ELSE
@@ -524,5 +425,5 @@
          IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
          CALL iom_open ( 'ndeposition.orca.nc', numdep )
-         CALL iom_get  ( numdep, jpdom_data, 'ndep', ndepo(:,:), jpan )
+         CALL iom_get  ( numdep, jpdom_data, 'ndep', ndepo(:,:), jpyr )
          CALL iom_close( numdep )
       ELSE
@@ -537,5 +438,5 @@
          IF(lwp) WRITE(numout,*) '       from bathy.orca.nc file '
          CALL iom_open ( 'bathy.orca.nc', numbath )
-         CALL iom_get  ( numbath, jpdom_data, 'bathy', cmask(:,:,:), jpan )
+         CALL iom_get  ( numbath, jpdom_data, 'bathy', cmask(:,:,:), jpyr )
          CALL iom_close( numbath )
          !
@@ -546,5 +447,5 @@
                      zmaskt = tmask(ji+1,jj,jk) * tmask(ji-1,jj,jk) * tmask(ji,jj+1,jk)    &
                         &                       * tmask(ji,jj-1,jk) * tmask(ji,jj,jk+1)
-                     IF( zmaskt == 0. )   cmask(ji,jj,jk ) = 0.1
+                     IF( zmaskt == 0. )   cmask(ji,jj,jk ) = MAX( 0.1, cmask(ji,jj,jk) ) 
                   ENDIF
                END DO
@@ -567,20 +468,19 @@
 
 
-      ! Number of seconds per year and per month
-      ryyss = nyear_len(1) * rday
-      rmtss = ryyss / raamo
+      !                                    ! Number of seconds per year and per month
+      ryyss  = nyear_len(1) * rday
+      rmtss  = ryyss / raamo
+      rday1  = 1. / rday
+      ryyss1 = 1. / ryyss
+      !                                    ! ocean surface cell
+      e1e2t(:,:) = e1t(:,:) * e2t(:,:)
 
       ! total atmospheric supply of Si
       ! ------------------------------
       sumdepsi = 0.e0
-      DO jm = 1, jpmois
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               sumdepsi = sumdepsi + zdustmo(ji,jj,jm) / (12.*rmtss) * 8.8        &
-                  &     * 0.075/28.1 * e1t(ji,jj) * e2t(ji,jj) * tmask(ji,jj,1) * tmask_i(ji,jj)
-            END DO
-         END DO
-      END DO
-      IF( lk_mpp )  CALL mpp_sum( sumdepsi )  ! sum over the global domain
+      DO jm = 1, jpmth
+         zcoef = 1. / ( 12. * rmtss ) * 8.8 * 0.075 / 28.1        
+         sumdepsi = sumdepsi + glob_sum( dustmo(:,:,jm) * e1e2t(:,:) ) * zcoef
+      ENDDO
 
       ! N/P and Si releases due to coastal rivers
@@ -588,5 +488,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            zcoef = ryyss * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1) * tmask(ji,jj,1) * tmask_i(ji,jj)
+            zcoef = ryyss * e1e2t(ji,jj)  * fse3t(ji,jj,1) * tmask(ji,jj,1) 
             cotdep(ji,jj) =  river(ji,jj)                  *1E9 / ( 12. * zcoef + rtrn )
             rivinp(ji,jj) = (river(ji,jj)+riverdoc(ji,jj)) *1E9 / ( 31.6* zcoef + rtrn )
@@ -597,20 +497,7 @@
       CALL lbc_lnk( cotdep , 'T', 1. )  ;  CALL lbc_lnk( rivinp , 'T', 1. )  ;  CALL lbc_lnk( nitdep , 'T', 1. )
 
-      rivpo4input = 0.e0
-      rivalkinput = 0.e0
-      nitdepinput = 0.e0
-      DO jj = 2 , jpjm1
-         DO ji = fs_2, fs_jpim1
-            zcoef = cvol(ji,jj,1) * ryyss
-            rivpo4input = rivpo4input + rivinp(ji,jj) * zcoef
-            rivalkinput = rivalkinput + cotdep(ji,jj) * zcoef
-            nitdepinput = nitdepinput + nitdep(ji,jj) * zcoef
-         END DO
-     END DO
-      IF( lk_mpp ) THEN
-         CALL mpp_sum( rivpo4input )  ! sum over the global domain
-         CALL mpp_sum( rivalkinput )  ! sum over the global domain
-         CALL mpp_sum( nitdepinput )  ! sum over the global domain
-      ENDIF
+      rivpo4input = glob_sum( rivinp(:,:) * cvol(:,:,1) ) * ryyss
+      rivalkinput = glob_sum( cotdep(:,:) * cvol(:,:,1) ) * ryyss
+      nitdepinput = glob_sum( nitdep(:,:) * cvol(:,:,1) ) * ryyss
 
 
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zsink.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zsink.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zsink.F90	(revision 2528)
@@ -19,5 +19,6 @@
    PRIVATE
 
-   PUBLIC   p4z_sink    ! called in p4zbio.F90
+   PUBLIC   p4z_sink         ! called in p4zbio.F90
+   PUBLIC   p4z_sink_init    ! called in trcsms_pisces.F90
 
    !! * Shared module variables
@@ -31,7 +32,4 @@
      sinkcal, sinksil,    &    !: CaCO3 and BSi sinking fluxes
      sinkfer                   !: Small BFe sinking flux
-
-   REAL(wp) ::   &
-     xstep , xstep2            !: Time step duration for biology
 
    INTEGER  :: &
@@ -71,7 +69,7 @@
 #  include "top_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -97,5 +95,5 @@
       REAL(wp) :: zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5
       REAL(wp) :: zval1, zval2, zval3, zval4
-#if defined key_trc_diaadd
+#if defined key_diatrc
       REAL(wp) :: zrfact2
       INTEGER  :: ik1
@@ -106,12 +104,6 @@
       !!---------------------------------------------------------------------
 
-      IF( ( kt * jnt ) == nittrc000  )  THEN 
-          CALL p4z_sink_init   ! Initialization (first time-step only)
-          xstep  = rfact2 / rday      ! Time step duration for biology
-          xstep2 = rfact2 /  2.
-      ENDIF
-
-!     Initialisation of variables used to compute Sinking Speed
-!     ---------------------------------------------------------
+      !     Initialisation of variables used to compute Sinking Speed
+      !     ---------------------------------------------------------
 
        znum3d(:,:,:) = 0.e0
@@ -120,6 +112,6 @@
        zval3 = 1. + xkr_eta
 
-!     Computation of the vertical sinking speed : Kriest et Evans, 2000
-!     -----------------------------------------------------------------
+     !     Computation of the vertical sinking speed : Kriest et Evans, 2000
+     !     -----------------------------------------------------------------
 
       DO jk = 1, jpkm1
@@ -128,9 +120,9 @@
                IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
                   znum = trn(ji,jj,jk,jppoc) / ( trn(ji,jj,jk,jpnum) + rtrn ) / xkr_massp
-! -------------- To avoid sinking speed over 50 m/day -------
+                  ! -------------- To avoid sinking speed over 50 m/day -------
                   znum  = MIN( xnumm(jk), znum )
                   znum  = MAX( 1.1      , znum )
                   znum3d(ji,jj,jk) = znum
-!------------------------------------------------------------
+                  !------------------------------------------------------------
                   zeps  = ( zval1 * znum - 1. )/ ( znum - 1. )
                   zfm   = xkr_frac**( 1. - zeps )
@@ -150,7 +142,6 @@
       wscal(:,:,:) = MAX( wsbio3(:,:,:), 50. )
 
-
-!   INITIALIZE TO ZERO ALL THE SINKING ARRAYS
-!   -----------------------------------------
+      !   INITIALIZE TO ZERO ALL THE SINKING ARRAYS
+      !   -----------------------------------------
 
       sinking (:,:,:) = 0.e0
@@ -160,7 +151,6 @@
       sinksil (:,:,:) = 0.e0
 
-!   Compute the sedimentation term using p4zsink2 for all
-!   the sinking particles
-!   -----------------------------------------------------
+     !   Compute the sedimentation term using p4zsink2 for all the sinking particles
+     !   -----------------------------------------------------
 
       CALL p4z_sink2( wsbio3, sinking , jppoc )
@@ -170,7 +160,6 @@
       CALL p4z_sink2( wscal , sinkcal , jpcal )
 
-!  Exchange between organic matter compartments due to
-!  coagulation/disaggregation
-!  ---------------------------------------------------
+     !  Exchange between organic matter compartments due to coagulation/disaggregation
+     !  ---------------------------------------------------
 
       zval1 = 1. + xkr_zeta
@@ -185,8 +174,8 @@
 
                   znum = trn(ji,jj,jk,jppoc)/(trn(ji,jj,jk,jpnum)+rtrn) / xkr_massp
-! -------------- To avoid sinking speed over 50 m/day -------
+                  !-------------- To avoid sinking speed over 50 m/day -------
                   znum  = min(xnumm(jk),znum)
                   znum  = MAX( 1.1,znum)
-!------------------------------------------------------------
+                  !------------------------------------------------------------
                   zeps  = ( zval1 * znum - 1.) / ( znum - 1.)
                   zdiv  = MAX( 1.e-4, ABS( zeps - zval3) ) * SIGN( 1., zeps - zval3 )
@@ -199,6 +188,6 @@
                   zsm   = xkr_frac**xkr_eta
 
-!    Part I : Coagulation dependant on turbulence
-!    ----------------------------------------------
+                  !    Part I : Coagulation dependant on turbulence
+                  !    ----------------------------------------------
 
                   zagg1 = ( 0.163 * trn(ji,jj,jk,jpnum)**2               &
@@ -207,5 +196,5 @@
                      &            * (zfm*xkr_mass_max**2-xkr_mass_min**2)                  &
                      &            * (zeps-1.)**2/(zdiv2*zdiv3))            &
-# if defined key_off_degrad
+# if defined key_degrad
                      &                 *facvol(ji,jj,jk)       &
 # endif
@@ -219,5 +208,5 @@
                      &                    -zfm*xkr_mass_max**3*(1.+3.*((zeps-1.)/           &
                      &                    (zeps-2.)+(zeps-1.)/zdiv3)+(zeps-1.)/zdiv1))      &
-#    if defined key_off_degrad
+#    if defined key_degrad
                      &                 *facvol(ji,jj,jk)             &
 #    endif
@@ -225,5 +214,5 @@
 
                   zagg3 = (  0.163*trn(ji,jj,jk,jpnum)**2*zfm**2*8. * xkr_mass_max**3   &
-#    if defined key_off_degrad
+#    if defined key_degrad
                      &                 *facvol(ji,jj,jk)             &
 #    endif
@@ -232,7 +221,7 @@
                   zaggsh = ( zagg1 + zagg2 + zagg3 ) * rfact2 * xdiss(ji,jj,jk) / 1000.
 
-!    Aggregation of small into large particles
-!    Part II : Differential settling
-!    ----------------------------------------------
+                 !    Aggregation of small into large particles
+                 !    Part II : Differential settling
+                 !    ----------------------------------------------
 
                   zagg4 = (  2.*3.141*0.125*trn(ji,jj,jk,jpnum)**2*                       &
@@ -242,5 +231,5 @@
                      &                 ((zfm*zfm*xkr_mass_max**2*zsm-xkr_mass_min**2)     &
                      &                 *xkr_eta)/(zdiv*zdiv3*zdiv5) )                     &
-# if defined key_off_degrad
+# if defined key_degrad
                      &                 *facvol(ji,jj,jk)        &
 # endif
@@ -252,5 +241,5 @@
                      &                 /zdiv3-(xkr_mass_min**2-zfm*zsm*xkr_mass_max**2)    &
                      &                 /zdiv)                   &
-# if defined key_off_degrad
+# if defined key_degrad
                      &                 *facvol(ji,jj,jk)        &
 # endif
@@ -261,10 +250,10 @@
                   zagg = 0.5 * xkr_stick * ( zaggsh + zaggsi )
 
-!     Aggregation of DOC to small particles
-!     --------------------------------------
+                  !     Aggregation of DOC to small particles
+                  !     --------------------------------------
 
                   zaggdoc = ( 0.4 * trn(ji,jj,jk,jpdoc)               &
                      &        + 1018.  * trn(ji,jj,jk,jppoc)  ) * xstep    &
-# if defined key_off_degrad
+# if defined key_degrad
                      &        * facvol(ji,jj,jk)                              &
 # endif
@@ -281,5 +270,5 @@
       END DO
 
-#if defined key_trc_diaadd
+#if defined key_diatrc
       zrfact2 = 1.e3 * rfact2r
       ik1 = iksed + 1
@@ -332,5 +321,5 @@
       !!
       !! ** Method  :   Read the nampiskrs namelist and check the parameters
-      !!      called at the first timestep (nittrc000)
+      !!      called at the first timestep 
       !!
       !! ** input   :   Namelist nampiskrs
@@ -473,6 +462,6 @@
       REAL(wp) ::   zagg1, zagg2, zagg3, zagg4
       REAL(wp) ::   zagg , zaggfe, zaggdoc, zaggdoc2
-      REAL(wp) ::   zfact, zwsmax
-#if defined key_trc_dia3d
+      REAL(wp) ::   zfact, zwsmax, zstep
+#if defined key_diatrc
       REAL(wp) ::   zrfact2
       INTEGER  ::   ik1
@@ -481,16 +470,11 @@
       !!---------------------------------------------------------------------
 
-      IF( ( kt * jnt ) == nittrc000  )  THEN
-          xstep  = rfact2 / rday      ! Timestep duration for biology
-          xstep2 = rfact2 /  2.
-      ENDIF
-
-!    Sinking speeds of detritus is increased with depth as shown
-!    by data and from the coagulation theory
-!    -----------------------------------------------------------
+      !    Sinking speeds of detritus is increased with depth as shown
+      !    by data and from the coagulation theory
+      !    -----------------------------------------------------------
       DO jk = 1, jpkm1
          DO jj = 1, jpj
             DO ji=1,jpi
-               zfact = MAX( 0., fsdepw(ji,jj,jk+1)-hmld(ji,jj) ) / 4000.
+               zfact = MAX( 0., fsdepw(ji,jj,jk+1) - hmld(ji,jj) ) / 4000.
                wsbio4(ji,jj,jk) = wsbio2 + ( 200.- wsbio2 ) * zfact
             END DO
@@ -498,15 +482,13 @@
       END DO
 
-!      LIMIT THE VALUES OF THE SINKING SPEEDS 
-!      TO AVOID NUMERICAL INSTABILITIES
-
+      ! limit the values of the sinking speeds to avoid numerical instabilities  
       wsbio3(:,:,:) = wsbio
-!
-! OA Below, this is garbage. the ideal would be to find a time-splitting
-! OA algorithm that does not increase the computing cost by too much
-! OA In ROMS, I have included a time-splitting procedure. But it is 
-! OA too expensive as the loop is computed globally. Thus, a small e3t
-! OA at one place determines the number of subtimesteps globally
-! OA AWFULLY EXPENSIVE !! Not able to find a better approach. Damned !!
+      !
+      ! OA Below, this is garbage. the ideal would be to find a time-splitting 
+      ! OA algorithm that does not increase the computing cost by too much
+      ! OA In ROMS, I have included a time-splitting procedure. But it is 
+      ! OA too expensive as the loop is computed globally. Thus, a small e3t
+      ! OA at one place determines the number of subtimesteps globally
+      ! OA AWFULLY EXPENSIVE !! Not able to find a better approach. Damned !!
 
       DO jk = 1,jpkm1
@@ -522,6 +504,6 @@
       wscal(:,:,:) = wsbio4(:,:,:)
 
-!   INITIALIZE TO ZERO ALL THE SINKING ARRAYS
-!   -----------------------------------------
+      !  Initializa to zero all the sinking arrays 
+      !   -----------------------------------------
 
       sinking (:,:,:) = 0.e0
@@ -532,7 +514,6 @@
       sinkfer2(:,:,:) = 0.e0
 
-!   Compute the sedimentation term using p4zsink2 for all
-!   the sinking particles
-!   -----------------------------------------------------
+      !   Compute the sedimentation term using p4zsink2 for all the sinking particles
+      !   -----------------------------------------------------
 
       CALL p4z_sink2( wsbio3, sinking , jppoc )
@@ -543,45 +524,33 @@
       CALL p4z_sink2( wscal , sinkcal , jpcal )
 
-!  Exchange between organic matter compartments due to
-!  coagulation/disaggregation
-!  ---------------------------------------------------
+      !  Exchange between organic matter compartments due to coagulation/disaggregation
+      !  ---------------------------------------------------
 
       DO jk = 1, jpkm1
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zfact = xstep * xdiss(ji,jj,jk)
+# if defined key_degrad
+               zstep = xstep * facvol(ji,jj,jk)
+# else
+               zstep = xstep 
+# endif
+               zfact = zstep * xdiss(ji,jj,jk)
                !  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)
                zagg2 = 1.054e4 * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
-# endif
 
                ! 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)
-               zagg4 = 0.e0 * xstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
-# endif
+               zagg3 = 0.66 * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
+               zagg4 = 0.e0 * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
+
                zagg   = zagg1 + zagg2 + zagg3 + zagg4
                zaggfe = zagg * trn(ji,jj,jk,jpsfe) / ( trn(ji,jj,jk,jppoc) + rtrn )
 
                ! Aggregation of DOC to small particles
-#if defined key_off_degrad
-               zaggdoc = ( 80.* trn(ji,jj,jk,jpdoc) + 698. * trn(ji,jj,jk,jppoc) )       &
-                  &      * facvol(ji,jj,jk)  * zfact * trn(ji,jj,jk,jpdoc)
-               zaggdoc2 = 1.05e4 * zfact * trn(ji,jj,jk,jpgoc)   &
-                  &      * 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)
+               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
@@ -595,5 +564,5 @@
       END DO
 
-#if defined key_trc_diaadd
+#if defined key_diatrc
       zrfact2 = 1.e3 * rfact2r
       ik1  = iksed + 1
@@ -623,4 +592,10 @@
    END SUBROUTINE p4z_sink
 
+   SUBROUTINE p4z_sink_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE p4z_sink_init  ***
+      !!----------------------------------------------------------------------
+   END SUBROUTINE p4z_sink_init
+
 #endif
 
@@ -641,5 +616,5 @@
       !!
       INTEGER  ::   ji, jj, jk, jn
-      REAL(wp) ::   zigma,zew,zign, zflx
+      REAL(wp) ::   zigma,zew,zign, zflx, zstep
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  ztraz, zakz
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwsink2
@@ -647,9 +622,11 @@
 
 
+      zstep = rfact2 / 2.
+
       ztraz(:,:,:) = 0.e0
       zakz (:,:,:) = 0.e0
 
       DO jk = 1, jpkm1
-# if defined key_off_degrad
+# if defined key_degrad
          zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rday * tmask(:,:,jk+1) * facvol(:,:,jk)
 # else
@@ -693,7 +670,7 @@
             DO jj = 1, jpj      
                DO ji = 1, jpi    
-                  zigma = zwsink2(ji,jj,jk+1) * xstep2 / fse3w(ji,jj,jk+1)
+                  zigma = zwsink2(ji,jj,jk+1) * zstep / fse3w(ji,jj,jk+1)
                   zew   = zwsink2(ji,jj,jk+1)
-                  psinkflx(ji,jj,jk+1) = -zew * ( trn(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * xstep2
+                  psinkflx(ji,jj,jk+1) = -zew * ( trn(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
                END DO
             END DO
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/par_pisces.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/par_pisces.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/par_pisces.F90	(revision 2528)
@@ -6,7 +6,7 @@
    !! History :   2.0  !  2007-12  (C. Ethe, G. Madec)  revised architecture
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
    USE par_lobster, ONLY : jp_lobster      !: number of tracers in LOBSTER
@@ -17,8 +17,8 @@
    IMPLICIT NONE
 
-   INTEGER, PARAMETER ::   jp_lp      = jp_lobster      !: cumulative number of already defined TRC
-   INTEGER, PARAMETER ::   jp_lp_2d   = jp_lobster_2d   !:
-   INTEGER, PARAMETER ::   jp_lp_3d   = jp_lobster_3d   !:
-   INTEGER, PARAMETER ::   jp_lp_trd  = jp_lobster_trd  !:
+   INTEGER, PUBLIC, PARAMETER ::   jp_lp      = jp_lobster      !: cumulative number of already defined TRC
+   INTEGER, PUBLIC, PARAMETER ::   jp_lp_2d   = jp_lobster_2d   !:
+   INTEGER, PUBLIC, PARAMETER ::   jp_lp_3d   = jp_lobster_3d   !:
+   INTEGER, PUBLIC, PARAMETER ::   jp_lp_trd  = jp_lobster_trd  !:
 
 #if defined key_pisces  &&  defined key_kriest
@@ -29,6 +29,6 @@
    LOGICAL, PUBLIC, PARAMETER ::   lk_kriest     = .TRUE.  !: Kriest flag 
    INTEGER, PUBLIC, PARAMETER ::   jp_pisces     =  23     !: number of passive tracers
-   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_2d  =  13     !: additional 2d output ('key_trc_diaadd')
-   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_3d  =  18     !: additional 3d output ('key_trc_diaadd')
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_2d  =  13     !: additional 2d output ('key_diatrc')
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_3d  =  18     !: additional 3d output ('key_diatrc')
    INTEGER, PUBLIC, PARAMETER ::   jp_pisces_trd =   1     !: number of sms trends for PISCES
 
@@ -67,6 +67,6 @@
    LOGICAL, PUBLIC, PARAMETER ::   lk_kriest     = .FALSE. !: Kriest flag 
    INTEGER, PUBLIC, PARAMETER ::   jp_pisces     = 24      !: number of PISCES passive tracers
-   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_2d  = 13      !: additional 2d output ('key_trc_diaadd')
-   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_3d  = 11      !: additional 3d output ('key_trc_diaadd')
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_2d  = 13      !: additional 2d output ('key_diatrc')
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_3d  = 11      !: additional 3d output ('key_diatrc')
    INTEGER, PUBLIC, PARAMETER ::   jp_pisces_trd =  1      !: number of sms trends for PISCES
 
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90	(revision 2528)
@@ -23,4 +23,5 @@
    REAL(wp) ::   rfact , rfactr    !: ???
    REAL(wp) ::   rfact2, rfact2r   !: ???
+   REAL(wp) ::   xstep             !: Time step duration for biology
 
    !!*  Biological parameters 
@@ -62,5 +63,5 @@
    REAL(wp), DIMENSION(jpi,jpj,jpk) ::   xlimbac    !: ??
    REAL(wp), DIMENSION(jpi,jpj,jpk) ::   xdiss      !: ??
-#if defined key_trc_dia3d
+#if defined key_diatrc
    REAL(wp), DIMENSION(jpi,jpj,jpk) ::   prodcal    !: Calcite production
    REAL(wp), DIMENSION(jpi,jpj,jpk) ::   grazing    !: Total zooplankton grazing
@@ -91,7 +92,7 @@
    
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 3.2 , LOCEAN-IPSL (2009) 
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!======================================================================   
 END MODULE sms_pisces    
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcctl_pisces.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcctl_pisces.F90	(revision 2524)
+++ 	(revision )
@@ -1,60 +1,0 @@
-MODULE trcctl_pisces
-   !!======================================================================
-   !!                         ***  trcctl_pisces.F90  ***
-   !! TOP :                Control of PISCES biogeochemical model
-   !!======================================================================
-   !!----------------------------------------------------------------------
-   !! History :   1.0  !  2000-12 (C. Ethe) assign a parameter to name individual tracers
-   !!----------------------------------------------------------------------
-
-#if defined key_pisces
-
-   USE oce_trc
-   USE trc
-
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC trc_ctl_pisces     ! called by ???
-
-
-   !!----------------------------------------------------------------------
-   !! NEMO/TOP 1.0 , LOCEAN-IPSL (2005) 
-   !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
-CONTAINS
-
-   SUBROUTINE trc_ctl_pisces
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE trc_ctl_pisces  ***
-      !!
-      !! ** Purpose :   control the cpp options, namelist and files 
-      !!----------------------------------------------------------------------
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' use PISCES biological model '
-
-   ! Check number of tracers
-   ! -----------------------
-#if  defined key_kriest
-      IF( jp_pisces /= 23) CALL ctl_stop( ' PISCES must have 23 passive tracers. Change jp_pisces in par_pisces.F90' )
-#else
-      IF( jp_pisces /= 24) CALL ctl_stop( ' PISCES must have 24 passive tracers. Change jp_pisces in par_pisces.F90' )
-#endif
-
-   END SUBROUTINE trc_ctl_pisces
-
-#else
-   !!----------------------------------------------------------------------
-   !!  Empty module :                                            No PISCES
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE trc_ctl_pisces                 ! Dummy routine
-   END SUBROUTINE trc_ctl_pisces
-#endif
-
-   !!======================================================================
-END MODULE trcctl_pisces
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 2528)
@@ -40,7 +40,7 @@
 #  include "top_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -55,17 +55,20 @@
 
 
+      !  Control consitency
+      CALL trc_ctl_pisces
+
+
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) ' trc_ini_pisces :   PISCES biochemical model initialisation'
       IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
 
-
       !                                            ! Time-step
-      rfact   = rdttra(1) * float(ndttrc)          ! ---------
+      rfact   = rdttrc(1)                          ! ---------
       rfactr  = 1. / rfact
-      rfact2  = rfact / float(nrdttrc)
+      rfact2  = rfact / FLOAT( nrdttrc )
       rfact2r = 1. / rfact2
 
-      IF(lwp) WRITE(numout,*) '    Tracer  time step    rfact  = ', rfact, ' rdt = ', rdt
-      IF(lwp) write(numout,*) '    Biology time step    rfact2 = ', rfact2
+      IF(lwp) WRITE(numout,*) '    Passive Tracer  time step    rfact  = ', rfact, ' rdt = ', rdttra(1)
+      IF(lwp) write(numout,*) '    PISCES  Biology time step    rfact2 = ', rfact2
 
 
@@ -80,6 +83,4 @@
 
       CALL p4z_che        ! initialize the chemical constants
-
-      ndayflxtr = 0      !  Initialize a counter for the computation of chemistry
 
       ! Initialization of tracer concentration in case of  no restart 
@@ -128,5 +129,25 @@
       !
    END SUBROUTINE trc_ini_pisces
-   
+ 
+   SUBROUTINE trc_ctl_pisces
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_ctl_pisces  ***
+      !!
+      !! ** Purpose :   control the cpp options, namelist and files 
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' use PISCES biological model '
+
+   ! Check number of tracers
+   ! -----------------------
+#if  defined key_kriest
+      IF( jp_pisces /= 23) CALL ctl_stop( ' PISCES must have 23 passive tracers. Change jp_pisces in par_pisces.F90' )
+#else
+      IF( jp_pisces /= 24) CALL ctl_stop( ' PISCES must have 24 passive tracers. Change jp_pisces in par_pisces.F90' )
+#endif
+
+   END SUBROUTINE trc_ctl_pisces
+  
 #else
    !!----------------------------------------------------------------------
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trclsm_pisces.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trclsm_pisces.F90	(revision 2524)
+++ 	(revision )
@@ -1,207 +1,0 @@
-MODULE trclsm_pisces
-   !!======================================================================
-   !!                      ***  MODULE trclsm_lobster  ***
-   !! TOP :   initialisation of some run parameters for PISCES bio-model
-   !!======================================================================
-   !! History :    -   !  1999-10 (M.A. Foujols, M. Levy) original code
-   !!              -   !  2000-01 (L. Bopp) hamocc3, p3zd
-   !!             1.0  !  2003-08 (C. Ethe)  module F90
-   !!             2.0  !  2007-12  (C. Ethe, G. Madec) from trclsm.pisces.h90
-   !!----------------------------------------------------------------------
-#if defined key_pisces
-   !!----------------------------------------------------------------------
-   !!   'key_pisces'   :                                   PISCES bio-model
-   !!----------------------------------------------------------------------
-   !! trc_lsm_pisces       : PISCES model namelist read
-   !!----------------------------------------------------------------------
-   USE oce_trc         ! Ocean variables
-   USE par_trc         ! TOP parameters
-   USE trc             ! TOP variables
-   USE sms_pisces      ! sms trends
-   USE in_out_manager  ! I/O manager
-
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   trc_lsm_pisces   ! called by trclsm.F90 module
-
-
-   !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
-   !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
-CONTAINS
-
-   SUBROUTINE trc_lsm_pisces
-      !!----------------------------------------------------------------------
-      !!                     ***  trc_lsm_pisces  ***  
-      !!
-      !! ** Purpose :   read PISCES namelist
-      !!
-      !! ** input   :   file 'namelist.trc.sms' containing the following
-      !!             namelist: natext, natbio, natsms
-      !!                       natkriest ("key_kriest")
-      !!----------------------------------------------------------------------
-      !!
-#if defined key_trc_diaadd && ! defined key_iomput
-      INTEGER ::  jl, jn
-      ! definition of additional diagnostic as a structure
-      TYPE DIAG
-         CHARACTER(len = 20)  :: snamedia   !: short name
-         CHARACTER(len = 80 ) :: lnamedia   !: long name
-         CHARACTER(len = 20 ) :: unitdia    !: unit
-      END TYPE DIAG
-
-      TYPE(DIAG) , DIMENSION(jp_pisces_2d) :: pisdia2d
-      TYPE(DIAG) , DIMENSION(jp_pisces_3d) :: pisdia3d
-#endif
-
-      NAMELIST/nampisbio/ part, nrdttrc, wsbio, xkmort, ferat3, wsbio2
-#if defined key_kriest
-      NAMELIST/nampiskrp/ xkr_eta, xkr_zeta, xkr_mass_min, xkr_mass_max
-#endif
-#if defined key_trc_diaadd && ! defined key_iomput
-      NAMELIST/nampisdia/ nwritedia, pisdia3d, pisdia2d     ! additional diagnostics
-#endif
-      NAMELIST/nampisdmp/ ln_pisdmp, ln_pisclo
-
-      !!----------------------------------------------------------------------
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' trc_lsm_pisces : read PISCES namelists'
-      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
-
-
-      !                               ! Open the namelist file
-      !                               ! ----------------------
-      CALL ctl_opn( numnat, 'namelist_pisces', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
-
-      REWIND( numnat )                    
-      READ  ( numnat, nampisbio )
-
-      IF(lwp) THEN                         ! control print
-         WRITE(numout,*) ' Namelist : nampisbio'
-         WRITE(numout,*) '    part of calcite not dissolved in guts     part      =', part
-         WRITE(numout,*) '    frequence pour la biologie                nrdttrc   =', nrdttrc
-         WRITE(numout,*) '    POC sinking speed                         wsbio     =', wsbio
-         WRITE(numout,*) '    half saturation constant for mortality    xkmort    =', xkmort
-         WRITE(numout,*) '    Fe/C in zooplankton                       ferat3    =', ferat3
-         WRITE(numout,*) '    Big particles sinking speed               wsbio2    =', wsbio2
-      ENDIF
-
-#if defined key_kriest
-
-      !                               ! nampiskrp : kriest parameters
-      !                               ! -----------------------------
-      xkr_eta      = 0.62        
-      xkr_zeta     = 1.62        
-      xkr_mass_min = 0.0002     
-      xkr_mass_max = 1.      
-
-      REWIND( numnat )                     ! read natkriest
-      READ  ( numnat, nampiskrp )
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) ' Namelist : nampiskrp'
-         WRITE(numout,*) '    Sinking  exponent                        xkr_eta      = ', xkr_eta
-         WRITE(numout,*) '    N content exponent                       xkr_zeta     = ', xkr_zeta
-         WRITE(numout,*) '    Minimum mass for Aggregates              xkr_mass_min = ', xkr_mass_min
-         WRITE(numout,*) '    Maximum mass for Aggregates              xkr_mass_max = ', xkr_mass_max
-         WRITE(numout,*)
-     ENDIF
-
-
-     ! Computation of some variables
-     xkr_massp = 5.7E-6 * 7.6 * xkr_mass_min**xkr_zeta
-
-#endif
-      !
-#if defined key_trc_diaadd && ! defined key_iomput
-
-      ! Namelist namlobdia
-      ! -------------------
-      nwritedia = 10                   ! default values
-
-      DO jl = 1, jp_pisces_2d
-         jn = jp_pcs0_2d + jl - 1
-         WRITE(ctrc2d(jn),'("2D_",I1)') jn                      ! short name
-         WRITE(ctrc2l(jn),'("2D DIAGNOSTIC NUMBER ",I2)') jn    ! long name
-         ctrc2u(jn) = ' '                                       ! units
-      END DO
-      !                                 ! 3D output arrays
-      DO jl = 1, jp_pisces_3d
-         jn = jp_pcs0_3d + jl - 1
-         WRITE(ctrc3d(jn),'("3D_",I1)') jn                      ! short name
-         WRITE(ctrc3l(jn),'("3D DIAGNOSTIC NUMBER ",I2)') jn    ! long name
-         ctrc3u(jn) = ' '                                       ! units
-      END DO
-
-      REWIND( numnat )               ! read natrtd
-      READ  ( numnat, nampisdia )
-
-      DO jl = 1, jp_pisces_2d
-         jn = jp_pcs0_2d + jl - 1
-         ctrc2d(jn) = pisdia2d(jl)%snamedia
-         ctrc2l(jn) = pisdia2d(jl)%lnamedia
-         ctrc2u(jn) = pisdia2d(jl)%unitdia
-      END DO
-
-      DO jl = 1, jp_pisces_3d
-         jn = jp_pcs0_3d + jl - 1
-         ctrc3d(jn) = pisdia3d(jl)%snamedia
-         ctrc3l(jn) = pisdia3d(jl)%lnamedia
-         ctrc3u(jn) = pisdia3d(jl)%unitdia
-      END DO
-
-      IF(lwp) THEN                   ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) ' Namelist : natadd'
-         WRITE(numout,*) '    frequency of outputs for additional arrays nwritedia = ', nwritedia
-         DO jl = 1, jp_pisces_3d
-            jn = jp_pcs0_3d + jl - 1
-            WRITE(numout,*) '   3d output field No : ',jn
-            WRITE(numout,*) '   short name         : ', TRIM(ctrc3d(jn))
-            WRITE(numout,*) '   long name          : ', TRIM(ctrc3l(jn))
-            WRITE(numout,*) '   unit               : ', TRIM(ctrc3u(jn))
-            WRITE(numout,*) ' '
-         END DO
-
-         DO jl = 1, jp_pisces_2d
-            jn = jp_pcs0_2d + jl - 1
-            WRITE(numout,*) '   2d output field No : ',jn
-            WRITE(numout,*) '   short name         : ', TRIM(ctrc2d(jn))
-            WRITE(numout,*) '   long name          : ', TRIM(ctrc2l(jn))
-            WRITE(numout,*) '   unit               : ', TRIM(ctrc2u(jn))
-            WRITE(numout,*) ' '
-         END DO
-      ENDIF
-#endif
-
-      REWIND( numnat )
-      READ  ( numnat, nampisdmp )
-
-      IF(lwp) THEN                         ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) ' Namelist : nampisdmp'
-         WRITE(numout,*) '    Relaxation of tracer to glodap mean value            ln_pisdmp      =', ln_pisdmp
-         WRITE(numout,*) '    Restoring of tracer to initial value  on closed seas  ln_pisclo      =', ln_pisclo
-         WRITE(numout,*) ' '
-      ENDIF
-
-   END SUBROUTINE trc_lsm_pisces
-
-#else
-   !!----------------------------------------------------------------------
-   !!  Dummy module :                                   No PISCES bio-model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE trc_lsm_pisces                      ! Empty routine
-   END  SUBROUTINE  trc_lsm_pisces
-#endif  
-
-   !!======================================================================
-END MODULE trclsm_pisces
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcnam_pisces.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcnam_pisces.F90	(revision 2528)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcnam_pisces.F90	(revision 2528)
@@ -0,0 +1,207 @@
+MODULE trcnam_pisces
+   !!======================================================================
+   !!                      ***  MODULE trcnam_lobster  ***
+   !! TOP :   initialisation of some run parameters for PISCES bio-model
+   !!======================================================================
+   !! History :    -   !  1999-10 (M.A. Foujols, M. Levy) original code
+   !!              -   !  2000-01 (L. Bopp) hamocc3, p3zd
+   !!             1.0  !  2003-08 (C. Ethe)  module F90
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec) from trcnam.pisces.h90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'   :                                   PISCES bio-model
+   !!----------------------------------------------------------------------
+   !! trc_nam_pisces       : PISCES model namelist read
+   !!----------------------------------------------------------------------
+   USE oce_trc         ! Ocean variables
+   USE par_trc         ! TOP parameters
+   USE trc             ! TOP variables
+   USE sms_pisces      ! sms trends
+   USE in_out_manager  ! I/O manager
+
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   trc_nam_pisces   ! called by trcnam.F90 module
+
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Id$ 
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE trc_nam_pisces
+      !!----------------------------------------------------------------------
+      !!                     ***  trc_nam_pisces  ***  
+      !!
+      !! ** Purpose :   read PISCES namelist
+      !!
+      !! ** input   :   file 'namelist.trc.sms' containing the following
+      !!             namelist: natext, natbio, natsms
+      !!                       natkriest ("key_kriest")
+      !!----------------------------------------------------------------------
+      !!
+#if defined key_diatrc && ! defined key_iomput
+      INTEGER ::  jl, jn
+      ! definition of additional diagnostic as a structure
+      TYPE DIAG
+         CHARACTER(len = 20)  :: snamedia   !: short name
+         CHARACTER(len = 80 ) :: lnamedia   !: long name
+         CHARACTER(len = 20 ) :: unitdia    !: unit
+      END TYPE DIAG
+
+      TYPE(DIAG) , DIMENSION(jp_pisces_2d) :: pisdia2d
+      TYPE(DIAG) , DIMENSION(jp_pisces_3d) :: pisdia3d
+#endif
+
+      NAMELIST/nampisbio/ part, nrdttrc, wsbio, xkmort, ferat3, wsbio2
+#if defined key_kriest
+      NAMELIST/nampiskrp/ xkr_eta, xkr_zeta, xkr_mass_min, xkr_mass_max
+#endif
+#if defined key_diatrc && ! defined key_iomput
+      NAMELIST/nampisdia/ nwritedia, pisdia3d, pisdia2d     ! additional diagnostics
+#endif
+      NAMELIST/nampisdmp/ ln_pisdmp, ln_pisclo
+
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' trc_nam_pisces : read PISCES namelists'
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+
+
+      !                               ! Open the namelist file
+      !                               ! ----------------------
+      CALL ctl_opn( numnat, 'namelist_pisces', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
+
+      REWIND( numnat )                    
+      READ  ( numnat, nampisbio )
+
+      IF(lwp) THEN                         ! control print
+         WRITE(numout,*) ' Namelist : nampisbio'
+         WRITE(numout,*) '    part of calcite not dissolved in guts     part      =', part
+         WRITE(numout,*) '    frequence pour la biologie                nrdttrc   =', nrdttrc
+         WRITE(numout,*) '    POC sinking speed                         wsbio     =', wsbio
+         WRITE(numout,*) '    half saturation constant for mortality    xkmort    =', xkmort
+         WRITE(numout,*) '    Fe/C in zooplankton                       ferat3    =', ferat3
+         WRITE(numout,*) '    Big particles sinking speed               wsbio2    =', wsbio2
+      ENDIF
+
+#if defined key_kriest
+
+      !                               ! nampiskrp : kriest parameters
+      !                               ! -----------------------------
+      xkr_eta      = 0.62        
+      xkr_zeta     = 1.62        
+      xkr_mass_min = 0.0002     
+      xkr_mass_max = 1.      
+
+      REWIND( numnat )                     ! read natkriest
+      READ  ( numnat, nampiskrp )
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) ' Namelist : nampiskrp'
+         WRITE(numout,*) '    Sinking  exponent                        xkr_eta      = ', xkr_eta
+         WRITE(numout,*) '    N content exponent                       xkr_zeta     = ', xkr_zeta
+         WRITE(numout,*) '    Minimum mass for Aggregates              xkr_mass_min = ', xkr_mass_min
+         WRITE(numout,*) '    Maximum mass for Aggregates              xkr_mass_max = ', xkr_mass_max
+         WRITE(numout,*)
+     ENDIF
+
+
+     ! Computation of some variables
+     xkr_massp = 5.7E-6 * 7.6 * xkr_mass_min**xkr_zeta
+
+#endif
+      !
+#if defined key_diatrc && ! defined key_iomput
+
+      ! Namelist namlobdia
+      ! -------------------
+      nwritedia = 10                   ! default values
+
+      DO jl = 1, jp_pisces_2d
+         jn = jp_pcs0_2d + jl - 1
+         WRITE(ctrc2d(jn),'("2D_",I1)') jn                      ! short name
+         WRITE(ctrc2l(jn),'("2D DIAGNOSTIC NUMBER ",I2)') jn    ! long name
+         ctrc2u(jn) = ' '                                       ! units
+      END DO
+      !                                 ! 3D output arrays
+      DO jl = 1, jp_pisces_3d
+         jn = jp_pcs0_3d + jl - 1
+         WRITE(ctrc3d(jn),'("3D_",I1)') jn                      ! short name
+         WRITE(ctrc3l(jn),'("3D DIAGNOSTIC NUMBER ",I2)') jn    ! long name
+         ctrc3u(jn) = ' '                                       ! units
+      END DO
+
+      REWIND( numnat )               ! read natrtd
+      READ  ( numnat, nampisdia )
+
+      DO jl = 1, jp_pisces_2d
+         jn = jp_pcs0_2d + jl - 1
+         ctrc2d(jn) = pisdia2d(jl)%snamedia
+         ctrc2l(jn) = pisdia2d(jl)%lnamedia
+         ctrc2u(jn) = pisdia2d(jl)%unitdia
+      END DO
+
+      DO jl = 1, jp_pisces_3d
+         jn = jp_pcs0_3d + jl - 1
+         ctrc3d(jn) = pisdia3d(jl)%snamedia
+         ctrc3l(jn) = pisdia3d(jl)%lnamedia
+         ctrc3u(jn) = pisdia3d(jl)%unitdia
+      END DO
+
+      IF(lwp) THEN                   ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) ' Namelist : natadd'
+         WRITE(numout,*) '    frequency of outputs for additional arrays nwritedia = ', nwritedia
+         DO jl = 1, jp_pisces_3d
+            jn = jp_pcs0_3d + jl - 1
+            WRITE(numout,*) '   3d output field No : ',jn
+            WRITE(numout,*) '   short name         : ', TRIM(ctrc3d(jn))
+            WRITE(numout,*) '   long name          : ', TRIM(ctrc3l(jn))
+            WRITE(numout,*) '   unit               : ', TRIM(ctrc3u(jn))
+            WRITE(numout,*) ' '
+         END DO
+
+         DO jl = 1, jp_pisces_2d
+            jn = jp_pcs0_2d + jl - 1
+            WRITE(numout,*) '   2d output field No : ',jn
+            WRITE(numout,*) '   short name         : ', TRIM(ctrc2d(jn))
+            WRITE(numout,*) '   long name          : ', TRIM(ctrc2l(jn))
+            WRITE(numout,*) '   unit               : ', TRIM(ctrc2u(jn))
+            WRITE(numout,*) ' '
+         END DO
+      ENDIF
+#endif
+
+      REWIND( numnat )
+      READ  ( numnat, nampisdmp )
+
+      IF(lwp) THEN                         ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) ' Namelist : nampisdmp'
+         WRITE(numout,*) '    Relaxation of tracer to glodap mean value            ln_pisdmp      =', ln_pisdmp
+         WRITE(numout,*) '    Restoring of tracer to initial value  on closed seas  ln_pisclo      =', ln_pisclo
+         WRITE(numout,*) ' '
+      ENDIF
+
+   END SUBROUTINE trc_nam_pisces
+
+#else
+   !!----------------------------------------------------------------------
+   !!  Dummy module :                                   No PISCES bio-model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_nam_pisces                      ! Empty routine
+   END  SUBROUTINE  trc_nam_pisces
+#endif  
+
+   !!======================================================================
+END MODULE trcnam_pisces
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcrst_pisces.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcrst_pisces.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcrst_pisces.F90	(revision 2528)
@@ -21,4 +21,6 @@
    USE iom
    USE trcdta
+   USE lib_mpp
+   USE lib_fortran
 
    IMPLICIT NONE
@@ -118,35 +120,20 @@
       IF(lwp)  WRITE(numout,*)
 
-      IF( cp_cfg == "orca" .AND. .NOT. lk_trc_c1d ) THEN      ! ORCA condiguration (not 1D) !
+      IF( cp_cfg == "orca" .AND. .NOT. lk_c1d ) THEN      ! ORCA condiguration (not 1D) !
          !                                                    ! --------------------------- !
          ! set total alkalinity, phosphate, nitrate & silicate
 
-         zalksum = 0.e0
-         zpo4sum = 0.e0
-         zno3sum = 0.e0
-         zsilsum = 0.e0
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zvol = cvol(ji,jj,jk)
-#  if defined key_off_degrad
-                  zvol = zvol * facvol(ji,jj,jk)
-#  endif
-                  zalksum = zalksum + trn(ji,jj,jk,jptal) * zvol
-                  zpo4sum = zpo4sum + trn(ji,jj,jk,jppo4) * zvol
-                  zno3sum = zno3sum + trn(ji,jj,jk,jpno3) * zvol
-                  zsilsum = zsilsum + trn(ji,jj,jk,jpsil) * zvol
-               END DO
-            END DO
-         END DO
-         IF( lk_mpp )   CALL mpp_sum( zalksum )     ! sum over the global domain
-         IF( lk_mpp )   CALL mpp_sum( zpo4sum )     ! sum over the global domain
-         IF( lk_mpp )   CALL mpp_sum( zno3sum )     ! sum over the global domain
-         IF( lk_mpp )   CALL mpp_sum( zsilsum )     ! sum over the global domain
          zarea   = 1. / areatot * 1.e6
-         zalksum = zalksum * zarea
-         zpo4sum = zpo4sum * zarea / 122.
-         zno3sum = zno3sum * zarea / 7.6
-         zsilsum = zsilsum * zarea
+# if defined key_degrad
+         zalksum = glob_sum( trn(:,:,:,jptal) * cvol(:,:,:) * facvol(:,:,:) ) * zarea
+         zpo4sum = glob_sum( trn(:,:,:,jppo4) * cvol(:,:,:) * facvol(:,:,:) ) * zarea / 122.
+         zno3sum = glob_sum( trn(:,:,:,jpno3) * cvol(:,:,:) * facvol(:,:,:) ) * zarea / 7.6
+         zsilsum = glob_sum( trn(:,:,:,jpsil) * cvol(:,:,:) * facvol(:,:,:) ) * zarea
+# else
+         zalksum = glob_sum( trn(:,:,:,jptal) * cvol(:,:,:)  ) * zarea
+         zpo4sum = glob_sum( trn(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea / 122.
+         zno3sum = glob_sum( trn(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea / 7.6
+         zsilsum = glob_sum( trn(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
+# endif
 
          IF(lwp) WRITE(numout,*) '       TALK mean : ', zalksum
@@ -263,5 +250,5 @@
 #if defined key_dtatrc
       ! Restore close seas values to initial data
-      CALL trc_dta( nittrc000 ) 
+      CALL trc_dta( nit000 ) 
       DO jn = 1, jptra
          IF( lutini(jn) ) THEN
Index: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcsms_pisces.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcsms_pisces.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/trcsms_pisces.F90	(revision 2528)
@@ -22,10 +22,20 @@
    USE p4zche          ! 
    USE p4zbio          ! 
+   USE p4zsink         ! 
+   USE p4zopt          ! 
+   USE p4zlim          ! 
+   USE p4zprod         !
+   USE p4zmort         !
+   USE p4zmicro        ! 
+   USE p4zmeso         ! 
+   USE p4zrem          ! 
    USE p4zsed          ! 
    USE p4zlys          ! 
    USE p4zflx          ! 
 
-   USE trdmld_trc_oce
-   USE trdmld_trc
+   USE prtctl_trc
+
+   USE trdmod_oce
+   USE trdmod_trc
 
    USE sedmodel
@@ -37,7 +47,7 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -59,11 +69,16 @@
       INTEGER ::   jnt, jn
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztrpis   ! used for pisces sms trends
+      CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
 
-      IF( kt == nittrc000  .AND. .NOT. ln_rsttr )   CALL trc_sms_pisces_init    ! Initialization (first time-step only)
+      IF( kt == nit000 )   CALL trc_sms_pisces_init    ! Initialization (first time-step only)
 
-      IF( ndayflxtr /= nday ) THEN      ! New days
+      IF( ndayflxtr /= nday_year ) THEN      ! New days
          !
-         ndayflxtr = nday
+         ndayflxtr = nday_year
+
+         IF(lwp) write(numout,*)
+         IF(lwp) write(numout,*) ' New chemical constants and various rates for biogeochemistry at new day : ', nday_year
+         IF(lwp) write(numout,*) '~~~~~~'
 
          CALL p4z_che          ! computation of chemical constants
@@ -71,5 +86,4 @@
          !
       ENDIF
-
 
       DO jnt = 1, nrdttrc          ! Potential time splitting if requested
@@ -91,8 +105,9 @@
       END DO
 
+
       IF( l_trdtrc ) THEN
           DO jn = jp_pcs0, jp_pcs1
             ztrpis(:,:,:) = tra(:,:,:,jn)
-            CALL trd_mod_trc( ztrpis, jn, jptrc_trd_sms, kt )   ! save trends
+            CALL trd_mod_trc( ztrpis, jn, jptra_trd_sms, kt )   ! save trends
           END DO
       END IF
@@ -121,18 +136,38 @@
       REAL(wp) ::  ztmas, ztmas1
 
-      ! Initialization of chemical variables of the carbon cycle
-      ! --------------------------------------------------------
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ztmas   = tmask(ji,jj,jk)
-               ztmas1  = 1. - tmask(ji,jj,jk)
-               zcaralk = trn(ji,jj,jk,jptal) - borat(ji,jj,jk) / (  1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) )  )
-               zco3    = ( zcaralk - trn(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1
-               zbicarb = ( 2. * trn(ji,jj,jk,jpdic) - zcaralk )
-               hi(ji,jj,jk) = ( ak23(ji,jj,jk) * zbicarb / zco3 ) * ztmas + 1.e-9 * ztmas1
+      IF( .NOT. ln_rsttr ) THEN
+         ! Initialization of chemical variables of the carbon cycle
+         ! --------------------------------------------------------
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztmas   = tmask(ji,jj,jk)
+                  ztmas1  = 1. - tmask(ji,jj,jk)
+                  zcaralk = trn(ji,jj,jk,jptal) - borat(ji,jj,jk) / (  1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) )  )
+                  zco3    = ( zcaralk - trn(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1
+                  zbicarb = ( 2. * trn(ji,jj,jk,jpdic) - zcaralk )
+                  hi(ji,jj,jk) = ( ak23(ji,jj,jk) * zbicarb / zco3 ) * ztmas + 1.e-9 * ztmas1
+               END DO
             END DO
          END DO
-      END DO
+         !
+      END IF
+
+      ! Time step duration for biology
+      xstep = rfact2 / rday
+
+      CALL p4z_sink_init      ! vertical flux of particulate organic matter
+      CALL p4z_opt_init       ! Optic: PAR in the water column
+      CALL p4z_lim_init       ! co-limitations by the various nutrients
+      CALL p4z_prod_init      ! phytoplankton growth rate over the global ocean. 
+      CALL p4z_rem_init       ! remineralisation
+      CALL p4z_mort_init      ! phytoplankton mortality
+      CALL p4z_micro_init     ! microzooplankton
+      CALL p4z_meso_init      ! mesozooplankton
+      CALL p4z_sed_init       ! sedimentation
+      CALL p4z_lys_init       ! calcite saturation
+      CALL p4z_flx_init       ! gas exchange
+
+      ndayflxtr = 0
 
    END SUBROUTINE trc_sms_pisces_init
