Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zbio.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zbio.F90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zbio.F90	(revision 858)
@@ -27,4 +27,5 @@
    USE p4zmeso         ! 
    USE p4zrem          ! 
+   USE prtctl_trc
    
    IMPLICIT NONE
@@ -43,5 +44,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_bio
+   SUBROUTINE p4z_bio ( kt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_bio  ***
@@ -53,14 +54,8 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kt
       INTEGER  ::   ji, jj, jk, jn
-      REAL(wp) ::   zdenom, ztemp
-      REAL(wp) ::   zprodt, zprodca
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdenom1
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zfracal
-#if defined key_kriest
-      REAL(wp) ::   znumpoc, znumdoc
-#else
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdenom2
-#endif
+      CHARACTER (len=25) :: charout
+
       !!---------------------------------------------------------------------
 
@@ -68,5 +63,5 @@
       !     OF PHYTOPLANKTON AND DETRITUS
 
-      zdiss(:,:,:) = 0.01
+      xdiss(:,:,:) = 0.01
 
 !!gm the use of nmld should be better here?
@@ -74,44 +69,9 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( fsdepw(ji,jj,jk+1) .le. hmld(ji,jj) )   zdiss(ji,jj,jk) = 1.e0
+               IF( fsdepw(ji,jj,jk+1) .le. hmld(ji,jj) )   xdiss(ji,jj,jk) = 1.e0
             END DO 
          END DO
       END DO
 
-      ! Compute de different ratios for scavenging of iron
-      ! --------------------------------------------------
-
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-#if ! defined key_kriest
-               zdenom = 1. / (  trn(ji,jj,jk,jppoc) + trn(ji,jj,jk,jpgoc)           &
-                  &           + trn(ji,jj,jk,jpdsi) + trn(ji,jj,jk,jpcal) + rtrn )
-               zdenom1(ji,jj,jk) = trn(ji,jj,jk,jppoc) * zdenom
-               zdenom2(ji,jj,jk) = trn(ji,jj,jk,jpgoc) * zdenom
-#else
-               zdenom = 1. / (  trn(ji,jj,jk,jppoc)                                 &
-                  &           + trn(ji,jj,jk,jpdsi) + trn(ji,jj,jk,jpcal) + rtrn )
-               zdenom1(ji,jj,jk) = trn(ji,jj,jk,jppoc) * zdenom
-#endif
-            END DO
-         END DO
-      END DO
-
-      ! Compute the fraction of nanophytoplankton that is made of calcifiers
-      ! --------------------------------------------------------------------
-
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ztemp = MAX( 0., tn(ji,jj,jk) )
-               zfracal(ji,jj,jk) = caco3r * xlimphy(ji,jj,jk)   &
-                  &                       * MAX( 0.0001, ztemp / ( 2.+ ztemp ) )   &
-                  &                       * MAX( 1., trn(ji,jj,jk,jpphy) * 1.e6 / 2. )
-               zfracal(ji,jj,jk) = MIN( 0.8 , zfracal(ji,jj,jk) )
-               zfracal(ji,jj,jk) = MAX( 0.01, zfracal(ji,jj,jk) )
-            END DO
-         END DO
-      END DO
 
       ! computation of the vertical flux of particulate organic matter
@@ -122,8 +82,11 @@
       ENDIF
 
+
       ! compute the PAR in the water column
       ! -----------------------------------
 
       CALL p4z_opt       ! Optical 
+
+
 
       ! compute the co-limitations by the various nutrients
@@ -132,9 +95,12 @@
       CALL p4z_lim
 
+
       ! compute phytoplankton growth rate over the global ocean. 
       ! -------------------------------------------------------
       ! (Growth rates for each element is computed (C, Si, Fe, Chl))
 
-      CALL p4z_prod
+      CALL p4z_prod ( kt )
+
+
 
       ! phytoplankton mortality (Mortality losses for each elements are computed (C, Fe, Si, Chl) )
@@ -143,5 +109,7 @@
       CALL p4z_nano       ! nanophytoplankton
       
+
       CALL p4z_diat       ! diatoms
+
 
       ! zooplankton sources/sinks routines (each elements are computed (C, Fe, Si, Chl) )
@@ -149,6 +117,8 @@
 
       CALL p4z_micro      ! microzooplankton
+
       
       CALL p4z_meso       ! mesozooplankton
+
 
       ! computation of remineralization terms of organic matter + scavenging of Fe
@@ -158,154 +128,14 @@
 
 
-      ! Determination of tracers concentration as a function of biological sources and sinks
-      ! ------------------------------------------------------------------------------------
-
-      DO jk = 1, jpkm1
-
-         ! Evolution of PO4
-         ! ----------------    
-         trn(:,:,jk,jppo4) = trn(:,:,jk,jppo4) - prorca(:,:,jk) - prorca2(:,:,jk)            &
-            &                                  + olimi (:,:,jk) + grarem (:,:,jk) * sigma1   &
-            &                                  + denitr(:,:,jk) + grarem2(:,:,jk) * sigma2
-
-         ! Evolution of NO3 and NH4
-         ! ------------------------
-         trn(:,:,jk,jpno3) = trn(:,:,jk,jpno3) - pronew(:,:,jk) - pronew2(:,:,jk)            &
-            &                                  + onitr (:,:,jk) - denitr (:,:,jk) * rdenit
-
-         trn(:,:,jk,jpnh4) = trn(:,:,jk,jpnh4) - proreg(:,:,jk) - proreg2(:,:,jk)            &
-            &                                  + olimi (:,:,jk) + grarem (:,:,jk) * sigma1   &
-            &                                                   + grarem2(:,:,jk) * sigma2   &
-            &                                  - onitr (:,:,jk) + denitr (:,:,jk)
-
-         ! Evolution of Phytoplankton
-         ! --------------------------
-         trn(:,:,jk,jpphy) = trn(:,:,jk,jpphy) + prorca (:,:,jk) * ( 1.- excret ) - tortp(:,:,jk)   &
-            &                                  - grazp  (:,:,jk) - grazn(:,:,jk)  - respp(:,:,jk)
-
-         trn(:,:,jk,jpnch) = trn(:,:,jk,jpnch) + prorca6(:,:,jk) * ( 1.- excret ) - tortnch(:,:,jk)   &
-            &                                  - grazpch(:,:,jk) - graznch(:,:,jk)- respnch(:,:,jk)
-
-         ! Evolution of Diatoms
-         ! --------------------
-         trn(:,:,jk,jpdia) = trn(:,:,jk,jpdia) + prorca2(:,:,jk) * ( 1.- excret2 ) - tortp2(:,:,jk)   &
-            &                                  - respp2 (:,:,jk) - grazd(:,:,jk)   - grazsd(:,:,jk)
-
-         trn(:,:,jk,jpdch) = trn(:,:,jk,jpdch) + prorca7(:,:,jk) * ( 1.- excret2 ) - tortdch(:,:,jk)   &
-            &                                  - respdch(:,:,jk) - grazdch(:,:,jk) - grazsch(:,:,jk)
-
-         ! Evolution of Zooplankton
-         ! ------------------------   
-         trn(:,:,jk,jpzoo) = trn(:,:,jk,jpzoo) + epsher * ( grazp(:,:,jk) + grazm(:,:,jk) + grazsd(:,:,jk) )   &
-            &                                  - grazz(:,:,jk) - tortz(:,:,jk) - respz(:,:,jk)
-
-         ! Evolution of Mesozooplankton
-         ! ------------------------
-         trn(:,:,jk,jpmes) = trn(:,:,jk,jpmes) + epsher2 * (  grazd  (:,:,jk) + grazz  (:,:,jk) + grazn(:,:,jk)   &
-            &                                               + grazpoc(:,:,jk) + grazffe(:,:,jk)  )                &
-            &                                  - tortz2(:,:,jk) - respz2(:,:,jk)
-    
-         ! Evolution of O2
-         ! ---------------     
-         trn(:,:,jk,jpoxy) = trn(:,:,jk,jpoxy) + o2ut * (  proreg(:,:,jk) + proreg2(:,:,jk) - olimi(:,:,jk)        &
-            &                                             -grarem(:,:,jk) * sigma1 - grarem2(:,:,jk) * sigma2  )   &
-            &                                  + ( o2ut + o2nit ) * ( pronew(:,:,jk) + pronew2(:,:,jk) )           &
-            &                                  - o2nit * onitr(:,:,jk)
-
-         ! Evolution of IRON
-         ! -----------------
-         trn(:,:,jk,jpfer) = trn(:,:,jk,jpfer) + ( excret  - 1.) * prorca5(:,:,jk) - xaggdfe (:,:,jk)   &
-            &                                  + ( excret2 - 1.) * prorca4(:,:,jk) - xbactfer(:,:,jk)   &
-            &                                  + grafer(:,:,jk) + grafer2(:,:,jk)                       &
-            &                                  + ofer  (:,:,jk) - xscave (:,:,jk)
-         !
-      END DO
-
-
-#if defined key_kriest
-
-#include "p4zbio_kriest.h90"
-
-#else
-
-#include "p4zbio_std.h90"
-
-#endif
-
-
-      DO jk = 1, jpkm1
-
-         ! Evolution of biogenic Silica
-         ! ----------------------------
-         trn(:,:,jk,jpbsi) = trn(:,:,jk,jpbsi) + prorca3(:,:,jk) * ( 1.- excret2 ) - grazss(:,:,jk)   &
-            &                                  - tortds (:,:,jk) - respds(:,:,jk)  - grazs (:,:,jk)
-
-         ! Evolution of sinking biogenic silica
-         ! ------------------------------------
-         trn(:,:,jk,jpdsi) = trn(:,:,jk,jpdsi) + tortds (:,:,jk) + respds(:,:,jk)  + grazs(:,:,jk)   &
-            &                                  - osil   (:,:,jk) + grazss(:,:,jk)
-
-         ! Evolution of biogenic diatom Iron
-         ! ---------------------------------
-         trn(:,:,jk,jpdfe) = trn(:,:,jk,jpdfe) + prorca4(:,:,jk) * ( 1.- excret2 ) - grazsf(:,:,jk)   &
-            &                                  - tortdf (:,:,jk) - respdf(:,:,jk)  - grazf (:,:,jk)
-
-         ! Evolution of biogenic nanophytoplankton Iron
-         ! --------------------------------------------
-         trn(:,:,jk,jpnfe) = trn(:,:,jk,jpnfe) + prorca5(:,:,jk) * ( 1.- excret )  - graznf(:,:,jk)   &
-            &                                  - tortnf (:,:,jk) - respnf(:,:,jk)  - grazpf(:,:,jk)
-
-         ! Evolution of dissolved Silica
-         ! -----------------------------
-         trn(:,:,jk,jpsil) = trn(:,:,jk,jpsil) - ( 1.- excret2 ) * prorca3(:,:,jk) + osil(:,:,jk)
-
-      END DO
-     
-      ! Evolution of calcite and silicates as a function of the two tracers
-      ! -------------------------------------------------------------------     
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-
-               zprodt  = prorca(ji,jj,jk) + prorca2(ji,jj,jk) - olimi(ji,jj,jk) - grarem(ji,jj,jk) * sigma1   &
-                  &                       - grarem2(ji,jj,jk) * sigma2 - denitr(ji,jj,jk)
-               zprodca = pronew(ji,jj,jk) + pronew2(ji,jj,jk) - onitr(ji,jj,jk) + rdenit * denitr(ji,jj,jk)
-    
-               ! potential production of calcite and biogenic silicate
-               ! ------------------------------------------------------     
-               prcaca(ji,jj,jk) = zfracal(ji,jj,jk)   &
-                  &             * (  part * (  unass*grazp(ji,jj,jk) + unass2*grazn(ji,jj,jk)  )   &
-                  &                + tortp(ji,jj,jk) + respp(ji,jj,jk)  )
-
-               ! Consumption of Total (12C)O2
-               ! ----------------------------     
-               trn(ji,jj,jk,jpdic) = trn(ji,jj,jk,jpdic) - zprodt - prcaca(ji,jj,jk)
-
-               ! Consumption of alkalinity due to ca++ uptake and increase of
-               !  alkalinity due to nitrate consumption during organic soft tissue production
-               ! ---------------------------------------------------------  
-               trn(ji,jj,jk,jptal) = trn(ji,jj,jk,jptal) + rno3 * zprodca - 2.* prcaca(ji,jj,jk)
-               !
-            END DO
-         END DO
-      END DO
-
-
-      ! Production of calcite due to biological production
-      ! --------------------------------------------------   
-      DO jk = 1, jpkm1
-         trn(:,:,jk,jpcal) = trn(:,:,jk,jpcal) + prcaca(:,:,jk)
-      END DO
-
-
       ! Loop to test if tracers concentrations fall below 0.
       ! ----------------------------------------------------
 
-      znegtr(:,:,:) = 1.e0
-      DO jn = 1, jptra
+      xnegtr(:,:,:) = 1.e0
+      DO jn = jp_pcs0, jp_pcs1
          DO jk = 1, jpk
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  IF( trn(ji,jj,jk,jn) < 0.e0 )   znegtr(ji,jj,jk) = 0.e0
+                  IF( ( trn(ji,jj,jk,jn) + tra(ji,jj,jk,jn) ) < 0.e0 )   &
+             &      xnegtr(ji,jj,jk) = 0.e0
                END DO
             END DO
@@ -314,24 +144,38 @@
       !                                ! where at least 1 tracer concentration becomes negative
       !                                ! all tracer tendancy are set to zero (i.e. trn = trb)
-      DO jn = 1, jptra
-         trn(:,:,:,jn) = trb(:,:,:,jn) + znegtr(:,:,:) * ( trn(:,:,:,jn) - trb(:,:,:,jn) )
+      DO jn = jp_pcs0, jp_pcs1
+         trn(:,:,:,jn) = trn(:,:,:,jn) + xnegtr(:,:,:) * tra(:,:,:,jn)
       END DO
 
-# if defined key_trc_dia3d
+
+      tra(:,:,:,:) = 0.0
+
+#if defined key_kriest
+        DO jk = 1,jpkm1
+           DO jj = 1,jpj
+              DO ji = 1,jpi
+                 trn(ji,jj,jk,jpnum) = MAX( trn(ji,jj,jk,jpnum),  &
+     &                                      trn(ji,jj,jk,jppoc) / xkr_massp / xnumm(jk) )
+
+                 trn(ji,jj,jk,jpnum) = MIN( trn(ji,jj,jk,jpnum),  &
+     &                                      trn(ji,jj,jk,jppoc) / xkr_massp / 1.1 )
+
+          END DO
+        END DO
+      END DO
+#endif
+
+
+# if defined key_trc_dia3d && defined key_kriest
 !!gm potential bug  hard coded index on trc3d
-          trc3d(:,:,:, 4) = etot(:,:,:)
-          trc3d(:,:,:, 5) = prorca (:,:,:) * znegtr(:,:,:) * 1.e3 * rfact2r
-          trc3d(:,:,:, 6) = prorca2(:,:,:) * znegtr(:,:,:) * 1.e3 * rfact2r
-          trc3d(:,:,:, 7) = pronew (:,:,:) * znegtr(:,:,:) * 1.e3 * rfact2r
-          trc3d(:,:,:, 8) = pronew2(:,:,:) * znegtr(:,:,:) * 1.e3 * rfact2r
-          trc3d(:,:,:, 9) = prorca3(:,:,:) * znegtr(:,:,:) * 1.e3 * rfact2r
-          trc3d(:,:,:,10) = prorca4(:,:,:) * znegtr(:,:,:) * 1.e3 * rfact2r
-#  if ! defined key_kriest
-          trc3d(:,:,:,11) = prorca5(:,:,:) * znegtr(:,:,:) * 1.e3 * rfact2r
-#  else
-          trc3d(:,:,:,11) = prcaca (:,:,:) * znegtr(:,:,:) * 1.e3 * rfact2r
-#  endif
+     trc3d(:,:,:,11) = tra(:,:,:,jpcal) * xnegtr(:,:,:) * 1.e3 * rfact2r
 # endif
       !
+       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+         WRITE(charout, FMT="('bio ')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc(tab4d=trn, mask=tmask, clinfo=ctrcnm)
+       ENDIF
+
    END SUBROUTINE p4z_bio
 
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zche.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zche.F90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zche.F90	(revision 858)
@@ -25,4 +25,155 @@
 
    PUBLIC   p4z_che    ! called in p4zprg.F90
+
+   !! * Module variables
+
+   REAL(wp) :: &
+      salchl = 1./1.80655 ! conversion factor for salinity --> chlorinity (Wooster et al. 1969)
+
+   REAL(wp) :: &            ! coeff. for apparent solubility equilibrium 
+      akcc1 = -171.9065 , &    ! Millero et al. 1995 from Mucci 1983
+      akcc2 = -0.077993 , &  
+      akcc3 = 2839.319  , &  
+      akcc4 = 71.595    , &  
+      akcc5 = -0.77712  , &  
+      akcc6 = 0.0028426 , &  
+      akcc7 = 178.34    , &  
+      akcc8 = -0.07711  , &  
+      akcc9 = 0.0041249
+
+
+   REAL(wp) :: &             ! universal gas constants
+      rgas = 83.143, &
+      oxyco = 1./22.4144
+
+   REAL(wp) :: &             ! borat constants
+      bor1 = 0.00023, &
+      bor2 = 1./10.82
+
+   REAL(wp) :: &              !
+      ca0 = -162.8301  , &
+      ca1 = 218.2968   , &
+      ca2 = 90.9241    , &
+      ca3 = -1.47696   , &
+      ca4 = 0.025695   , &
+      ca5 = -0.025225  , &
+      ca6 = 0.0049867
+
+   REAL(wp) :: &              ! coeff. for 1. dissoc. of carbonic acid (Edmond and Gieskes, 1970)   
+      c10 = -3670.7   , &
+      c11 = 62.008    , &
+      c12 = -9.7944   , &
+      c13 = 0.0118    , &
+      c14 = -0.000116
+
+   REAL(wp) :: &              ! coeff. for 2. dissoc. of carbonic acid (Millero, 1995)   
+      c20 = -1394.7   , &
+      c21 = -4.777    , &
+      c22 = 0.0184    , &
+      c23 = -0.000118
+
+   REAL(wp) :: &             ! constants for calculate concentrations 
+      st1  = 0.14     , &    ! for sulfate (Morris & Riley 1966)
+      st2  = 1./96.062, &
+      ks0  = 141.328  , &
+      ks1  = -4276.1  , &
+      ks2  = -23.093  , &
+      ks3  = -13856.  , &
+      ks4  = 324.57   , &
+      ks5  = -47.986  , &
+      ks6  = 35474.   , &
+      ks7  = -771.54  , &
+      ks8  = 114.723  , &
+      ks9  = -2698.   , &
+      ks10 = 1776.    , &
+      ks11 = 1.       , &
+      ks12 = -0.001005 
+
+   REAL(wp) :: &             ! constants for calculate concentrations 
+      ft1  = 0.000067   , &  ! fluorides (Dickson & Riley 1979 )
+      ft2  = 1./18.9984 , &
+      kf0  = -12.641    , &
+      kf1  = 1590.2     , &
+      kf2  = 1.525      , &
+      kf3  = 1.0        , &
+      kf4  =-0.001005
+
+
+   REAL(wp) :: &              ! coeff. for 1. dissoc. of boric acid (Dickson and Goyet, 1994)
+      cb0  = -8966.90, &
+      cb1  = -2890.53, &
+      cb2  = -77.942 , &
+      cb3  = 1.728   , &
+      cb4  = -0.0996 , &
+      cb5  = 148.0248, &
+      cb6  = 137.1942, &
+      cb7  = 1.62142 , &
+      cb8  = -24.4344, &
+      cb9  = -25.085 , &
+      cb10 = -0.2474 , &
+      cb11 = 0.053105
+
+   REAL(wp) :: &             ! coeff. for dissoc. of water (Dickson and Riley, 1979 )
+      cw0 = -13847.26  , &
+      cw1 = 148.9652   , &
+      cw2 = -23.6521   , &
+      cw3 = 118.67     , &
+      cw4 = -5.977     , &
+      cw5 = 1.0495     , &
+      cw6 = -0.01615
+ 
+   REAL(wp) :: &             ! coeff. for dissoc. of phosphate (Millero (1974)
+      cp10 = 115.54    , &
+      cp11 = -4576.752 , &
+      cp12 = -18.453   , &
+      cp13 = -106.736  , &
+      cp14 = 0.69171   , &
+      cp15 = -0.65643  , &
+      cp16 = -0.01844  , &
+      cp20 = 172.1033  , &
+      cp21 = -8814.715 , &
+      cp22 = -27.927   , &
+      cp23 = -160.340  , &
+      cp24 = 1.3566    , &
+      cp25 = 0.37335   , &
+      cp26 = -0.05778  , &
+      cp30 = -18.126   , &
+      cp31 = -3070.75  , &
+      cp32 = 17.27039  , &
+      cp33 = 2.81197   , &
+      cp34 = -44.99486 , &
+      cp35 = -0.09984
+
+   REAL(wp) :: &             ! coeff. for dissoc. of silicates (Millero (1974)  
+      cs10 = 117.385   , &  
+      cs11 = -8904.2   , & 
+      cs12 = -19.334   , & 
+      cs13 = -458.79   , & 
+      cs14 = 3.5913    , & 
+      cs15 = 188.74    , & 
+      cs16 = -1.5998   , & 
+      cs17 = -12.1652  , & 
+      cs18 = 0.07871   , & 
+      cs19 = 0.        , & 
+      cs20 = 1.        , & 
+      cs21 = -0.001005
+
+   REAL(wp) :: &              ! volumetric solubility constants for o2 in ml/l (Weiss, 1974)
+      ox0 = -58.3877   , &
+      ox1 = 85.8079    , &
+      ox2 = 23.8439    , &
+      ox3 = -0.034892  , &
+      ox4 = 0.015568   , &
+      ox5 = -0.0019387 
+
+   REAL(wp), DIMENSION(5)  :: &  ! coeff. for seawater pressure correction 
+      devk1, devk2, devk3,    &  ! (millero 95)
+      devk4, devk5
+
+   DATA devk1 / -25.5    , -15.82    , -29.48  , -25.60     , -48.76    /   
+   DATA devk2 / 0.1271   , -0.0219   , 0.1622  , 0.2324     , 0.5304    /   
+   DATA devk3 / 0.       , 0.        , 2.608E-3,  -3.6246E-3, 0.        /   
+   DATA devk4 / -3.08E-3 , 1.13E-3   , -2.84E-3, -5.13E-3   , -11.76E-3 /   
+   DATA devk5 / 0.0877E-3, -0.1475E-3,  0.     , 0.0794E-3  , 0.3692E-3 /
 
    !!* Substitution
@@ -46,16 +197,16 @@
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   ztkel, zsal , zqtt  , zbuf1 , zbuf2
-      REAL(wp) ::   zpres, ztc  , zcl   , zcpexp, zcek0, zoxy  , zcpexp2
+      REAL(wp) ::   zpres, ztc  , zcl   , zcpexp, zoxy  , zcpexp2
       REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1
       REAL(wp) ::   zlqtt, zqtt2, zsal15, zis   , zis2 , zisqrt
       REAL(wp) ::   zckb , zck1 , zck2  , zckw  , zak1 , zak2  , zakb , zaksp0, zakw
-      REAL(wp) ::   zckp1, zckp2, zckp3 , zcksi , zakp1, zakp2 , zakp3, zaksi
-      REAL(wp) ::   zst  , zft  , zcks  , zckf  , zaks , zakf  , zaksp1
+      REAL(wp) ::   zst  , zft  , zcks  , zckf  , zaksp1
       !!---------------------------------------------------------------------
 
       ! CHEMICAL CONSTANTS - SURFACE LAYER
       ! ----------------------------------
-
+!CDIR NOVERRCHK
       DO jj = 1, jpj
+!CDIR NOVERRCHK
          DO ji = 1, jpi
 
@@ -69,5 +220,4 @@
             !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980)
             !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS
-            zcek0 = c00 + c01 / zqtt + c02 * zlqtt + zsal * ( c03 + c04 * zqtt + c05 * zqtt2 )
             zcek1 = ca0 + ca1 / zqtt + ca2 * zlqtt + ca3 * zqtt2 + zsal*( ca4 + ca5 * zqtt + ca6 * zqtt2 )
 
@@ -76,7 +226,6 @@
 
             !                             ! SET SOLUBILITIES OF O2 AND CO2
-            chemc(ji,jj,1) = EXP( zcek0 ) * 1.e-6 * rhop(ji,jj,1) / 1000.
+            chemc(ji,jj,1) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000.
             chemc(ji,jj,2) = EXP( zoxy  ) * oxyco
-            chemc(ji,jj,3) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000.
 
          END DO
@@ -85,7 +234,9 @@
       ! CHEMICAL CONSTANTS - DEEP OCEAN
       ! -------------------------------
-
+!CDIR NOVERRCHK
       DO jk = 1, jpk
+!CDIR NOVERRCHK
          DO jj = 1, jpj
+!CDIR NOVERRCHK
             DO ji = 1, jpi
 
@@ -136,11 +287,4 @@
                zckw    = cw0 * ztr + cw1 + cw2 * zlogt + ( cw3 * ztr + cw4 + cw5 * zlogt ) * zsqrt + cw6 * zsal
 
-               ! DISSOCIATION CONSTANT FOR PHOSPHATE AND SILICATE (seawater scale)
-               zckp1 = cp10 + cp11 * ztr + cp12 * zlogt + zsqrt * ( cp13 * ztr + cp14 ) + zsal * ( cp15 * ztr + cp16 )
-               zckp2 = cp20 + cp21 * ztr + cp22 * zlogt + zsqrt * ( cp23 * ztr + cp24 ) + zsal * ( cp25 * ztr + cp26 )
-               zckp3 = cp30 + cp31 * ztr                + zsqrt * ( cp32 * ztr + cp33 ) + zsal * ( cp34 * ztr + cp35 )
-               zcksi = cs10 + cs11 * ztr + cs12 * zlogt + zisqrt* ( cs13 * ztr + cs14 ) + zis * ( cs15 * ztr + cs16 )  &
-                  &                                    + zis2  * ( cs17 * ztr + cs18 ) + LOG( 1.   + cs19 * zsal )     &
-                  &                                                                    + LOG( cs20 + cs21 * zsal )
 
                ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER
@@ -153,12 +297,6 @@
                zak2    = 10**(zck2)
                zakb    = EXP( zckb  )
-               zakp1   = EXP( zckp1 )
-               zakp2   = EXP( zckp2 )
-               zakp3   = EXP( zckp3 )
-               zaksi   = EXP( zcksi )
                zakw    = EXP( zckw )
                zaksp1  = 10**(zaksp0)
-               zaks    = exp( zcks )
-               zakf    = exp( zckf )
 
                ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970)
@@ -178,56 +316,29 @@
                !        CORRECTION AFTER CULBERSON AND PYTKOWICZ (1968)
                !        (CF. BROECKER ET AL., 1982)
+
+               zbuf1  = -(devk1(1)+devk2(1)*ztc+devk3(1)*ztc*ztc)
+               zbuf2  = 0.5*(devk4(1)+devk5(1)*ztc)
+               ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  =     - ( devk1(2) + devk2(2) * ztc + devk3(2) * ztc * ztc )
+               zbuf2  = 0.5 * ( devk4(2) + devk5(2) * ztc )
+               ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
                zbuf1  =     - ( devk1(3) + devk2(3) * ztc + devk3(3) * ztc * ztc )
                zbuf2  = 0.5 * ( devk4(3) + devk5(3) * ztc )
                akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
 
-               zbuf1  = -(devk1(1)+devk2(1)*ztc+devk3(1)*ztc*ztc)
-               zbuf2  = 0.5*(devk4(1)+devk5(1)*ztc)
-               ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk1(2) + devk2(2) * ztc + devk3(2) * ztc * ztc )
-               zbuf2  = 0.5 * ( devk4(2) + devk5(2) * ztc )
-               ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
                zbuf1  =     - ( devk1(4) + devk2(4) * ztc + devk3(4) * ztc * ztc )
                zbuf2  = 0.5 * ( devk4(4) + devk5(4) * ztc )
-               akp13(ji,jj,jk) = zakp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk1(5) + devk2(5) * ztc + devk3(5) * ztc * ztc )
-               zbuf2  = 0.5 * ( devk4(5) + devk5(5) * ztc )
-               akp23(ji,jj,jk) = zakp2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk1(6) + devk2(6) * ztc + devk3(6) * ztc * ztc )
-               zbuf2  = 0.5 * ( devk4(6) + devk5(6) * ztc )
-               akp33(ji,jj,jk) = zakp3 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk1(7) + devk2(7) * ztc + devk3(7) * ztc * ztc )
-               zbuf2  = 0.5 * ( devk4(7) + devk5(7) * ztc )
                akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
 
-               !  Ksi
-               !            aksi3(ji,jj,jk) = zaksi
-               !
-               !  Or using coefficient of borates (cf millero 95+ corrected version html doc co2sys)
-               !  "deltaVsi and deltaKsi have been estimated from the value of boric acid"
-               zbuf1  =     - ( devk1(3) + devk2(3) * ztc + devk3(3) * ztc * ztc )
-               zbuf2  = 0.5 * ( devk4(3) + devk5(3) * ztc )
-               aksi3(ji,jj,jk) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
 
                ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE 
                !        AS FUNCTION OF PRESSURE FOLLOWING MILLERO
                !        (P. 1285) AND BERNER (1976)
-               zbuf1  =     - ( devk1(8) + devk2(8) * ztc + devk3(8) * ztc * ztc )
-               zbuf2  = 0.5 * ( devk4(8) + devk5(8) * ztc )
+               zbuf1  =     - ( devk1(5) + devk2(5) * ztc + devk3(5) * ztc * ztc )
+               zbuf2  = 0.5 * ( devk4(5) + devk5(5) * ztc )
                aksp(ji,jj,jk) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
 
-               !  Pressure correction for sulfate and fluoride
-               zbuf1  =     - ( devk1(9) + devk2(9) * ztc + devk3(9) * ztc * ztc )
-               zbuf2  = 0.5 * ( devk4(9) + devk5(9) * ztc )
-               aks3(ji,jj,jk) = zaks   * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk1(10) + devk2(10) * ztc + devk3(10) * ztc * ztc )
-               zbuf2  = 0.5 * ( devk4(10) + devk5(10) * ztc )
-               akf3(ji,jj,jk) = zakf   * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
 
                ! TOTAL BORATE CONCENTR. [MOLES/L]
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zdiat.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zdiat.F90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zdiat.F90	(revision 858)
@@ -16,4 +16,5 @@
    USE trp_trc         ! 
    USE sms             ! 
+   USE prtctl_trc
 
    IMPLICIT NONE
@@ -40,6 +41,9 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zfact, zstep, zcompadi
+      INTEGER  ::  ji, jj, jk
+      REAL(wp) ::  zfactfe,zfactsi,zfactch, zstep, zcompadi
+      REAL(wp) ::  zrespp2, ztortp2, zmortp2
+      CHARACTER (len=25) :: charout
+ 
       !!---------------------------------------------------------------------
 
@@ -57,5 +61,4 @@
 
                zcompadi = MAX( ( trn(ji,jj,jk,jpdia) - 1e-8), 0. )
-               zfact    = 1. / ( trn(ji,jj,jk,jpdia) + rtrn )
 
 !    Aggregation term for diatoms is increased in case of nutrient
@@ -64,34 +67,52 @@
 !     ------------------------------------------------------------
 
-               respp2 (ji,jj,jk) = 1.e6 * zstep * (  wchl + wchld * ( 1.- xlimdia(ji,jj,jk) )  )    &
+               zrespp2  = 1.e6 * zstep * (  wchl + wchld * ( 1.- xlimdia(ji,jj,jk) )  )    &
 # if defined key_off_degrad
-                  &              * facvol(ji,jj,jk)       &
+                  &       * facvol(ji,jj,jk)       &
 # endif
-                  &              * zdiss(ji,jj,jk) * zcompadi * trn(ji,jj,jk,jpdia)
+                  &       * xdiss(ji,jj,jk) * zcompadi * trn(ji,jj,jk,jpdia)
                                                                                
-               respds (ji,jj,jk) = respp2(ji,jj,jk) * trn(ji,jj,jk,jpbsi) * zfact
-
-               respdf (ji,jj,jk) = respp2(ji,jj,jk) * trn(ji,jj,jk,jpdfe) * zfact
-                                                                               
-               respdch(ji,jj,jk) = respp2(ji,jj,jk) * trn(ji,jj,jk,jpdch) * zfact
 
 !     Phytoplankton mortality. 
 !     ------------------------
-               tortp2 (ji,jj,jk) = mprat2 * zstep * trn(ji,jj,jk,jpdia)     &
+               ztortp2  = mprat2 * zstep * trn(ji,jj,jk,jpdia)     &
 # if defined key_off_degrad
-                  &              * facvol(ji,jj,jk)       &
+                  &        * facvol(ji,jj,jk)       &
 # endif
-                  &              / ( xkmort + trn(ji,jj,jk,jpdia) ) * zcompadi
+                  &      / ( xkmort + trn(ji,jj,jk,jpdia) ) * zcompadi
 
-               tortds (ji,jj,jk) = tortp2(ji,jj,jk) * trn(ji,jj,jk,jpbsi) * zfact
+                zmortp2 = zrespp2 + ztortp2
 
-               tortdf (ji,jj,jk) = tortp2(ji,jj,jk) * trn(ji,jj,jk,jpdfe) * zfact
+!   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 )
+               zfactsi = trn(ji,jj,jk,jpbsi) / ( trn(ji,jj,jk,jpdia) + rtrn )
 
-               tortdch(ji,jj,jk) = tortp2(ji,jj,jk) * trn(ji,jj,jk,jpdch) * zfact
-
+               tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) - zmortp2 
+               tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) - zmortp2 * zfactch
+               tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) - zmortp2 * zfactfe
+               tra(ji,jj,jk,jpbsi) = tra(ji,jj,jk,jpbsi) - zmortp2 * zfactsi
+               tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) + zmortp2 * zfactsi
+#if defined key_kriest
+               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zmortp2  
+               tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + ztortp2 * xkr_ndiat + zrespp2 * xkr_naggr
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zmortp2 * zfactfe
+#else
+               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zrespp2 + 0.5 * ztortp2
+               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + 0.5 * ztortp2
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 0.5 * ztortp2 * zfactfe
+               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + ( zrespp2 + 0.5 * ztortp2 ) * zfactfe
+#endif
             END DO
          END DO
       END DO
       !
+        IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+         WRITE(charout, FMT="('diat')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+       ENDIF
+             
    END SUBROUTINE p4z_diat
 
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zflx.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zflx.F90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zflx.F90	(revision 858)
@@ -19,4 +19,5 @@
    USE trp_trc
    USE sms
+   USE prtctl_trc
 
    IMPLICIT NONE
@@ -24,4 +25,8 @@
 
    PUBLIC   p4z_flx    ! called in p4zprg.F90
+
+   REAL(wp) :: &  ! pre-industrial atmospheric [co2] (ppm) 	
+     atcox = 0.20946
+
 
    !!* Substitution
@@ -35,5 +40,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_flx
+   SUBROUTINE p4z_flx ( kt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_flx  ***
@@ -43,9 +48,12 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jrorr
-      REAL(wp) ::   zpdtan, zttc, zws
+      INTEGER, INTENT(in) :: kt
+      INTEGER  ::   ji, jj, jrorr, nspyr
+      REAL(wp) ::   zttc, zws
       REAL(wp) ::   zfld, zflu, zoxy16, zflu16, zfact
       REAL(wp) ::   zph, zah2, zbot, zdic, zalk, zschmitto2, zalka, zschmittco2
-      REAL(wp), DIMENSION(jpi,jpj) ::   zkgco2, zkgo2, zh2co3
+      REAL(wp), DIMENSION(jpi,jpj) ::   zkgco2, zkgo2, zh2co3, zqcumtemp 
+      CHARACTER (len=25) :: charout
+
       !!---------------------------------------------------------------------
 
@@ -55,5 +63,5 @@
       ! -----------------------------------------------------
 
-      zpdtan = raass / rdt
+      nspyr = INT( raass / rdt )
 
       ! SURFACE CHEMISTRY (PCO2 AND [H+] IN
@@ -63,5 +71,7 @@
       DO jrorr = 1, 10
 
+!CDIR NOVERRCHK
          DO jj = 1, jpj
+!CDIR NOVERRCHK
             DO ji = 1, jpi
 
@@ -94,5 +104,7 @@
       ! -------------------------------------------
 
-      DO jj = 1, jpj
+!CDIR NOVERRCHK
+      DO jj = 1, jpj
+!CDIR NOVERRCHK
          DO ji = 1, jpi
 
@@ -140,12 +152,15 @@
       END DO
 
+      zqcumtemp(:,:) = 0.
       DO jj = 1, jpj
          DO ji = 1, jpi
 
             ! Compute CO2 flux for the sea and air
-            zfld = atcco2 * tmask(ji,jj,1) * chemc(ji,jj,3) * zkgco2(ji,jj)
+            zfld = atcco2 * tmask(ji,jj,1) * chemc(ji,jj,1) * zkgco2(ji,jj)
             zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj)
             tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) / fse3t(ji,jj,1)
 
+            zqcumtemp(ji,jj) = ( zfld - zflu ) * rfact &
+               &             * e1t(ji,jj) * e2t(ji,jj) * tmask(ji,jj,1) * 1000.
             ! Compute O2 flux 
             zoxy16 = trn(ji,jj,1,jpoxy)
@@ -158,9 +173,33 @@
             trc2d(ji,jj,2) = zflu16 * 1000.
             trc2d(ji,jj,3) = zkgco2(ji,jj)
-            trc2d(ji,jj,4) = atcco2 - zh2co3(ji,jj) / ( chemc(ji,jj,3) + rtrn )
+            trc2d(ji,jj,4) = atcco2 - zh2co3(ji,jj) / ( chemc(ji,jj,1) + rtrn )
 # endif
          END DO
       END DO
       !
+
+      DO jj = 1, jpj
+        DO ji = 1, jpi
+           qcumul(1) = qcumul(1) + zqcumtemp(ji,jj) * tmask_i(ji,jj)
+        END DO
+      END DO
+
+      IF( MOD( kt, nspyr ) == 0 ) THEN
+        WRITE(numout,*) ' Atmospheric pCO2    :'
+        WRITE(numout,*) '-------------------- : ',kt,'  ',atcco2
+        WRITE(numout,*) '(ppm)'
+        WRITE(numout,*) 'Total Flux of Carbon :'
+        WRITE(numout,*) '-------------------- : ',qcumul(1)*12./1E15
+        WRITE(numout,*) '(GtC/an)'
+        qcumul(1) = 0.
+      ENDIF
+
+       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+         WRITE(charout, FMT="('flx ')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+       ENDIF
+
+
    END SUBROUTINE p4z_flx
 
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zlim.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zlim.F90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zlim.F90	(revision 858)
@@ -44,4 +44,5 @@
       REAL(wp) ::   zlim1, zlim2, zlim3, zlim4, zno3, zferlim
       REAL(wp) ::   zconctemp, zconctemp2, zconctempn, zconctempn2
+      REAL(wp) ::   ztemp, zdenom
       !!---------------------------------------------------------------------
 
@@ -70,12 +71,12 @@
             DO ji = 1, jpi
                zconctemp   = MAX( 0.e0 , trn(ji,jj,jk,jpdia)-5e-7 )
-               zconctemp2  = MIN( 5.e-7, trn(ji,jj,jk,jpdia)      )
+               zconctemp2  = trn(ji,jj,jk,jpdia) - zconctemp
                zconctempn  = MAX( 0.e0 , trn(ji,jj,jk,jpphy)-1e-6 )
-               zconctempn2 = MIN( 1.e-6, trn(ji,jj,jk,jpphy)      )
+               zconctempn2 = trn(ji,jj,jk,jpphy) - zconctempn
                concdfe(ji,jj,jk) = ( zconctemp2 * conc3 + 0.4e-9 * zconctemp)   &
-                  &              / ( zconctemp2 + zconctemp + rtrn )
+                  &              / ( trn(ji,jj,jk,jpdia) + rtrn )
                concdfe(ji,jj,jk) = MAX( conc3, concdfe(ji,jj,jk) )
                concnfe(ji,jj,jk) = ( zconctempn2 * conc2 + 0.08e-9 * zconctempn)   &
-                  &              / ( zconctempn2 + zconctempn + rtrn )
+                  &              / ( trn(ji,jj,jk,jpphy) + rtrn )
                concnfe(ji,jj,jk) = MAX( conc2, concnfe(ji,jj,jk) )
             END DO
@@ -90,11 +91,9 @@
 !      Small flagellates
 !      -----------------------------------------------
+              zdenom = 1. / &
+                  & ( conc0 * concnnh4 + concnnh4 * trn(ji,jj,jk,jpno3) + conc0 * trn(ji,jj,jk,jpnh4) )
+               xnanono3(ji,jj,jk) = trn(ji,jj,jk,jpno3) * concnnh4 * zdenom
+               xnanonh4(ji,jj,jk) = trn(ji,jj,jk,jpnh4) * conc0    * zdenom
 
-               xnanono3(ji,jj,jk) = trn(ji,jj,jk,jpno3) * concnnh4           &
-                  &               / ( conc0 * concnnh4 + concnnh4 * trn(ji,jj,jk,jpno3)   &
-                  &                                    + conc0    * trn(ji,jj,jk,jpnh4) )
-               xnanonh4(ji,jj,jk) = trn(ji,jj,jk,jpnh4) * conc0              &
-                  &               / ( conc0 * concnnh4 + concnnh4 * trn(ji,jj,jk,jpno3)   &
-                  &                                    + conc0    * trn(ji,jj,jk,jpnh4) )
                zlim1 = xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk)
                zlim2 = trn(ji,jj,jk,jppo4) / ( trn(ji,jj,jk,jppo4) + concnnh4          ) 
@@ -116,11 +115,9 @@
 !   Michaelis-Menten Limitation term for nutrients Diatoms
 !   ----------------------------------------------
+              zdenom = 1. / &
+                  & ( conc1 * concdnh4 + concdnh4 * trn(ji,jj,jk,jpno3) + conc1 * trn(ji,jj,jk,jpnh4) )
 
-               xdiatno3(ji,jj,jk) = trn(ji,jj,jk,jpno3) * concdnh4                         &
-                  &               / ( conc1  * concdnh4 + concdnh4 * trn(ji,jj,jk,jpno3)   &
-                  &                                     + conc1    * trn(ji,jj,jk,jpnh4) )
-               xdiatnh4(ji,jj,jk) = trn(ji,jj,jk,jpnh4) * conc1                            &
-                  &               / ( conc1  * concdnh4 + concdnh4 * trn(ji,jj,jk,jpno3)   &
-                  &                                     + conc1   * trn(ji,jj,jk,jpnh4) )
+               xdiatno3(ji,jj,jk) = trn(ji,jj,jk,jpno3) * concdnh4 * zdenom
+               xdiatnh4(ji,jj,jk) = trn(ji,jj,jk,jpnh4) * conc1    * zdenom 
 
                zlim1 = xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk)
@@ -130,4 +127,21 @@
                xlimdia(ji,jj,jk) = MIN( zlim1, zlim2, zlim3, zlim4 )
 
+            END DO
+         END DO
+      END DO
+
+
+      ! Compute the fraction of nanophytoplankton that is made of calcifiers
+      ! --------------------------------------------------------------------
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ztemp = MAX( 0., tn(ji,jj,jk) )
+               xfracal(ji,jj,jk) = caco3r * xlimphy(ji,jj,jk)   &
+                  &                       * MAX( 0.0001, ztemp / ( 2.+ ztemp ) )   &
+                  &                       * MAX( 1., trn(ji,jj,jk,jpphy) * 1.e6 / 2. )
+               xfracal(ji,jj,jk) = MIN( 0.8 , xfracal(ji,jj,jk) )
+               xfracal(ji,jj,jk) = MAX( 0.01, xfracal(ji,jj,jk) )
             END DO
          END DO
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zlys.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zlys.F90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zlys.F90	(revision 858)
@@ -19,4 +19,5 @@
    USE trp_trc
    USE sms
+   USE prtctl_trc
 
    IMPLICIT NONE
@@ -24,4 +25,9 @@
 
    PUBLIC   p4z_lys    ! called in p4zprg.F90
+
+   !! * Module variables
+
+   REAL(wp) :: &
+      calcon = 1.03E-2        ! mean calcite concentration [Ca2+] in sea water [mole/kg solution]
 
    !!----------------------------------------------------------------------
@@ -48,4 +54,5 @@
       REAL(wp) ::   zomegaca, zexcess, zexcess0
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zco3
+      CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
 
@@ -57,6 +64,9 @@
       DO jn = 1, 5                               !  BEGIN OF ITERATION
          !
+!CDIR NOVERRCHK
          DO jk = 1, jpkm1
+!CDIR NOVERRCHK
             DO jj = 1, jpj
+!CDIR NOVERRCHK
                DO ji = 1, jpi
 
@@ -131,9 +141,15 @@
 
 # if defined key_trc_dia3d
-         trc3d(:,:,:,1) = rhop(:,:,:)
+         trc3d(:,:,:,1) = hi(:,:,:)
          trc3d(:,:,:,2) = zco3(:,:,:)
          trc3d(:,:,:,3) = aksp(:,:,:) / calcon
 # endif
       !
+       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+         WRITE(charout, FMT="('lys ')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+       ENDIF
+
    END SUBROUTINE p4z_lys
 
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zmeso.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zmeso.F90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zmeso.F90	(revision 858)
@@ -16,4 +16,5 @@
    USE trp_trc         ! 
    USE sms             ! 
+   USE prtctl_trc
 
    IMPLICIT NONE
@@ -40,7 +41,15 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zcompadi, zcompaph, zcompapoc, zcompaz
-      REAL(wp) ::   zfact, zstep, zcompam, zdenom, zgraze2
+      INTEGER  :: ji, jj, jk
+      REAL(wp) :: zcompadi, zcompaph, zcompapoc, zcompaz
+      REAL(wp) :: zfact, zstep, zcompam, zdenom, zgraze2
+      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
+      CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
 
@@ -60,6 +69,6 @@
 !     Respiration rates of both zooplankton
 !     -------------------------------------
-               respz2(ji,jj,jk) = resrat2 * zfact * ( 1. + 3. * nitrfac(ji,jj,jk) )        &
-                  &             * trn(ji,jj,jk,jpmes) / ( xkmort + trn(ji,jj,jk,jpmes) )
+               zrespz2(ji,jj,jk)  = 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
@@ -67,5 +76,5 @@
 !     mimic predation.
 !     ---------------------------------------------------------------
-               tortz2(ji,jj,jk) = mzrat2 * 1.e6 * zfact * trn(ji,jj,jk,jpmes)
+               ztortz2(ji,jj,jk) = mzrat2 * 1.e6 * zfact * trn(ji,jj,jk,jpmes)
                !
             END DO
@@ -73,8 +82,8 @@
       END DO
 
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-
+
+      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 )
@@ -95,32 +104,28 @@
                   &     * trn(ji,jj,jk,jpmes)
 
-               grazd  (ji,jj,jk) = zgraze2 * xprefc   * zcompadi
-               grazz  (ji,jj,jk) = zgraze2 * xprefz   * zcompaz
-               grazn  (ji,jj,jk) = zgraze2 * xprefp   * zcompaph
-               grazpoc(ji,jj,jk) = zgraze2 * xprefpoc * zcompapoc
-
-               graznf (ji,jj,jk) = grazn  (ji,jj,jk) * trn(ji,jj,jk,jpnfe) / (trn(ji,jj,jk,jpphy) + rtrn)
-
-               graznch(ji,jj,jk) = grazn  (ji,jj,jk) * trn(ji,jj,jk,jpnch) / (trn(ji,jj,jk,jpphy) + rtrn)
-
-               grazs  (ji,jj,jk) = grazd  (ji,jj,jk) * trn(ji,jj,jk,jpbsi) / (trn(ji,jj,jk,jpdia) + rtrn)
-
-               grazf  (ji,jj,jk) = grazd  (ji,jj,jk) * trn(ji,jj,jk,jpdfe) / (trn(ji,jj,jk,jpdia) + rtrn)
-
-               grazdch(ji,jj,jk) = grazd  (ji,jj,jk) * trn(ji,jj,jk,jpdch) / (trn(ji,jj,jk,jpdia) + rtrn)
-
-               grazpof(ji,jj,jk) = grazpoc(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
-
+               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
-               grazffe(ji,jj,jk) = 5.e3 * zstep * wsbio4(ji,jj,jk)          &
+               zgrazffe(ji,jj,jk) = 5.e3 * zstep * wsbio4(ji,jj,jk)          &
 #  if defined key_off_degrad
                   &     * facvol(ji,jj,jk)          &
@@ -128,9 +133,9 @@
                   &     * tgfunc2(ji,jj,jk) * trn(ji,jj,jk,jpgoc) * trn(ji,jj,jk,jpmes)
 
-               grazfff(ji,jj,jk) = grazffe(ji,jj,jk)         &
+               zgrazfff(ji,jj,jk) = zgrazffe(ji,jj,jk)       &
                   &     * trn(ji,jj,jk,jpbfe) / (trn(ji,jj,jk,jpgoc) + rtrn)
 # else
 ! KRIEST3
-               grazffe(ji,jj,jk) = 0.5 * 1.3e-2 / 5.5e-7 * 0.3 * zstep * wsbio3(ji,jj,jk)     &
+               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
@@ -139,5 +144,5 @@
                   &     /  (trn(ji,jj,jk,jppoc) * 1.e7 + 0.1)
 
-!!C        grazffe(ji,jj,jk) = 5.e3 * zstep * wsbio3(ji,jj,jk)
+!!C        zgrazffe(ji,jj,jk) = 5.e3 * zstep * wsbio3(ji,jj,jk)
 !!C     &     * tgfunc2(ji,jj,jk) * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpmes)
 !!C#    if defined key_off_degrad
@@ -145,7 +150,55 @@
 !!C#    endif
 
-               grazfff(ji,jj,jk) = grazffe(ji,jj,jk)       &
+               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
+      
+
+      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_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.)  )
+#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
+
+               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,jpoxy) = tra(ji,jj,jk,jpoxy) - o2ut * zgrarem2 * sigma2
+               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zgrafer2
+               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zgrarem2 * sigma2
+               
+#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_nmeso
+#else
+               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zgrapoc2
+#endif
             END DO
          END DO
@@ -155,44 +208,62 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-
-!    Mesozooplankton efficiency
-!    --------------------------
-               grarem2(ji,jj,jk) = (grazd(ji,jj,jk) + grazz(ji,jj,jk)              &
-                  &     + grazn(ji,jj,jk) + grazpoc(ji,jj,jk) + grazffe(ji,jj,jk))   &
-                  &     * (1.-epsher2-unass2)
-#if ! defined key_kriest
-               grafer2(ji,jj,jk) = (grazf(ji,jj,jk) + graznf(ji,jj,jk)                       &
-                  &     + grazz(ji,jj,jk) * ferat3 + grazpof(ji,jj,jk)                       &
-                  &     + grazfff(ji,jj,jk)) * (1.-epsher2-unass2)                       &
-                  &     + epsher2 * (grazd(ji,jj,jk) * MAX(                       &
-                  &    (trn(ji,jj,jk,jpdfe) / (trn(ji,jj,jk,jpdia) + rtrn)                       &
-                  &    -ferat3),0.) + grazn(ji,jj,jk) * MAX(                       &
-                  &    (trn(ji,jj,jk,jpnfe) / (trn(ji,jj,jk,jpphy) + rtrn)                       &
-                  &    -ferat3),0.) + grazpoc(ji,jj,jk) * MAX(                       &
-                  &    (trn(ji,jj,jk,jpsfe) / (trn(ji,jj,jk,jppoc) + rtrn)                       &
-                  &    -ferat3),0.) + grazffe(ji,jj,jk) * MAX(                       &
-                  &    (trn(ji,jj,jk,jpbfe) / (trn(ji,jj,jk,jpgoc) + rtrn)                       &
-                  &    -ferat3),0.) )
-#else
-               grafer2(ji,jj,jk) = (grazf(ji,jj,jk) + graznf(ji,jj,jk)                       &
-                  &     + grazz(ji,jj,jk) * ferat3 + grazpof(ji,jj,jk)                       &
-                  &     + grazfff(ji,jj,jk)) * (1.-epsher2-unass2)                       &
-                  &     + epsher2 * (grazd(ji,jj,jk) * MAX(                       &
-                  &    (trn(ji,jj,jk,jpdfe) / (trn(ji,jj,jk,jpdia) + rtrn)                       &
-                  &    -ferat3),0.) + grazn(ji,jj,jk) * MAX(                       &
-                  &    (trn(ji,jj,jk,jpnfe) / (trn(ji,jj,jk,jpphy) + rtrn)                       &
-                  &    -ferat3),0.) + grazpoc(ji,jj,jk) * MAX(                       &
-                  &    (trn(ji,jj,jk,jpsfe) / (trn(ji,jj,jk,jppoc) + rtrn)                       &
-                  &    -ferat3),0.) + grazffe(ji,jj,jk) * MAX(                       &
-                  &    (trn(ji,jj,jk,jpsfe) / (trn(ji,jj,jk,jppoc) + rtrn)                       &
-                  &    -ferat3),0.) )
-#endif
-               grapoc2(ji,jj,jk) = (grazd(ji,jj,jk) + grazz(ji,jj,jk)                     &
-                  &     + grazn(ji,jj,jk) + grazpoc(ji,jj,jk) + grazffe(ji,jj,jk)) * unass2
                !
+               !   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) * part * unass2 * zgrazn(ji,jj,jk)
+
+               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprcaca
+               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2. * zprcaca
+               tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) + zprcaca
+#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_nmeso &
+                  &                  - zgrazffe(ji,jj,jk)   * 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)
+#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,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)
+#endif
+
             END DO
          END DO
       END DO
       !
+       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+         WRITE(charout, FMT="('meso')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+       ENDIF
+
    END SUBROUTINE p4z_meso
 
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zmicro.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zmicro.F90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zmicro.F90	(revision 858)
@@ -16,4 +16,5 @@
    USE trp_trc         ! 
    USE sms             ! 
+   USE prtctl_trc
 
    IMPLICIT NONE
@@ -40,8 +41,14 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zcompadi, zcompadi2, zcompaz , zcompaph, zcompapoc
-      REAL(wp) ::   zgraze  , zdenom  , zdenom2
-      REAL(wp) ::   zfact   , zstep   , zinano , zidiat, zipoc
+      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) :: 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
+      CHARACTER (len=25) :: charout
+
       !!---------------------------------------------------------------------
 
@@ -62,5 +69,5 @@
 !     -------------------------------------
 
-               respz(ji,jj,jk) = resrat * zfact  * ( 1.+ 3.* nitrfac(ji,jj,jk) )     &
+               zrespz(ji,jj,jk) = resrat * zfact  * ( 1.+ 3.* nitrfac(ji,jj,jk) )     &
                   &            * trn(ji,jj,jk,jpzoo) / ( xkmort + trn(ji,jj,jk,jpzoo) )
 
@@ -69,22 +76,22 @@
 !     mimic predation.
 !     ---------------------------------------------------------------
-
-               tortz(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
-
+               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
                zcompadi  = MAX( ( trn(ji,jj,jk,jpdia) - 1.e-8 ), 0.e0 )
                zcompadi2 = MIN( zcompadi, 5.e-7 )
                zcompaph  = MAX( ( trn(ji,jj,jk,jpphy) - 2.e-7 ), 0.e0 )
                zcompapoc = MAX( ( trn(ji,jj,jk,jppoc) - 1.e-8 ), 0.e0 )
-
-!     Microzooplankton grazing
-!     ------------------------
+               
+               !     Microzooplankton grazing
+               !     ------------------------
                zdenom2 = 1./ ( zprefp * zcompaph + zprefc * zcompapoc + zprefd * zcompadi2 + rtrn )
 
@@ -101,23 +108,52 @@
                zdenom = 1./ ( xkgraz + zinano * zcompaph + zipoc * zcompapoc + zidiat * zcompadi2 )
 
-               grazp (ji,jj,jk) = zgraze * zinano * zcompaph * zdenom
-               grazm (ji,jj,jk) = zgraze * zipoc  * zcompapoc * zdenom
-               grazsd(ji,jj,jk) = zgraze * zidiat * zcompadi2 * zdenom
-
-               grazpf (ji,jj,jk) = grazp (ji,jj,jk) * trn(ji,jj,jk,jpnfe) / (trn(ji,jj,jk,jpphy) + rtrn)
-
-               grazpch(ji,jj,jk) = grazp (ji,jj,jk) * trn(ji,jj,jk,jpnch) / (trn(ji,jj,jk,jpphy) + rtrn)
-
-               grazmf (ji,jj,jk) = grazm (ji,jj,jk) * trn(ji,jj,jk,jpsfe) / (trn(ji,jj,jk,jppoc) + rtrn)
-
-               grazsf (ji,jj,jk) = grazsd(ji,jj,jk) * trn(ji,jj,jk,jpdfe) / (trn(ji,jj,jk,jpdia) + rtrn)
-
-               grazss (ji,jj,jk) = grazsd(ji,jj,jk) * trn(ji,jj,jk,jpbsi) / (trn(ji,jj,jk,jpdia) + rtrn)
-
-               grazsch(ji,jj,jk) = grazsd(ji,jj,jk) * trn(ji,jj,jk,jpdch) / (trn(ji,jj,jk,jpdia) + rtrn)
-
-            END DO
-         END DO
-      END DO
+               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
+      
+
+      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
+
+               !  Update of the TRA arrays
+               !  ------------------------
+
+               tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zgrarem * sigma1
+               tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zgrarem * sigma1
+               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zgrarem * (1.-sigma1)
+               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,jpdic) = tra(ji,jj,jk,jpdic) + zgrarem * sigma1
+#if defined key_kriest
+               tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zgrapoc * xkr_ndiat
+#endif
+            END DO
+         END DO
+      END DO
+
+!
+!   Update the arrays TRA which contain the biological sources and sinks
+!   --------------------------------------------------------------------
 
       DO jk = 1, jpkm1
@@ -125,25 +161,40 @@
             DO ji = 1, jpi
 
-!    Various remineralization and excretion terms
-!    --------------------------------------------
-
-               grarem(ji,jj,jk) = (  grazp(ji,jj,jk) + grazm (ji,jj,jk)        &
-                  &                                  + grazsd(ji,jj,jk)  ) * ( 1.- epsher - unass )
-
-               grafer(ji,jj,jk) = (  grazpf(ji,jj,jk) + grazsf(ji,jj,jk)      &
-                  &                                   + grazmf(ji,jj,jk)  ) * ( 1.- epsher - unass )   &
-                  &             + (  grazm (ji,jj,jk) * MAX( (trn(ji,jj,jk,jpsfe) /     &
-                  &                                          (trn(ji,jj,jk,jppoc) + rtrn) - ferat3), 0.e0 )   &
-                  &                + grazp (ji,jj,jk) * MAX( (trn(ji,jj,jk,jpnfe) /   &
-                  &                                          (trn(ji,jj,jk,jpphy) + rtrn) - ferat3), 0.e0 )   &
-                  &                + grazsd(ji,jj,jk) * MAX( (trn(ji,jj,jk,jpdfe) /   &
-                  &                                         (trn(ji,jj,jk,jpdia) + rtrn) - ferat3), 0.e0 )  ) * epsher 
-
-               grapoc(ji,jj,jk) = (  grazp(ji,jj,jk) + grazm(ji,jj,jk) + grazsd(ji,jj,jk)  ) * unass
-
+               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) * part * unass * zgrazp(ji,jj,jk)
+               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprcaca
+               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal)- 2. * zprcaca
+               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_ndiat
+#endif
             END DO
          END DO
       END DO
       !
+       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+         WRITE(charout, FMT="('micro')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+       ENDIF
+
    END SUBROUTINE p4z_micro
 
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4znano.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4znano.F90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4znano.F90	(revision 858)
@@ -16,4 +16,5 @@
    USE trp_trc         ! 
    USE sms             ! 
+   USE prtctl_trc
 
    IMPLICIT NONE
@@ -40,6 +41,9 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zfact, zstep, zcompaph
+      INTEGER  :: ji, jj, jk
+      REAL(wp) :: zstep, zcompaph
+      REAL(wp) :: zfactfe,zfactch,zprcaca,zfracal
+      REAL(wp) :: ztortp,zrespp,zmortp
+      CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
 
@@ -53,18 +57,13 @@
 
                zcompaph = MAX( ( trn(ji,jj,jk,jpphy) - 1e-8 ), 0.e0 )
-               zfact    = 1./ ( trn(ji,jj,jk,jpphy) + rtrn )
 
 !     Squared mortality of Phyto similar to a sedimentation term during
 !     blooms (Doney et al. 1996)
 !     -----------------------------------------------------------------
-               respp(ji,jj,jk) = wchl * 1.e6 * zstep * zdiss(ji,jj,jk)   &
+               zrespp = wchl * 1.e6 * zstep * xdiss(ji,jj,jk)   &
 # if defined key_off_degrad
-                  &            * facvol(ji,jj,jk)     &
+                  &        * facvol(ji,jj,jk)     &
 # endif
-                  &            * zcompaph * trn(ji,jj,jk,jpphy)
-                                                                               
-               respnf (ji,jj,jk) = respp(ji,jj,jk) * trn(ji,jj,jk,jpnfe) * zfact
-                                                                               
-               respnch(ji,jj,jk) = respp(ji,jj,jk) * trn(ji,jj,jk,jpnch) * zfact
+                  &        * zcompaph * trn(ji,jj,jk,jpphy)
 
 !     Phytoplankton mortality. This mortality loss is slightly
@@ -72,18 +71,46 @@
 !     as observed for instance in case of iron limitation.
 !     ----------------------------------------------------------
-               tortp  (ji,jj,jk) = mprat * zstep * trn(ji,jj,jk,jpphy)          &
+               ztortp = mprat * zstep * trn(ji,jj,jk,jpphy)          &
 # if defined key_off_degrad
-                  &              * facvol(ji,jj,jk)     &
+                  &          * facvol(ji,jj,jk)     &
 # endif
-                  &              / ( xkmort + trn(ji,jj,jk,jpphy) ) * zcompaph
+                  &   / ( xkmort + trn(ji,jj,jk,jpphy) ) * zcompaph
+
                                                                                
-               tortnf (ji,jj,jk) = tortp(ji,jj,jk) * trn(ji,jj,jk,jpnfe) * zfact
-                                                                               
-               tortnch(ji,jj,jk) = tortp(ji,jj,jk) * trn(ji,jj,jk,jpnch) * zfact
+               zmortp = zrespp + ztortp
 
+               !   Update the arrays TRA which contains the biological sources and sinks
+               
+               zfactfe = trn(ji,jj,jk,jpnfe)/(trn(ji,jj,jk,jpphy)+rtrn)
+               zfactch = trn(ji,jj,jk,jpnch)/(trn(ji,jj,jk,jpphy)+rtrn)
+               
+               tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) - zmortp
+               tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) - zmortp * zfactch
+               tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) - zmortp * zfactfe
+               zprcaca = xfracal(ji,jj,jk) * zmortp
+               zfracal = 0.5 * xfracal(ji,jj,jk)
+               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprcaca
+               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2. * zprcaca
+               tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) + zprcaca
+#if defined key_kriest
+               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zmortp
+               tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + ztortp * xkr_nnano + zrespp * xkr_ndiat
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zmortp * zfactfe
+#else
+               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zfracal * zmortp
+               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + ( 1. - zfracal ) * zmortp
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ( 1. - zfracal ) * zmortp * zfactfe
+               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zfracal * zmortp * zfactfe
+#endif
             END DO
          END DO
       END DO
       !
+       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+         WRITE(charout, FMT="('nano')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+       ENDIF
+
    END SUBROUTINE p4z_nano
 
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zopt.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zopt.F90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zopt.F90	(revision 858)
@@ -91,5 +91,7 @@
       END DO
 
+!CDIR NOVERRCHK
       DO jj = 1,jpj
+!CDIR NOVERRCHK
          DO ji = 1,jpi
 
@@ -108,6 +110,9 @@
       END DO
 
+!CDIR NOVERRCHK
       DO jk = 2, jpkm1
+!CDIR NOVERRCHK
           DO jj = 1, jpj
+!CDIR NOVERRCHK
             DO ji = 1, jpi
 
@@ -131,4 +136,7 @@
 
       etot(:,:,:) = ze1(:,:,:) + ze2(:,:,:) + ze3(:,:,:)
+      enano(:,:,:) = 2.1 * ze1(:,:,:) + 0.42 * ze2(:,:,:) + 0.4 * ze3(:,:,:)
+      ediat(:,:,:) = 1.6 * ze1(:,:,:) + 0.69 * ze2(:,:,:) + 0.7 * ze3(:,:,:)
+
 
       IF( ln_qsr_sms ) THEN
@@ -137,5 +145,7 @@
 !   ------------------------------------------------------------------------------
 
-         DO jj = 1, jpj
+!CDIR NOVERRCHK
+         DO jj = 1, jpj
+!CDIR NOVERRCHK
             DO ji = 1, jpi
 
@@ -155,6 +165,9 @@
          END DO
 
+!CDIR NOVERRCHK
          DO jk = 2, jpkm1
+!CDIR NOVERRCHK
             DO jj = 1, jpj
+!CDIR NOVERRCHK
                DO ji = 1, jpi
 
@@ -162,7 +175,7 @@
 !     -------------------------------------------------
 
-                  zblight1 = zekb(ji,jj,jk-1) * fse3t(ji,jj,jk-1)
-                  zglight1 = zekg(ji,jj,jk-1) * fse3t(ji,jj,jk-1)
-                  zrlight1 = zekr(ji,jj,jk-1) * fse3t(ji,jj,jk-1)
+                  zblight = zekb(ji,jj,jk-1) * fse3t(ji,jj,jk-1)
+                  zglight = zekg(ji,jj,jk-1) * fse3t(ji,jj,jk-1)
+                  zrlight = zekr(ji,jj,jk-1) * fse3t(ji,jj,jk-1)
 
                   ze3lum(ji,jj,jk) = ze3lum(ji,jj,jk-1) * EXP( -zblight )
@@ -182,15 +195,15 @@
 !     ---------------------------------
     
-      zmeu(:,:) = 300.e0
+      heup(:,:) = 300.e0
 
       DO jk = 2, jpkm1
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( etot(ji,jj,jk) >= 0.0043 * qsr(ji,jj) )   zmeu(ji,jj) = fsdepw(ji,jj,jk+1)
-            END DO
-         END DO
-      END DO
-
-      zmeu(:,:) = MIN( 300., zmeu(:,:) )
+               IF( etot(ji,jj,jk) >= 0.0043 * qsr(ji,jj) )   heup(ji,jj) = fsdepw(ji,jj,jk+1)
+            END DO
+         END DO
+      END DO
+
+      heup(:,:) = MIN( 300., heup(:,:) )
 
 !    Computation of the mean light over the mixed layer depth
@@ -224,6 +237,7 @@
       END DO
 
+
 # if defined key_trc_diaadd
-      trc2d(:,:,11) = zmeu(:,:)
+      trc2d(:,:,11) = heup(:,:)
 # endif
       !
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zprg.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zprg.F90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zprg.F90	(revision 858)
@@ -25,5 +25,5 @@
    USE p4zlys          ! 
    USE p4zflx          ! 
-   
+
    IMPLICIT NONE
    PRIVATE
@@ -54,4 +54,6 @@
       INTEGER ::   jnt, jn
       INTEGER ::   iyy, imm, idd
+
+
       !!---------------------------------------------------------------------
 
@@ -74,7 +76,9 @@
       DO jnt = 1, nrdttrc             ! ???
          !
-         CALL p4z_bio          ! Compute soft tissue production (POC)
+         CALL p4z_bio ( kt )         ! Compute soft tissue production (POC)
+
 
          CALL p4z_sed          ! compute soft tissue remineralisation
+
          !
          trb(:,:,:,:) = trn(:,:,:,:)
@@ -84,7 +88,8 @@
       CALL p4z_lys             ! Compute CaCO3 saturation
 
-      CALL p4z_flx             ! Compute surface fluxes
+      CALL p4z_flx( kt )             ! Compute surface fluxes
 
-      DO jn = 1, jptra
+
+      DO jn = jp_pcs0, jp_pcs1
         CALL lbc_lnk( trn(:,:,:,jn), 'T', 1. )
         CALL lbc_lnk( trb(:,:,:,jn), 'T', 1. )
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zprod.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zprod.F90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zprod.F90	(revision 858)
@@ -17,5 +17,6 @@
    USE sms             ! 
    USE p4zday          !
-   
+   USE prtctl_trc
+
    IMPLICIT NONE
    PRIVATE
@@ -23,4 +24,5 @@
    PUBLIC   p4z_prod    ! called in p4zbio.F90
 
+   
    !!* Substitution
 #  include "domzgr_substitute.h90"
@@ -33,5 +35,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_prod
+   SUBROUTINE p4z_prod( kt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_prod  ***
@@ -42,15 +44,36 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zsilfac, zfact
+      INTEGER, INTENT(in) :: kt
+      INTEGER  ::   ji, jj, jk, nspyr
+      REAL(wp) ::   zsilfac, zfact, zrfact2
       REAL(wp) ::   zprdiachl, zprbiochl, zsilim, ztn, zadap, zadap2
-      REAL(wp) ::   zlim, zsilfac2, zsiborn, zprod
-      REAL(wp) ::   zmxltst, zmxlday, zlim1
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zmixnano , zmixdiat
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpislopen, zpislope2n
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsopt    , zpislopead
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zprdia   , zprbio
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zetot2   , zpislopead2
+      REAL(wp) ::   zlim, zsilfac2, zsiborn, zprod, zetot2, zmax, zproreg, zproreg2
+      REAL(wp) ::   zmxltst, zmxlday, zlim1, zexcret, zexcret2
+      REAL(wp) ::   zpislopen  , zpislope2n
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zmixnano   , zmixdiat
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpislopead , zpislopead2
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zprdia     , zprbio, zysopt
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zprorca    , zprorca2, zprorca4
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zprorca5   , zprorca6, zprorca7
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpronew    , zpronew2
+      CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
+
+      zprorca (:,:,:) = 0.0
+      zprorca2(:,:,:) = 0.0
+      zprorca4(:,:,:) = 0.0
+      zprorca5(:,:,:) = 0.0
+      zprorca6(:,:,:) = 0.0
+      zprorca7(:,:,:) = 0.0
+      zpronew (:,:,:) = 0.0
+      zpronew2(:,:,:) = 0.0
+      zprdia  (:,:,:) = 0.0
+      zprbio  (:,:,:) = 0.0
+      zysopt  (:,:,:) = 0.0
+
+      nspyr  = INT( raass / rdt )
+
+      zexcret  = 1. - excret
+      zexcret2 = 1. - excret2
 
 !     Computation of the optimal production
@@ -65,51 +88,49 @@
       CALL p4z_day       ! Computation of the day length
 
-
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
+!CDIR NOVERRCHK
+      DO jk = 1, jpkm1
+!CDIR NOVERRCHK
+         DO jj = 1, jpj
+!CDIR NOVERRCHK
             DO ji = 1, jpi
 
 !      Computation of the P-I slope for nanos and diatoms
 !      --------------------------------------------------
-
-               ztn    = MAX( 0., tn(ji,jj,jk) - 15. )
-               zadap  = 1.+ 2.* ztn / ( 2.+ ztn )
-               zadap2 = 1.e0
-
-               zfact  = EXP( -0.21 * emoy(ji,jj,jk) )
-
-               zpislopead (ji,jj,jk) = pislope  * ( 1.+ zadap  * zfact )
-               zpislopead2(ji,jj,jk) = pislope2 * ( 1.+ zadap2 * zfact )
-
-               zpislopen(ji,jj,jk) = zpislopead(ji,jj,jk) * trn(ji,jj,jk,jpnch)               &
-                  &                / ( trn(ji,jj,jk,jpphy) * 12.                   + rtrn )   &
-                  &                / ( prmax(ji,jj,jk) * rjjss * xlimphy(ji,jj,jk) + rtrn )
-
-               zpislope2n(ji,jj,jk) = zpislopead2(ji,jj,jk) * trn(ji,jj,jk,jpdch)              &
-                  &                 / ( trn(ji,jj,jk,jpdia) * 12.                   + rtrn )   &
-                  &                 / ( prmax(ji,jj,jk) * rjjss * xlimdia(ji,jj,jk) + rtrn )
-
-            END DO
-         END DO
-      END DO
-
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
+                IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+                   ztn    = MAX( 0., tn(ji,jj,jk) - 15. )
+                   zadap  = 0.+ 1.* ztn / ( 2.+ ztn )
+                   zadap2 = 0.e0
+
+                   zfact  = EXP( -0.21 * emoy(ji,jj,jk) )
+
+                   zpislopead (ji,jj,jk) = pislope  * ( 1.+ zadap  * zfact )
+                   zpislopead2(ji,jj,jk) = pislope2 * ( 1.+ zadap2 * zfact )
+
+                   zpislopen = zpislopead(ji,jj,jk) * trn(ji,jj,jk,jpnch)               &
+                     &         / ( trn(ji,jj,jk,jpphy) * 12.                   + rtrn )   &
+                     &         / ( prmax(ji,jj,jk) * rjjss * xlimphy(ji,jj,jk) + rtrn )
+
+                   zpislope2n = zpislopead2(ji,jj,jk) * trn(ji,jj,jk,jpdch)              &
+                     &          / ( trn(ji,jj,jk,jpdia) * 12.                   + rtrn )   &
+                     &          / ( prmax(ji,jj,jk) * rjjss * xlimdia(ji,jj,jk) + rtrn )
 
 !     Computation of production function
 !     ----------------------------------
 
-               zprbio(ji,jj,jk) = prmax(ji,jj,jk) * (  1.- EXP( -zpislopen (ji,jj,jk) * etot(ji,jj,jk) )  )
-               zprdia(ji,jj,jk) = prmax(ji,jj,jk) * (  1.- EXP( -zpislope2n(ji,jj,jk) * etot(ji,jj,jk) )  )
-
-            END DO
-         END DO
-      END DO
-
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-
+                   zprbio(ji,jj,jk) = prmax(ji,jj,jk) * &
+                     &                (  1.- EXP( -zpislopen * enano(ji,jj,jk) )  )
+                   zprdia(ji,jj,jk) = prmax(ji,jj,jk) * &
+                     &                (  1.- EXP( -zpislope2n * ediat(ji,jj,jk) )  )
+               ENDIF
+            END DO
+         END DO
+      END DO
+
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+                IF( etot(ji,jj,jk) > 1.E-3 ) THEN
 !    Si/C of diatoms
 !    ------------------------
@@ -118,18 +139,18 @@
 !    to mimic the very high ratios observed in the Southern Ocean (silpot2)
 
-               zlim1  = trn(ji,jj,jk,jpsil) / ( trn(ji,jj,jk,jpsil) + xksi1 )
-               zlim   = xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk)
-
-               zsilim = MIN( zprdia(ji,jj,jk)    / ( rtrn + prmax(ji,jj,jk) ),                    &
+                  zlim1  = trn(ji,jj,jk,jpsil) / ( trn(ji,jj,jk,jpsil) + xksi1 )
+                  zlim   = xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk)
+
+                  zsilim = MIN( zprdia(ji,jj,jk)    / ( rtrn + prmax(ji,jj,jk) ),                    &
                   &          trn(ji,jj,jk,jpfer) / ( concdfe(ji,jj,jk) + trn(ji,jj,jk,jpfer) ),   &
                   &          trn(ji,jj,jk,jppo4) / ( concdnh4 + trn(ji,jj,jk,jppo4) ),            &
                   &          zlim )
-               zsilfac = 5.4 * EXP( -4.23 * zsilim ) * MAX( 0.e0, MIN( 1., 2.2 * ( zlim1 - 0.5 ) )  ) + 1.e0
-               zsiborn = MAX( 0.e0, ( trn(ji,jj,jk,jpsil) - 15.e-6 ) )
-               zsilfac2 = 1.+ 3.* zsiborn / ( zsiborn + xksi2 )
-               zsilfac = MIN( 6.4,zsilfac * zsilfac2)
-
-               zsopt(ji,jj,jk) = grosip * trn(ji,jj,jk,jpsil) / ( trn(ji,jj,jk,jpsil) + xksi1 ) * zsilfac
-
+                  zsilfac = 5.4 * EXP( -4.23 * zsilim ) * MAX( 0.e0, MIN( 1., 2.2 * ( zlim1 - 0.5 ) )  ) + 1.e0
+                  zsiborn = MAX( 0.e0, ( trn(ji,jj,jk,jpsil) - 15.e-6 ) )
+                  zsilfac2 = 1.+ 3.* zsiborn / ( zsiborn + xksi2 )
+                  zsilfac = MIN( 6.4,zsilfac * zsilfac2)
+                  zysopt(ji,jj,jk) = grosip * zlim1 * zsilfac
+
+              ENDIF
             END DO
          END DO
@@ -142,5 +163,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            zmxltst = MAX( 0.e0, hmld(ji,jj) - zmeu(ji,jj) )
+            zmxltst = MAX( 0.e0, hmld(ji,jj) - heup(ji,jj) )
             zmxlday = zmxltst**2 / rjjss
             zmixnano(ji,jj) = 1.- zmxlday / ( 1.+ zmxlday )
@@ -163,79 +184,157 @@
       END DO
 
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
 
 !      Computation of the maximum light intensity
 !      ------------------------------------------
-               zetot2(ji,jj,jk) = etot(ji,jj,jk) * 24. / ( strn(ji,jj) + rtrn )
-               IF( strn(ji,jj) < 1.e0 )   zetot2(ji,jj,jk) = etot(ji,jj,jk)
-
-            END DO
-         END DO
-      END DO
-
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-
+            IF( strn(ji,jj) < 1.e0 )   strn(ji,jj) = 24.
+         END DO
+      END DO
+
+!CDIR NOVERRCHK
+      DO jk = 1, jpkm1
+!CDIR NOVERRCHK
+         DO jj = 1, jpj
+!CDIR NOVERRCHK
+            DO ji = 1, jpi
+
+               IF( etot(ji,jj,jk) > 1.E-3 ) THEN
 !     Computation of the various production terms for nanophyto.
 !     ----------------------------------------------------------
-               zpislopen(ji,jj,jk) = zpislopead(ji,jj,jk)          &
-                  &                * trn(ji,jj,jk,jpnch) / ( rtrn + trn(ji,jj,jk,jpphy) * 12.)         &
-                  &                / ( prmax(ji,jj,jk) * rjjss * MAX( 0.1, xlimphy(ji,jj,jk) ) + rtrn )
-
-               zprbiochl = prmax(ji,jj,jk) * (  1.- EXP( -zpislopen(ji,jj,jk) * zetot2(ji,jj,jk) )  )
-
-               prorca(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * trn(ji,jj,jk,jpphy) * rfact2
-
-               pronew(ji,jj,jk) = prorca(ji,jj,jk) * xnanono3(ji,jj,jk)    &
+                  zetot2 = enano(ji,jj,jk) * 24. / ( strn(ji,jj) + rtrn )
+                  zmax = MAX( 0.1, xlimphy(ji,jj,jk) )
+                  zpislopen = zpislopead(ji,jj,jk)          &
+                  &         * trn(ji,jj,jk,jpnch) / ( rtrn + trn(ji,jj,jk,jpphy) * 12.)         &
+                  &         / ( prmax(ji,jj,jk) * rjjss * zmax + rtrn )
+
+                  zprbiochl = prmax(ji,jj,jk) * (  1.- EXP( -zpislopen * zetot2 )  )
+
+                  zprorca(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * trn(ji,jj,jk,jpphy) * rfact2
+
+                  zpronew(ji,jj,jk) = zprorca(ji,jj,jk) * xnanono3(ji,jj,jk)    &
                   &             / ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) + rtrn )
-               proreg(ji,jj,jk) = prorca(ji,jj,jk) - pronew(ji,jj,jk)
-
-               zprod = rjjss * prorca(ji,jj,jk) * zprbiochl * trn(ji,jj,jk,jpphy)   &
-                  &                             * MAX( 0.1, xlimphy(ji,jj,jk) )
-
-               prorca5(ji,jj,jk) = (fecnm)**2 * zprod / chlcnm            &
-                  &              / (  zpislopead(ji,jj,jk) * zetot2(ji,jj,jk) * trn(ji,jj,jk,jpnfe) + rtrn  )
-
-               prorca6(ji,jj,jk) = chlcnm * 144. * zprod                  &
-                  &              / (  zpislopead(ji,jj,jk) * zetot2(ji,jj,jk) * trn(ji,jj,jk,jpnch) + rtrn  )
-
-            END DO
-         END DO
-      END DO
-
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-
+                  zprod = rjjss * zprorca(ji,jj,jk) * zprbiochl * trn(ji,jj,jk,jpphy) *zmax
+
+                  zprorca5(ji,jj,jk) = (fecnm)**2 * zprod / chlcnm            &
+                  &              / (  zpislopead(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpnfe) + rtrn  )
+
+                  zprorca6(ji,jj,jk) = chlcnm * 144. * zprod                  &
+                  &              / (  zpislopead(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpnch) + rtrn  )
+               ENDIF
+            END DO
+         END DO
+      END DO
+
+!CDIR NOVERRCHK
+      DO jk = 1, jpkm1
+!CDIR NOVERRCHK
+         DO jj = 1, jpj
+!CDIR NOVERRCHK
+            DO ji = 1, jpi
+               IF( etot(ji,jj,jk) > 1.E-3 ) THEN
 !       Computation of the various production terms for diatoms
 !       -------------------------------------------------------
-               zpislope2n(ji,jj,jk) = zpislopead2(ji,jj,jk) * trn(ji,jj,jk,jpdch)        &
-                  &                 / ( rtrn + trn(ji,jj,jk,jpdia) * 12.)        &
-                  &                 / ( prmax(ji,jj,jk) * rjjss * MAX( 0.1, xlimdia(ji,jj,jk) ) + rtrn )
-
-               zprdiachl = prmax(ji,jj,jk) * (  1.- EXP( -zetot2(ji,jj,jk) * zpislope2n(ji,jj,jk) )  )
-
-               prorca2(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * trn(ji,jj,jk,jpdia) * rfact2
-
-               pronew2(ji,jj,jk) = prorca2(ji,jj,jk) * xdiatno3(ji,jj,jk)     &
+                  zetot2 = ediat(ji,jj,jk) * 24. / ( strn(ji,jj) + rtrn )
+                  zmax = MAX( 0.1, xlimdia(ji,jj,jk) )
+                  zpislope2n = zpislopead2(ji,jj,jk) * trn(ji,jj,jk,jpdch)        &
+                  &           / ( rtrn + trn(ji,jj,jk,jpdia) * 12.)        &
+                  &           / ( prmax(ji,jj,jk) * rjjss * zmax + rtrn )
+
+                  zprdiachl = prmax(ji,jj,jk) * (  1.- EXP( -zetot2 * zpislope2n )  )
+
+                  zprorca2(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * trn(ji,jj,jk,jpdia) * rfact2
+
+                  zpronew2(ji,jj,jk) = zprorca2(ji,jj,jk) * xdiatno3(ji,jj,jk)     &
                   &              / ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) + rtrn )
-               proreg2(ji,jj,jk) = prorca2(ji,jj,jk) - pronew2(ji,jj,jk)
-               prorca3(ji,jj,jk) = prorca2(ji,jj,jk) * zsopt(ji,jj,jk)
-
-               zprod=rjjss * prorca2(ji,jj,jk) * zprdiachl * trn(ji,jj,jk,jpdia) * MAX( 0.1, xlimdia(ji,jj,jk) )
-
-               prorca4(ji,jj,jk) = (fecdm)**2 * zprod / chlcdm                   &
-                  &              / ( zpislopead2(ji,jj,jk) * zetot2(ji,jj,jk) * trn(ji,jj,jk,jpdfe) + rtrn )
-
-               prorca7(ji,jj,jk) = chlcdm * 144. * zprod       &
-                  &              / ( zpislopead2(ji,jj,jk) * zetot2(ji,jj,jk) * trn(ji,jj,jk,jpdch) + rtrn )
-
+
+                  zprod = rjjss * zprorca2(ji,jj,jk) * zprdiachl * trn(ji,jj,jk,jpdia) * zmax
+
+                  zprorca4(ji,jj,jk) = (fecdm)**2 * zprod / chlcdm                   &
+                  &              / ( zpislopead2(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpdfe) + rtrn )
+
+                  zprorca7(ji,jj,jk) = chlcdm * 144. * zprod       &
+                  &              / ( zpislopead2(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpdch) + rtrn )
+
+               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
+              zproreg  = zprorca(ji,jj,jk) - zpronew(ji,jj,jk)
+              zproreg2 = zprorca2(ji,jj,jk) - zpronew2(ji,jj,jk)
+              tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - zprorca(ji,jj,jk) - zprorca2(ji,jj,jk)
+              tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) - zpronew(ji,jj,jk) - zpronew2(ji,jj,jk)
+              tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) - zproreg - zproreg2
+              tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) + zprorca(ji,jj,jk) * zexcret
+              tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) + zprorca6(ji,jj,jk) * zexcret
+              tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) + zprorca5(ji,jj,jk) * zexcret
+              tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) + zprorca2(ji,jj,jk) * zexcret2
+              tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) + zprorca7(ji,jj,jk) * zexcret2
+              tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) + zprorca4(ji,jj,jk) * zexcret2
+              tra(ji,jj,jk,jpbsi) = tra(ji,jj,jk,jpbsi) + zprorca2(ji,jj,jk) * zysopt(ji,jj,jk) * zexcret2
+              tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + &
+              &                     excret2 * zprorca2(ji,jj,jk) + excret * zprorca(ji,jj,jk)
+              tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + o2ut * ( zproreg + zproreg2) &
+              &                    + ( o2ut + o2nit ) * ( zpronew(ji,jj,jk) + zpronew2(ji,jj,jk) )
+              tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) &
+              &                     - zexcret * zprorca5(ji,jj,jk) - zexcret2 * zprorca4(ji,jj,jk)
+              tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) &
+              &                     - zexcret2 * zprorca2(ji,jj,jk) * zysopt(ji,jj,jk)
+              tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprorca(ji,jj,jk) - zprorca2(ji,jj,jk)
+              tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) &
+              &                    + rno3 * ( zpronew(ji,jj,jk) + zpronew2(ji,jj,jk) )
+          END DO
+        END DO
+     END DO
+
+
+     DO jk = 1, jpkm1
+        DO jj = 1, jpj
+          DO ji = 1, jpi
+             qcumul(2) = qcumul(2) + ( zprorca(ji,jj,jk) + zprorca2(ji,jj,jk) )  &
+             &                       * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) * tmask_i(ji,jj)
+          END DO
+        END DO
+      END DO
+
+      IF( MOD( kt, nspyr ) == 0 ) THEN
+        WRITE(numout,*) 'Total PP :'
+        WRITE(numout,*) '-------------------- : ',qcumul(2)*12./1E12
+        WRITE(numout,*) '(GtC/an)'
+        qcumul(2) = 0.
+      ENDIF
+
+#if defined key_trc_dia3d
+      zrfact2 = 1.e3 * rfact2r
+!   Supplementary diagnostics
+!   -------------------------
+      trc3d(:,:,:,4)  = etot(:,:,:)
+      trc3d(:,:,:,5)  = zprorca(:,:,:)  * zrfact2
+      trc3d(:,:,:,6)  = zprorca2(:,:,:) * zrfact2
+      trc3d(:,:,:,7)  = zpronew(:,:,:)  * zrfact2
+      trc3d(:,:,:,8)  = zpronew2(:,:,:) * zrfact2
+      trc3d(:,:,:,9)  = zprorca2(:,:,:) * zysopt(:,:,:) * zrfact2
+      trc3d(:,:,:,10) = zprorca4(:,:,:) * zrfact2
+#if ! defined key_kriest
+      trc3d(:,:,:,11) = zprorca5(:,:,:) * zrfact2
+#endif
+#endif
+
+       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+         WRITE(charout, FMT="('prod')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+       ENDIF
+
    END SUBROUTINE p4z_prod
 
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zrem.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zrem.F90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zrem.F90	(revision 858)
@@ -17,4 +17,5 @@
    USE trp_trc         ! 
    USE sms             ! 
+   USE prtctl_trc
 
    IMPLICIT NONE
@@ -45,7 +46,11 @@
       REAL(wp) ::   zkeq  , zfeequi, zsiremin
       REAL(wp) ::   zsatur, zsatur2, znusil
-      REAL(wp) ::   zlamfac, zstep
+      REAL(wp) ::   zbactfer, zorem, zorem2, zofer, zofer2
+      REAL(wp) ::   zosil, zdenom, zdenom1, zdenom2, zscave, zaggdfe
+      REAL(wp) ::   zlamfac, zstep, zonitr
       REAL(wp), DIMENSION(jpi,jpj)     ::   ztempbac
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdepbac, zfesatur
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdepbac, zfesatur, zolimi
+      CHARACTER (len=25) :: charout
+
       !!---------------------------------------------------------------------
 
@@ -85,4 +90,5 @@
       nitrfac(:,:,:) = MIN( 1., nitrfac(:,:,:) )
 
+
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -102,10 +108,10 @@
 !     Ammonification in oxic waters with oxygen consumption
 !     -----------------------------------------------------
-               olimi(ji,jj,jk) = MIN(  ( trn(ji,jj,jk,jpoxy) - rtrn ) / o2ut,                     &
+               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
 !     -------------------------------------------------------
-               denitr(ji,jj,jk) = MIN(  ( trn(ji,jj,jk,jpno3) - rtrn ) / rdenit,           &
+               denitr(ji,jj,jk) = MIN(  ( trn(ji,jj,jk,jpno3) - rtrn ) / rdenit,   &
                   &                     zremik * nitrfac(ji,jj,jk) * trn(ji,jj,jk,jpdoc)  )
             END DO
@@ -113,6 +119,6 @@
       END DO
 
-      olimi (:,:,:) = MAX( 0.e0, olimi (:,:,:) )
-      denitr(:,:,:) = MAX( 0.e0, denitr(:,:,:) )
+      zolimi (:,:,:) = MAX( 0.e0, zolimi (:,:,:) )
+      denitr (:,:,:) = MAX( 0.e0, denitr (:,:,:) )
 
       DO jk = 1, jpkm1
@@ -123,12 +129,28 @@
 !    below 2 umol/L. Inhibited at strong light 
 !    ----------------------------------------------------------
-               onitr(ji,jj,jk) = nitrif * zstep * trn(ji,jj,jk,jpnh4) / ( 1.+ emoy(ji,jj,jk) )     &
-# if defined key_off_degrad
-                  &            * facvol(ji,jj,jk)              &
-# endif
-                  &            * ( 1.- nitrfac(ji,jj,jk) )
-            END DO
-         END DO
-      END DO
+               zonitr  = nitrif * zstep * trn(ji,jj,jk,jpnh4) / ( 1.+ emoy(ji,jj,jk) )     &
+# if defined key_off_degrad
+                  &      * facvol(ji,jj,jk)              &
+# 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
+
+            END DO
+         END DO
+      END DO
+
+       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+         WRITE(charout, FMT="('rem1')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+       ENDIF
 
       DO jk = 1, jpkm1
@@ -141,12 +163,25 @@
 !    significant
 !    ----------------------------------------------------------
-               xbactfer(ji,jj,jk) = 15.e-6 * rfact2 * 4.* 0.4 * prmax(ji,jj,jk)           &
+               zbactfer = 15.e-6 * rfact2 * 4.* 0.4 * prmax(ji,jj,jk)           &
                   &               * ( xlimphy(ji,jj,jk) * zdepbac(ji,jj,jk))**2           &
                   &                  / ( xkgraz2 + zdepbac(ji,jj,jk) )                    &
                   &                  * ( 0.5 + SIGN( 0.5, trn(ji,jj,jk,jpfer) -2.e-11 )  )
 
-            END DO
-         END DO
-      END DO
+               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zbactfer
+#if defined key_kriest
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zbactfer
+#else
+               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zbactfer
+#endif
+
+            END DO
+         END DO
+      END DO
+
+       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+         WRITE(charout, FMT="('rem2')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+       ENDIF
 
       DO jk = 1, jpkm1
@@ -167,15 +202,38 @@
 !    means a disaggregation constant about 0.5 the value in oxic zones
 !    -----------------------------------------------------------------
-               orem (ji,jj,jk) = zremip * trn(ji,jj,jk,jppoc)
-               ofer (ji,jj,jk) = zremip * trn(ji,jj,jk,jpsfe)
+               zorem  = zremip * trn(ji,jj,jk,jppoc)
+               zofer  = zremip * trn(ji,jj,jk,jpsfe)
 #if ! defined key_kriest
-               orem2(ji,jj,jk) = zremip * trn(ji,jj,jk,jpgoc)
-               ofer2(ji,jj,jk) = zremip * trn(ji,jj,jk,jpbfe)
-#else
-               orem2(ji,jj,jk) = zremip * trn(ji,jj,jk,jpnum)
-#endif
-            END DO
-         END DO
-      END DO
+               zorem2 = zremip * trn(ji,jj,jk,jpgoc)
+               zofer2 = zremip * trn(ji,jj,jk,jpbfe)
+#else
+               zorem2 = zremip * trn(ji,jj,jk,jpnum)
+#endif
+
+!  Update the appropriate tracers trends
+!  -------------------------------------
+
+               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zorem
+               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zofer
+#if defined key_kriest
+               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zorem
+               tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) - zorem2
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zofer
+#else
+               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zorem2 - zorem
+               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) - zorem2
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zofer2 - zofer
+               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) - zofer2
+#endif
+
+            END DO
+         END DO
+      END DO
+
+       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+         WRITE(charout, FMT="('rem3')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+       ENDIF
 
       DO jk = 1, jpkm1
@@ -194,5 +252,9 @@
                zsiremin = xsirem * zstep * znusil
 #    endif
-               osil(ji,jj,jk) = zsiremin * trn(ji,jj,jk,jpdsi)
+               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
@@ -200,9 +262,32 @@
       END DO
 
+      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+         WRITE(charout, FMT="('rem4')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+       ENDIF
+
       zfesatur(:,:,:) = 0.6e-9
-
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
+!CDIR NOVERRCHK
+      DO jk = 1, jpkm1
+!CDIR NOVERRCHK
+         DO jj = 1, jpj
+!CDIR NOVERRCHK
+            DO ji = 1, jpi
+!
+!      Compute de different ratios for scavenging of iron
+!      --------------------------------------------------
+
+#if  defined key_kriest
+               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)  &
+           &            + 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
@@ -222,7 +307,7 @@
 #endif
 # if defined key_off_degrad
-               xscave(ji,jj,jk) = zfeequi * zlam1b * zstep  * facvol(ji,jj,jk)
+               zscave = zfeequi * zlam1b * zstep  * facvol(ji,jj,jk)
 # else
-               xscave(ji,jj,jk) = zfeequi * zlam1b * zstep
+               zscave = zfeequi * zlam1b * zstep
 # endif
 
@@ -236,22 +321,58 @@
                zlam1b = (  80.* ( trn(ji,jj,jk,jpdoc) + 35.e-6 )                           &
                   &     + 698.*   trn(ji,jj,jk,jppoc) + 1.05e4 * trn(ji,jj,jk,jpgoc)  )                    &
-                  &   * zdiss(ji,jj,jk) + 1E-4 * (1.-zlamfac)                &
+                  &   * xdiss(ji,jj,jk) + 1E-4 * (1.-zlamfac)                &
                   &   + xlam1 * MAX( 0.e0, ( trn(ji,jj,jk,jpfer) * 1.e9 - 1.)  )
 #else
                zlam1b = (  80.* (trn(ji,jj,jk,jpdoc) + 35E-6)           &
                   &     + 698.*  trn(ji,jj,jk,jppoc)  )                    &
-                  &   * zdiss(ji,jj,jk) + 1E-4 * (1.-zlamfac)           &
+                  &   * xdiss(ji,jj,jk) + 1E-4 * (1.-zlamfac)           &
                   &   + xlam1 * MAX( 0.e0, ( trn(ji,jj,jk,jpfer) * 1.e9 - 1.)  )
 #endif
 
 # if defined key_off_degrad
-               xaggdfe(ji,jj,jk) = zlam1b * zstep * 0.5 * ( trn(ji,jj,jk,jpfer) - zfeequi ) * facvol(ji,jj,jk)
+               zaggdfe = zlam1b * zstep * 0.5 * ( trn(ji,jj,jk,jpfer) - zfeequi ) * facvol(ji,jj,jk)
 # else
-               xaggdfe(ji,jj,jk) = zlam1b * zstep * 0.5 * ( trn(ji,jj,jk,jpfer) - zfeequi )
-# endif
+               zaggdfe = zlam1b * zstep * 0.5 * ( trn(ji,jj,jk,jpfer) - zfeequi )
+# endif
+
+               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfe
+
+#if defined key_kriest
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1
+#else
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1
+               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zscave * zdenom2
+#endif
+
             END DO
          END DO
       END DO
       !
+
+       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+         WRITE(charout, FMT="('rem5')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+       ENDIF
+
+!     Update the arrays TRA which contain the biological sources and sinks
+!     --------------------------------------------------------------------
+
+      DO jk = 1, jpkm1
+         tra(:,:,jk,jppo4) = tra(:,:,jk,jppo4) + zolimi(:,:,jk) + denitr(:,:,jk)
+         tra(:,:,jk,jpnh4) = tra(:,:,jk,jpnh4) + zolimi(:,:,jk) + denitr(:,:,jk)
+         tra(:,:,jk,jpno3) = tra(:,:,jk,jpno3) - denitr(:,:,jk) * rdenit
+         tra(:,:,jk,jpdoc) = tra(:,:,jk,jpdoc) - zolimi(:,:,jk) - denitr(:,:,jk)
+         tra(:,:,jk,jpoxy) = tra(:,:,jk,jpoxy) - zolimi(:,:,jk) * o2ut
+         tra(:,:,jk,jpdic) = tra(:,:,jk,jpdic) + zolimi(:,:,jk) + denitr(:,:,jk)
+         tra(:,:,jk,jptal) = tra(:,:,jk,jptal) + denitr(:,:,jk) * rno3 * rdenit
+     END DO
+
+       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+         WRITE(charout, FMT="('rem6')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+       ENDIF
+
    END SUBROUTINE p4z_rem
 
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zsed.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zsed.F90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zsed.F90	(revision 858)
@@ -17,4 +17,6 @@
    USE sms
    USE lib_mpp
+   USE prtctl_trc
+
 
    IMPLICIT NONE
@@ -50,4 +52,5 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   zsidep
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   znitrpot, zirondep
+      CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
 
@@ -221,5 +224,5 @@
             DO ji = 2, jpim1
                zdenitot = zdenitot + denitr(ji,jj,jk) * rdenit * e1t(ji,jj) * e2t(ji,jj)   &
-                  &    *fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) * znegtr(ji,jj,jk)
+                  &    *fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) * xnegtr(ji,jj,jk)
             END DO
          END DO
@@ -231,6 +234,9 @@
       ! -------------------------------------------------------------
 
+!CDIR NOVERRCHK
       DO jk = 1, jpk
-         DO jj = 1, jpj
+!CDIR NOVERRCHK
+         DO jj = 1, jpj
+!CDIR NOVERRCHK
             DO ji = 1, jpi
                zlim = ( 1.- xnanono3(ji,jj,jk) - xnanonh4(ji,jj,jk) )
@@ -265,5 +271,6 @@
             DO ji = 1, jpi
 # if ! defined key_cfg_1d && ( 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) * zdenitot / znitrpottot
+               zfact = znitrpot(ji,jj,jk) * 1.e-7
 # else
                zfact = znitrpot(ji,jj,jk) * 1.e-7
@@ -285,4 +292,10 @@
 # endif
       !
+       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+         WRITE(charout, FMT="('sed ')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc(tab4d=trn, mask=tmask, clinfo=ctrcnm)
+       ENDIF
+
    END SUBROUTINE p4z_sed
 
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zsink.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zsink.F90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zsink.F90	(revision 858)
@@ -18,4 +18,6 @@
    USE sms
    USE p4zsink2        !
+   USE prtctl_trc
+
 
    IMPLICIT NONE
@@ -46,4 +48,5 @@
       INTEGER  ::   iksed
       REAL(wp) ::   zagg1, zagg2, zagg3, zagg4
+      REAL(wp) ::   zagg , zaggfe, zaggdoc, zaggdoc2
       REAL(wp) ::   zfact, zstep, zwsmax
 #if defined key_trc_dia3d
@@ -53,4 +56,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsinkfer, zsinkfer2
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsinkcal, zsinksil
+      CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
 
@@ -120,5 +124,5 @@
             DO ji = 1, jpi
 
-               zfact = zstep * zdiss(ji,jj,jk)
+               zfact = zstep * xdiss(ji,jj,jk)
 
 !    Part I : Coagulation dependent on turbulence
@@ -155,21 +159,27 @@
 # endif
 
-               xagg  (ji,jj,jk) = zagg1 + zagg2 + zagg3 + zagg4
-               xaggfe(ji,jj,jk) = xagg(ji,jj,jk) * trn(ji,jj,jk,jpsfe) / ( trn(ji,jj,jk,jppoc) + rtrn )
+               zagg   = zagg1 + zagg2 + zagg3 + zagg4
+               zaggfe = zagg * trn(ji,jj,jk,jpsfe) / ( trn(ji,jj,jk,jppoc) + rtrn )
 
 !     Aggregation of DOC to small particles
 !     --------------------------------------
 
-               xaggdoc(ji,jj,jk) = ( 80.* trn(ji,jj,jk,jpdoc) + 698. * trn(ji,jj,jk,jppoc) )       &
-# if defined key_off_degrad
-                  &              * facvol(ji,jj,jk)                           &
-# endif
-                  &              * zfact * trn(ji,jj,jk,jpdoc)
-
-               xaggdoc2(ji,jj,jk) = 1.05e4 * zfact * trn(ji,jj,jk,jpgoc)   &
-# if defined key_off_degrad
-                  &               * facvol(ji,jj,jk)                            &
+               zaggdoc = ( 80.* trn(ji,jj,jk,jpdoc) + 698. * trn(ji,jj,jk,jppoc) )       &
+# if defined key_off_degrad
+                  &      * facvol(ji,jj,jk)                           &
+# endif
+                  &      * zfact * trn(ji,jj,jk,jpdoc)
+
+               zaggdoc2 = 1.05e4 * zfact * trn(ji,jj,jk,jpgoc)   &
+# if defined key_off_degrad
+                  &        * facvol(ji,jj,jk)                            &
 # endif      
-                  &               * trn(ji,jj,jk,jpdoc)
+                  &        * trn(ji,jj,jk,jpdoc)
+
+               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zagg + zaggdoc
+               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zagg + zaggdoc2
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zaggfe
+               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zaggfe
+               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc2
 
             END DO
@@ -187,4 +197,10 @@
 # endif
       !
+       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+         WRITE(charout, FMT="('sink')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+       ENDIF
+
    END SUBROUTINE p4z_sink
 
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zsink2.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zsink2.F90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zsink2.F90	(revision 858)
@@ -47,11 +47,11 @@
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   sinktemp   ! ???
       !!
-      INTEGER  ::   ji, jj, jk
+      INTEGER  ::   ji, jj, jk, jnt
       REAL(wp) ::   zigma,zew,zstep,zign
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztraz, zakz
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zkz  , wstmp2
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zkz  , zwstmp2
       !!---------------------------------------------------------------------
 
-      zstep  = rfact2
+      zstep  = rfact2 / 2.
 
       ztraz(:,:,:) = 0.e0
@@ -61,70 +61,92 @@
       DO jk = 1, jpkm1
 # if defined key_off_degrad
-         wstmp2(:,:,jk+1)=-wstmp(:,:,jk)/rjjss*tmask(:,:,jk+1)*facvol(:,:,jk)
+         zwstmp2(:,:,jk+1)=-wstmp(:,:,jk)/rjjss*tmask(:,:,jk+1)*facvol(:,:,jk)
 # else
-         wstmp2(:,:,jk+1)=-wstmp(:,:,jk)/rjjss*tmask(:,:,jk+1)
+         zwstmp2(:,:,jk+1)=-wstmp(:,:,jk)/rjjss*tmask(:,:,jk+1)
 
 # endif
       END DO
  
-      wstmp2(:,:,1) = 0.e0
+      zwstmp2(:,:,1) = 0.e0
 !
 ! Vertical advective flux
 !-------------------------------
-! ... first guess of the slopes
-!   ... interior values
-      DO jk = 2, jpkm1
-         ztraz(:,:,jk) = (trn(:,:,jk-1,jn) - trn(:,:,jk,jn)) *tmask(:,:,jk)
-      END DO
+
+      DO jnt = 1, 2
+
+! ... first guess of the slopes interior values
+
+         DO jk = 2, jpkm1
+            ztraz(:,:,jk) = (trn(:,:,jk-1,jn) - trn(:,:,jk,jn)) *tmask(:,:,jk)
+         END DO
+
+         ztraz(:,:,1  ) = 0.0
+         ztraz(:,:,jpk) = 0.0
 !
 ! slopes
-      DO jk=2,jpkm1
-         DO jj = 1,jpj
-            DO ji = 1, jpi
-               zign = 0.5*(sign(1.,ztraz(ji,jj,jk)*ztraz(ji,jj,jk+1))+1)
-               zakz(ji,jj,jk) = 0.5*(ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign
-            END DO
-          END DO
-        END DO        
-!
-! Slopes limitation
-      DO jk = 2, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zakz(ji,jj,jk) = sign(1.,zakz(ji,jj,jk)) *        &
-                  &             min(abs(zakz(ji,jj,jk)),         &
-                  &             2.*abs(ztraz(ji,jj,jk+1)),       &
-                  &             2.*abs(ztraz(ji,jj,jk)))
+         DO jk=2,jpkm1
+            DO jj = 1,jpj
+               DO ji = 1, jpi
+                  zign = 0.25 + SIGN( 0.25,ztraz(ji,jj,jk)*ztraz(ji,jj,jk+1) )
+                  zakz(ji,jj,jk) = (ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign
+               END DO
             END DO
          END DO
-      END DO        
-
-! vertical advective flux
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj      
-            DO ji = 1, jpi    
-               zigma = wstmp2(ji,jj,jk+1)*zstep/fse3w(ji,jj,jk+1)
-               zew   = wstmp2(ji,jj,jk+1)
-               sinktemp(ji,jj,jk+1) = -zew*(trn(ji,jj,jk,jn)                &
-                  &                 -0.5*(1+zigma)*zakz(ji,jj,jk))*zstep
+         !
+         ! Slopes limitation
+         DO jk = 2, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zakz(ji,jj,jk) = SIGN(1.,zakz(ji,jj,jk)) *        &
+                     &             MIN(ABS(zakz(ji,jj,jk)),         &
+                     &             2.*ABS(ztraz(ji,jj,jk+1)),       &
+                     &             2.*ABS(ztraz(ji,jj,jk)))
+               END DO
             END DO
          END DO
-      END DO 
-!
-! Boundary conditions
-      sinktemp(:,:,1  ) = 0.e0
-      sinktemp(:,:,jpk) = 0.e0
+         
+         ! vertical advective flux
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj      
+               DO ji = 1, jpi    
+                  zigma = zwstmp2(ji,jj,jk+1)*zstep/fse3w(ji,jj,jk+1)
+                  zew   = zwstmp2(ji,jj,jk+1)
+                  sinktemp(ji,jj,jk+1) = -zew*(trn(ji,jj,jk,jn)                &
+                     &                 -0.5*(1+zigma)*zakz(ji,jj,jk))*zstep
+               END DO
+            END DO
+         END DO
+         !
+         ! Boundary conditions
+         sinktemp(:,:,1  ) = 0.e0
+         sinktemp(:,:,jpk) = 0.e0
+         
+         DO jk=1,jpkm1
+            DO jj = 1,jpj
+               DO ji = 1, jpi
+                  trn(ji,jj,jk,jn) = trn(ji,jj,jk,jn)               &
+                     &        + (sinktemp(ji,jj,jk)-sinktemp(ji,jj,jk+1))     &
+                     &        /fse3t(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+
+      ENDDO
 
       DO jk=1,jpkm1
          DO jj = 1,jpj
             DO ji = 1, jpi
-               trn(ji,jj,jk,jn) = trn(ji,jj,jk,jn)               &
-                  &        + (sinktemp(ji,jj,jk)-sinktemp(ji,jj,jk+1))     &
-                  &        /fse3t(ji,jj,jk)
-            END DO
-         END DO
-      END DO
+               ! 
+               trb(ji,jj,jk,jn) = trb(ji,jj,jk,jn)      &
+               &        + 2.*(sinktemp(ji,jj,jk)-sinktemp(ji,jj,jk+1)) &
+               &        /fse3t(ji,jj,jk)
+!
+          ENDDO
+        ENDDO
+      ENDDO
+!
+        trn(:,:,:,jn)   = trb(:,:,:,jn)
+        sinktemp(:,:,:) = 2. * sinktemp(:,:,:)
 
-      trb(:,:,:,jn) = trn(:,:,:,jn)
       !
    END SUBROUTINE p4z_sink2
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zsink_kriest.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zsink_kriest.F90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/p4zsink_kriest.F90	(revision 858)
@@ -20,4 +20,5 @@
    USE sms
    USE p4zsink2
+   USE prtctl_trc
 
    IMPLICIT NONE
@@ -45,11 +46,12 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk
-      INTEGER  ::   iksed
-      REAL(wp) ::   zagg1, zagg2, zagg3, zagg4, zagg5, zaggsi, zaggsh
-      REAL(wp) ::   znum , zeps, zfm, zgm, zsm
-      REAL(wp) ::   zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5
-      REAL(wp) ::   zval1, zval2, zval3, zval4
-      REAL(wp) ::   zstep
+      INTEGER  :: ji, jj, jk
+      INTEGER  :: iksed
+      REAL(wp) :: zagg1, zagg2, zagg3, zagg4, zagg5, zaggsi, zaggsh
+      REAL(wp) :: zagg , zaggdoc, znumdoc
+      REAL(wp) :: znum , zeps, zfm, zgm, zsm
+      REAL(wp) :: zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5
+      REAL(wp) :: zval1, zval2, zval3, zval4
+      REAL(wp) :: zstep
 #if defined key_trc_dia3d
       REAL(wp) ::   zrfact2
@@ -59,4 +61,6 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   sinkfer
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   sinkcal, sinksil
+      CHARACTER (len=25) :: charout
+
       !!---------------------------------------------------------------------
 
@@ -184,5 +188,5 @@
                      &    )
 
-                  zaggsh = ( zagg1 + zagg2 + zagg3 ) * rfact2 * zdiss(ji,jj,jk) / 1000.
+                  zaggsh = ( zagg1 + zagg2 + zagg3 ) * rfact2 * xdiss(ji,jj,jk) / 1000.
 
 !    Aggregation of small into large particles
@@ -213,15 +217,20 @@
                   zaggsi = ( zagg4 + zagg5 ) * zstep / 10.
 
-                  xagg(ji,jj,jk) = 0.5 * xkr_stick * ( zaggsh + zaggsi )
+                  zagg = 0.5 * xkr_stick * ( zaggsh + zaggsi )
 
 !     Aggregation of DOC to small particles
 !     --------------------------------------
 
-                  xaggdoc(ji,jj,jk) = (     0.4 * trn(ji,jj,jk,jpdoc)               &
-                     &                 + 1018.  * trn(ji,jj,jk,jppoc)  ) * zstep    &
-# if defined key_off_degrad
-                     &              * facvol(ji,jj,jk)                              &
-# endif
-                     &              * zdiss(ji,jj,jk) * trn(ji,jj,jk,jpdoc)
+                  zaggdoc = ( 0.4 * trn(ji,jj,jk,jpdoc)               &
+                     &        + 1018.  * trn(ji,jj,jk,jppoc)  ) * zstep    &
+# if defined key_off_degrad
+                     &        * facvol(ji,jj,jk)                              &
+# endif
+                     &        * xdiss(ji,jj,jk) * trn(ji,jj,jk,jpdoc)
+                      
+                  znumdoc = trn(ji,jj,jk,jpnum) / ( trn(ji,jj,jk,jppoc) + rtrn )
+                  tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zaggdoc
+                  tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zaggdoc * znumdoc - zagg
+                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc
 
                ENDIF
@@ -246,4 +255,9 @@
 #    endif
       !
+       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+         WRITE(charout, FMT="('sink')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+       ENDIF
    END SUBROUTINE p4z_sink_kriest
 
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/sms_pisces.h90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/sms_pisces.h90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/sms_pisces.h90	(revision 858)
@@ -15,30 +15,15 @@
    !! Variable for chemistry of the CO2 cycle
    !! ---------------------------------------------------------------------
-   REAL(wp) ::   atcco2, atcox
+   REAL(wp) ::   atcco2
    !
    REAL(wp), DIMENSION(jpi,jpj,jpk) ::   akb3, ak13, ak23, aksp, akw3             !: ???
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   akp13, akp23, akp33, aksi3, aks3, akf3   !: ???
    REAL(wp), DIMENSION(jpi,jpj,jpk) ::   hi, borat                                !: ???
+   REAL, DIMENSION(2) :: qcumul
 
    !!----------------------------------------------------------------------
    !!  Variable for chemistry of the CO2 cycle
    !! ---------------------------------------------------------------------
-   REAL(wp), DIMENSION(10) ::   devk1, devk2, devk3, devk4, devk5     
    !
-   REAL(wp) ::   akcc1, akcc2, akcc3, akcc4,akcc5, akcc6, akcc7, akcc8, akcc9     !: ???
-   REAL(wp) ::   bor1, bor2, c00, c01, c02, c03, c04, c05, c10, c11               !: ???
-   REAL(wp) ::   c12, c13, c20, c21, c22, c23, cb0, cb1, cb2, cb3                 !: ???
-   REAL(wp) ::   cb4, cb5, cb6, cb7, cb8, cb9, cb10, cb11, c14                    !: ???
-   REAL(wp) ::   cw3, cw4, cw5, cw6, cw0, cw1, cw2, ox0, ox1, ox2, ox3, ox4,ox5   !: ???
-   REAL(wp) ::   salchl, rgas, oxyco, ca0, ca1, ca2, ca3, ca4, ca5, ca6           !: ???
-   REAL(wp) ::   cp10, cp11, cp12, cp13, cp14, cp15, cp16, cp20, cp21             !: ???
-   REAL(wp) ::   cp22, cp23, cp24, cp25, cp26, cp30, cp31, cp32, cp33             !: ???
-   REAL(wp) ::   cp34, cp35, cs10, cs11, cs12, cs13, cs14, cs15, cs16             !: ???
-   REAL(wp) ::   cs17, cs18, cs19, cs20, cs21                                     !: ???
-   REAL(wp) ::   st1, st2, ft1, ft2, ks0, ks1, ks2, ks3, ks4, ks5                 !: ???
-   REAL(wp) ::   ks6, ks7, ks8, ks9, ks10, ks11, ks12, kf0, kf1                   !: ???
-   REAL(wp) ::   kf2, kf3, kf4
-   !
-   REAL(wp), DIMENSION(jpi,jpj,3) ::   chemc   !: ???
+   REAL(wp), DIMENSION(jpi,jpj,2) ::   chemc   !: ???
 
    !!----------------------------------------------------------------------
@@ -68,5 +53,5 @@
    REAL(wp) ::   caco3r, kdca, nca, part, rno3, o2ut, po4r               !: ???
    REAL(wp) ::   sco2, dispo0, conc0,sumdepsi,rivalkinput,sedfeinput     !: ???
-   REAL(wp) ::   calcon, rivpo4input,nitdepinput,oxymin,spocri           !: ???
+   REAL(wp) ::   rivpo4input,nitdepinput,oxymin,spocri           !: ???
    REAL(wp) ::   nitrif,rdenit,o2nit,concnnh4,concdnh4                   !: ???
    REAL(wp) ::   pislope,excret,wsbio,resrat,mprat,wchl,wchld            !: ???
@@ -83,6 +68,7 @@
    !!---------------------------------------------
    REAL(wp), DIMENSION(3,61)        ::   xkrgb                 !: ???
-   REAL(wp), DIMENSION(jpi,jpj)     ::   zmeu                  !: ???
+   REAL(wp), DIMENSION(jpi,jpj)     ::   heup                  !: ???
    REAL(wp), DIMENSION(jpi,jpj,jpk) ::   etot, etot3, emoy     !: ???
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   enano, ediat
 
 !!----------------------------------------------------------
@@ -92,28 +78,6 @@
    !
    REAL(wp), DIMENSION(jpi,jpj,jpk) ::   prmax, tgfunc, tgfunc2                       !: ???
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   prcaca, prorca, prorca2, prorca3             !: ???
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   prorca4, prorca5, prorca6, prorca7           !: ???
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   pronew, pronew2, proreg, proreg2             !: ???
    REAL(wp), DIMENSION(jpi,jpj,jpk) ::   xnanono3, xdiatno3, xnanonh4, xdiatnh4       !: ???
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   xlimphy, xlimdia, concdfe, concnfe, znegtr   !: ???
-
-   !!------------------------------------------
-   !! Sinks for phytoplankton
-   !!------------------------------------------
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   tortp, tortnf, tortnch            !: ???
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   respp, respp2, respnch, respdch   !: ???
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   tortp2, tortdf, tortdch, tortds   !: ???
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   respds, respdf, respnf            !: ???
-
-   !!------------------------------------
-   !!  SMS for zooplankton
-   !!-------------------------------------
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   respz, tortz, grazp, grazpf                   !: ???
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   grazpch, grazm, grazmf, grazsd                !: ???
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   grazsf, grazss, grazsch, grarem               !: ???
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   grafer,respz2,tortz2,grazd, grazz,grazn       !: ???
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   grazpoc,graznf, graznch, grazs, grazf         !: ???
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   grazdch, grazpof, grarem2, grafer2, grapoc2   !: ???
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   grapoc, grazffe, grazfff                      !: ???
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   xlimphy, xlimdia, concdfe, concnfe, xnegtr   !: ???
 
    !!---------------------------------------------
@@ -122,7 +86,5 @@
    REAL(wp) ::   wsbio2
    !
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   xagg, xaggfe, zdiss, xaggdoc, xaggdfe, xbactfer   !: ???
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   xscave, olimi, orem, orem2, ofer, ofer2           !: ???
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   osil, xaggdoc2, nitrfac, xlimbac                  !: ???
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   xdiss, xfracal, nitrfac, xlimbac                  !: ???
    REAL(wp), DIMENSION(jpi,jpj,jpk) ::   wsbio4, wsbio3, wscal                             !: ???
 
@@ -136,5 +98,5 @@
    REAL(wp), DIMENSION(jpi,jpj)     ::   dust, cotdep, nitdep, rivinp   !: ???
    REAL(wp), DIMENSION(jpi,jpj,12)  ::   dustmo                         !: ???
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   onitr, denitr, ironsed         !: ???
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   denitr, ironsed         !: ???
 
 #if defined key_kriest
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/trcctl.pisces.h90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/trcctl.pisces.h90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/trcctl.pisces.h90	(revision 858)
@@ -18,7 +18,7 @@
    ! -----------------------
 #if  defined key_kriest
-      IF( jptra /= 23) THEN 
+      IF( jp_pisces /= 23) THEN 
 #else
-      IF( jptra /= 24) THEN
+      IF( jp_pisces /= 24) THEN
 #endif
           IF (lwp) THEN 
@@ -26,5 +26,5 @@
               WRITE (numout,*) ' =======   ============= '
               WRITE (numout,*)                               &
-              &   ' STOP, change jptra',               & 
+              &   ' STOP, change jp_pisces',               & 
               &   ' in parameter.passivetrc.pisces.h '  
           END IF 
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 857)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 858)
@@ -21,4 +21,6 @@
    USE oce_trc         ! ocean variables
    USE trp_trc         !
+   USE p4zche 
+   USE lbclnk
 
    USE iom
@@ -28,4 +30,5 @@
 
    PUBLIC   trc_ini_pisces   ! called by trcini.F90 module
+
 
 #  include "domzgr_substitute.h90"
@@ -295,252 +298,4 @@
       o2ut   = 140.     / 122.
 
-      !----------------------------------------------------------------------
-      ! Initialize chemical variables 
-      !----------------------------------------------------------------------
-
-      ! set pre-industrial atmospheric [co2] (ppm) and o2/n2 ratio
-      ! ----------------------------------------------------------
-      atcox = 0.20946
-
-      ! Set lower/upper limits for temperature and salinity
-      ! ---------------------------------------------------
-      salchl = 1.e0 / 1.80655
-      calcon = 1.03e-2
-
-      ! Set coefficients for apparent solubility equilibrium of calcite
-      ! Millero et al. 1995 from Mucci 1983
-      ! --------------------------------------------------------------
-      akcc1 = -171.9065
-      akcc2 =   -0.077993
-      akcc3 = 2839.319
-      akcc4 =   71.595
-      akcc5 =   -0.77712
-      akcc6 =    0.0028426
-      akcc7 =  178.34
-      akcc8 =   -0.07711
-      akcc9 =    0.0041249
-
-      ! Set coefficients for seawater pressure correction
-      ! -------------------------------------------------
-      devk1(1) = -25.5
-      devk2(1) =   0.1271
-      devk3(1) =   0.e0
-      devk4(1) =  -3.08E-3
-      devk5(1) =   0.0877E-3
-      !
-      devk1(2) = -15.82
-      devk2(2) =  -0.0219
-      devk3(2) =   0.e0
-      devk4(2) =   1.13E-3
-      devk5(2) =  -0.1475E-3
-      !
-      devk1(3) = -29.48
-      devk2(3) =   0.1622
-      devk3(3) =   2.608E-3
-      devk4(3) =  -2.84E-3
-      devk5(3) =   0.e0
-      !
-      devk1(4) = -14.51
-      devk2(4) =   0.1211
-      devk3(4) =  -0.321E-3
-      devk4(4) =  -2.67E-3
-      devk5(4) =   0.0427E-3
-      !
-      devk1(5) = -23.12
-      devk2(5) =   0.1758
-      devk3(5) =  -2.647E-3
-      devk4(5) =  -5.15E-3
-      devk5(5) =   0.09E-3
-      !
-      devk1(6) = -26.57
-      devk2(6) =   0.2020
-      devk3(6) =  -3.042E-3
-      devk4(6) =  -4.08E-3
-      devk5(6) =   0.0714E-3
-      !
-      devk1(7) = -25.60
-      devk2(7) =   0.2324
-      devk3(7) =  -3.6246E-3
-      devk4(7) =  -5.13E-3
-      devk5(7) =   0.0794E-3
-      !
-      ! For calcite with Edmond and Gieske 1970
-      !     devkst = 0.23
-      !     devks  = 35.4
-      ! Millero 95 takes this depth dependance for calcite
-      devk1(8) = -48.76
-      devk2(8) =   0.5304
-      devk3(8) =   0.e0
-      devk4(8) = -11.76E-3
-      devk5(8) =   0.3692E-3
-      !
-      ! Coefficients for sulfate and fluoride
-      devk1(9) = -18.03
-      devk2(9) =   0.0466
-      devk3(9) =   0.316e-3
-      devk4(9) =  -4.53e-3
-      devk5(9) =   0.09e-3
-
-      devk1(10) = -9.78
-      devk2(10) = -0.0090
-      devk3(10) = -0.942e-3
-      devk4(10) = -3.91e-3
-      devk5(10) =  0.054e-3
-
-
-      ! Set universal gas constants
-      ! ---------------------------
-      rgas  = 83.143
-      oxyco =  1.e0 / 22.4144
-
-      ! Set boron constants
-      ! -------------------
-      bor1 = 0.00023
-      bor2 = 1.e0 / 10.82
-
-      ! Set volumetric solubility constants for co2 in ml/l (Weiss, 1974)
-      ! -----------------------------------------------------------------
-      c00 = -60.2409
-      c01 =  93.4517
-      c02 =  23.3585
-      c03 =   0.023517
-      c04 =  -0.023656
-      c05 =   0.0047036
-      !
-      ca0 = -162.8301
-      ca1 =  218.2968
-      ca2 =   90.9241
-      ca3 =   -1.47696
-      ca4 =    0.025695
-      ca5 =   -0.025225
-      ca6 =    0.0049867
-
-      ! Set coeff. for 1. dissoc. of carbonic acid (Edmond and Gieskes, 1970)
-      ! ---------------------------------------------------------------------
-      c10 = -3670.7
-      c11 =    62.008
-      c12 =    -9.7944
-      c13 =     0.0118
-      c14 =    -0.000116
-
-      ! Set coeff. for 2. dissoc. of carbonic acid (Edmond and Gieskes, 1970)
-      ! ---------------------------------------------------------------------
-      c20 = -1394.7
-      c21 =    -4.777
-      c22 =     0.0184
-      c23 =    -0.000118
-
-      ! Set constants for calculate concentrations for sulfate and fluoride
-      ! sulfates (Morris & Riley 1966)
-      !----------------------------------------------------------------------
-      st1 = 0.14
-      st2 = 1.e0 / 96.062
-
-      ! fluoride
-      ! --------
-      ft1 = 0.000067
-      ft2 = 1.e0 / 18.9984
-
-      ! sulfates (Dickson 1990 change to mol:kg soln, idem OCMIP)
-      !----------------------------------------------------------
-      ks0  =    141.328
-      ks1  =  -4276.1
-      ks2  =    -23.093
-      ks3  = -13856.
-      ks4  =    324.57
-      ks5  =    -47.986
-      ks6  =  35474.
-      ks7  =   -771.54
-      ks8  =    114.723
-      ks9  =  -2698.
-      ks10 =   1776.
-      ks11 =      1.
-      ks12 =     -0.001005
-
-      ! fluorides (Dickson & Riley 1979 change to mol/kg soln)
-      !-------------------------------------------------------
-      kf0 =  -12.641
-      kf1 = 1590.2
-      kf2 =    1.525
-      kf3 =    1.0
-      kf4 =   -0.001005
-
-      ! Set coeff. for 1. dissoc. of boric acid (Edmond and Gieskes, 1970)
-      ! ------------------------------------------------------------------
-      cb0  = -8966.90
-      cb1  = -2890.53
-      cb2  =   -77.942
-      cb3  =     1.728
-      cb4  =    -0.0996
-      cb5  =   148.0248
-      cb6  =   137.1942
-      cb7  =     1.62142
-      cb8  =   -24.4344
-      cb9  =   -25.085
-      cb10 =    -0.2474
-      cb11 =     0.053105
-
-      ! Set coeff. for dissoc. of water (Dickson and Riley, 1979, 
-      !   eq. 7, coefficient cw2 corrected from 0.9415 to 0.09415 
-      !   after pers. commun. to B. Bacastow, 1988)
-      ! ---------------------------------------------------------
-      cw0 = -13847.26
-      cw1 =    148.9652
-      cw2 =    -23.6521
-      cw3 =    118.67
-      cw4 =     -5.977
-      cw5 =      1.0495
-      cw6 =     -0.01615
-
-      ! Set coeff. for dissoc. of phosphate (Millero (1974)
-      ! ---------------------------------------------------
-      cp10 =   115.54
-      cp11 = -4576.752
-      cp12 =   -18.453
-      cp13 =  -106.736
-      cp14 =     0.69171
-      cp15 =    -0.65643
-      cp16 =    -0.01844
-      !
-      cp20 =   172.1033
-      cp21 = -8814.715
-      cp22 =   -27.927
-      cp23 =  -160.340
-      cp24 =     1.3566
-      cp25 =     0.37335
-      cp26 =    -0.05778
-      !
-      cp30 =   -18.126
-      cp31 = -3070.75
-      cp32 =    17.27039
-      cp33 =     2.81197
-      cp34 =   -44.99486
-      cp35 =    -0.09984
-
-      ! Set coeff. for dissoc. of phosphate (Millero (1974)
-      ! ---------------------------------------------------
-      cs10 =   117.385
-      cs11 = -8904.2
-      cs12 =   -19.334
-      cs13 =  -458.79
-      cs14 =     3.5913
-      cs15 =   188.74
-      cs16 =    -1.5998
-      cs17 =   -12.1652
-      cs18 =     0.07871
-      cs19 =     0.e0
-      cs20 =     1.e0
-      cs21 =    -0.001005
-
-
-      ! Set volumetric solubility constants for o2 in ml/l (Weiss, 1970)
-      ! ----------------------------------------------------------------
-      ox0 = -58.3877
-      ox1 =  85.8079
-      ox2 =  23.8439
-      ox3 =  -0.034892
-      ox4 =   0.015568
-      ox5 =  -0.0019387
-
       !  FROM THE NEW BIOOPTIC MODEL PROPOSED JM ANDRE, WE READ HERE
       !  A PRECOMPUTED ARRAY CORRESPONDING TO THE ATTENUATION COEFFICIENT
@@ -554,5 +309,5 @@
 
 
-      CALL p4zche        ! initialize the chemical constants
+      CALL p4z_che        ! initialize the chemical constants
 
 
