Index: trunk/NEMO/TOP_SRC/PISCES/p4zbio.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zbio.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zbio.F90	(revision 935)
@@ -0,0 +1,155 @@
+MODULE p4zbio
+   !!======================================================================
+   !!                         ***  MODULE p4zbio  ***
+   !! TOP :   PISCES bio-model
+   !!======================================================================
+   !! History :   1.0  !  2004     (O. Aumont) Original code
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'                                       PISCES bio-model
+   !!----------------------------------------------------------------------
+   !!   p4z_bio        :   computes the interactions between the different
+   !!                      compartments of PISCES
+   !!----------------------------------------------------------------------
+   USE oce_trc         !
+   USE trp_trc         ! 
+   USE sms             ! 
+   USE p4zsink         ! 
+   USE p4zopt          ! 
+   USE p4zlim          ! 
+   USE p4zprod         !
+   USE p4zmort         !
+   USE p4zmicro        ! 
+   USE p4zmeso         ! 
+   USE p4zrem          ! 
+   USE prtctl_trc
+   
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   p4z_bio    ! called in p4zprg.F90
+
+   !! * Shared module variables
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !:
+      xnegtr            ! Array used to indicate negative tracer values
+
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Header:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE p4z_bio ( kt, jnt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_bio  ***
+      !!
+      !! ** Purpose :   Ecosystem model in the whole ocean: computes the
+      !!              different interactions between the different compartments
+      !!              of PISCES
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kt, jnt
+      INTEGER  ::  ji, jj, jk, jn
+#if defined key_kriest
+      REAL(wp) ::  zcoef1, zcoef2
+#endif
+      CHARACTER (len=25) :: charout
+
+      !!---------------------------------------------------------------------
+
+      !     ASSIGN THE SHEAR RATE THAT IS USED FOR AGGREGATION
+      !     OF PHYTOPLANKTON AND DETRITUS
+
+      xdiss(:,:,:) = 0.01
+
+!!gm the use of nmld should be better here?
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( fsdepw(ji,jj,jk+1) .le. hmld(ji,jj) )   xdiss(ji,jj,jk) = 1.e0
+            END DO 
+         END DO
+      END DO
+
+          
+      CALL p4z_sink ( kt, jnt )     ! vertical flux of particulate organic matter
+
+      CALL p4z_opt  ( kt, jnt )     ! Optic: PAR in the water column
+
+      CALL p4z_lim  ( kt, jnt )     ! co-limitations by the various nutrients
+
+      CALL p4z_prod ( kt, jnt )     ! phytoplankton growth rate over the global ocean. 
+      !                             ! (for each element : C, Si, Fe, Chl )
+
+      CALL p4z_rem  ( kt, jnt )     ! remineralization terms of organic matter+scavenging of Fe
+
+      CALL p4z_mort ( kt, jnt )     ! phytoplankton mortality
+
+      !                             ! zooplankton sources/sinks routines 
+      CALL p4z_micro( kt, jnt )           ! microzooplankton
+      CALL p4z_meso ( kt, jnt )           ! mesozooplankton
+
+
+      !                                                      ! test if tracers concentrations fall below 0.
+      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) + tra(ji,jj,jk,jn) ) < 0.e0 )    xnegtr(ji,jj,jk) = 0.e0
+               END DO
+            END DO
+         END DO
+      END DO
+      !                                ! where at least 1 tracer concentration becomes negative
+      !                                ! all tracer tendancy are set to zero (i.e. trn = trb)
+      DO jn = jp_pcs0, jp_pcs1
+         trn(:,:,:,jn) = trn(:,:,:,jn) + xnegtr(:,:,:) * tra(:,:,:,jn)
+      END DO
+
+
+      tra(:,:,:,:) = 0.e0
+
+#if defined key_kriest
+      ! 
+      zcoef1 = 1.e0 / xkr_massp 
+      zcoef2 = 1.e0 / xkr_massp / 1.1
+      DO jk = 1,jpkm1
+         trn(:,:,jk,jpnum) = MAX(  trn(:,:,jk,jpnum), trn(:,:,jk,jppoc) * zcoef1 / xnumm(jk)  )
+         trn(:,:,jk,jpnum) = MIN(  trn(:,:,jk,jpnum), trn(:,:,jk,jppoc) * zcoef2              )
+      END DO
+#endif
+
+
+# if defined key_trc_dia3d && defined key_kriest
+!!gm potential bug  hard coded index on trc3d
+      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
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No PISCES bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE p4z_bio                         ! Empty routine
+   END SUBROUTINE p4z_bio
+#endif 
+
+   !!======================================================================
+END MODULE  p4zbio
Index: trunk/NEMO/TOP_SRC/PISCES/p4zche.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zche.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zche.F90	(revision 935)
@@ -0,0 +1,375 @@
+MODULE p4zche
+   !!======================================================================
+   !!                         ***  MODULE p4zche  ***
+   !! TOP :   PISCES Sea water chemistry computed following OCMIP protocol
+   !!======================================================================
+   !! History :    -   !  1988     (E. Maier-Reimer)  Original code
+   !!              -   !  1998     (O. Aumont)  addition
+   !!              -   !  1999     (C. Le Quere)  modification
+   !!             1.0  !  2004     (O. Aumont)  modification
+   !!              -   !  2006     (R. Gangsto)  modification
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'                                       PISCES bio-model
+   !!----------------------------------------------------------------------
+   !!   p4z_che        :  Sea water chemistry computed following OCMIP protocol
+   !!----------------------------------------------------------------------
+   USE oce_trc         !
+   USE trp_trc         !
+   USE sms             ! 
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   p4z_che    ! called in p4zprg.F90
+
+   !! * Shared module variables
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !:
+      sio3eq, fekeq           !: chemistry of Fe and Si
+
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2)   ::   &   !:
+      chemc                   !: Solubilities of O2 and CO2
+
+   !! * 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
+#include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Header:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE p4z_che
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_che  ***
+      !!
+      !! ** Purpose :   Sea water chemistry computed following OCMIP protocol
+      !!
+      !! ** Method  : - ...
+      !!---------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk
+      REAL(wp) ::   ztkel, zsal , zqtt  , zbuf1 , zbuf2
+      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) ::   zst  , zft  , zcks  , zckf  , zaksp1
+      !!---------------------------------------------------------------------
+
+      ! CHEMICAL CONSTANTS - SURFACE LAYER
+      ! ----------------------------------
+!CDIR NOVERRCHK
+      DO jj = 1, jpj
+!CDIR NOVERRCHK
+         DO ji = 1, jpi
+
+            !                             ! SET ABSOLUTE TEMPERATURE
+            ztkel = tn(ji,jj,1) + 273.16
+            zqtt  = ztkel * 0.01
+            zqtt2 = zqtt * zqtt
+            zsal  = sn(ji,jj,1) + (1.- tmask(ji,jj,1) ) * 35.
+            zlqtt = LOG( zqtt )
+
+            !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980)
+            !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS
+            zcek1 = ca0 + ca1 / zqtt + ca2 * zlqtt + ca3 * zqtt2 + zsal*( ca4 + ca5 * zqtt + ca6 * zqtt2 )
+
+            !                             ! LN(K0) OF SOLUBILITY OF O2 and N2 (EQ. 4, WEISS, 1970)
+            zoxy  = ox0 + ox1 / zqtt + ox2 * zlqtt + zsal * ( ox3 + ox4 * zqtt + ox5 * zqtt2 )
+
+            !                             ! SET SOLUBILITIES OF O2 AND CO2
+            chemc(ji,jj,1) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000.
+            chemc(ji,jj,2) = EXP( zoxy  ) * oxyco
+
+         END DO
+      END DO
+
+      ! CHEMICAL CONSTANTS - DEEP OCEAN
+      ! -------------------------------
+!CDIR NOVERRCHK
+      DO jk = 1, jpk
+!CDIR NOVERRCHK
+         DO jj = 1, jpj
+!CDIR NOVERRCHK
+            DO ji = 1, jpi
+
+               ! SET PRESSION
+               zpres   = 1.025e-1 * fsdept(ji,jj,jk)
+
+               ! SET ABSOLUTE TEMPERATURE
+               ztkel   = tn(ji,jj,jk) + 273.16
+               zqtt    = ztkel * 0.01
+               zsal    = sn(ji,jj,jk) + ( 1.-tmask(ji,jj,jk) ) * 35.
+               zsqrt  = SQRT( zsal )
+               zsal15  = zsqrt * zsal
+               zlogt  = LOG( ztkel )
+               ztr    = 1. / ztkel
+               zis    = 19.924 * zsal / ( 1000.- 1.005 * zsal )
+               zis2   = zis * zis
+               zisqrt = SQRT( zis )
+               ztc     = tn(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20.
+
+               ! CHLORINITY (WOOSTER ET AL., 1969)
+               zcl     = zsal * salchl
+
+               ! TOTAL SULFATE CONCENTR. [MOLES/kg soln]
+               zst     = st1 * zcl * st2
+
+               ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln]
+               zft     = ft1 * zcl * ft2
+
+               ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990)
+               zcks    = EXP(  ks1 * ztr + ks0 + ks2 * zlogt                           &
+                  &                     + ( ks3 * ztr + ks4 + ks5 * zlogt ) * zisqrt   &
+                  &                     + ( ks6 * ztr + ks7 + ks8 * zlogt ) * zis      &
+                  &                     + ks9 * ztr * zis * zisqrt + ks10 * ztr *zis2 + LOG( ks11 + ks12 *zsal )  )
+
+               ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79)
+               zckf    = EXP(  kf1 * ztr + kf0 + kf2 * zisqrt + LOG( kf3 + kf4 * zsal )  )
+
+               ! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE
+               zckb    = ( cb0 + cb1 * zsqrt + cb2  * zsal + cb3 * zsal15 + cb4 * zsal * zsal ) * ztr   &
+                  &    + ( cb5 + cb6 * zsqrt + cb7  * zsal )                                            &
+                  &    + ( cb8 + cb9 * zsqrt + cb10 * zsal ) * zlogt + cb11 * zsqrt * ztkel             &
+                  &    + LOG(  ( 1.+ zst / zcks + zft / zckf ) / ( 1.+ zst / zcks )  )
+!!gm zsal**2 to be replaced by a *...
+               zck1    = c10 * ztr + c11 + c12 * zlogt + c13 * zsal + c14 * zsal**2
+               zck2    = c20 * ztr + c21 + c22 * zsal   + c23 * zsal**2
+
+               ! PKW (H2O) (DICKSON AND RILEY, 1979)
+               zckw    = cw0 * ztr + cw1 + cw2 * zlogt + ( cw3 * ztr + cw4 + cw5 * zlogt ) * zsqrt + cw6 * zsal
+
+
+               ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER
+               !       (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983)
+               zaksp0  = akcc1 + akcc2 * ztkel + akcc3 * ztr + akcc4 * LOG10( ztkel )   &
+                  &   + ( akcc5 + akcc6 * ztkel + akcc7 * ztr ) * zsqrt + akcc8 * zsal + akcc9 * zsal15
+
+               ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
+               zak1    = 10**(zck1)
+               zak2    = 10**(zck2)
+               zakb    = EXP( zckb  )
+               zakw    = EXP( zckw )
+               zaksp1  = 10**(zaksp0)
+
+               ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970)
+               !        (REFERENCE TO CULBERSON & PYTKOQICZ (1968) AS MADE
+               !        IN BROECKER ET AL. (1982) IS INCORRECT; HERE RGAS IS
+               !        TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres  IN
+               !        DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS
+               !        MULTIPLIED BY LN(10.) TO ALLOW USE OF EXP-FUNCTION
+               !        WITH BASIS E IN THE FORMULA FOR AKSPP (CF. EDMOND
+               !        & GIESKES (1970), P. 1285-1286 (THE SMALL
+               !        FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE
+               !        SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285))
+               zcpexp  = zpres /(rgas*ztkel)
+               zcpexp2 = zpres * zpres/(rgas*ztkel)
+
+               ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE
+               !        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(4) + devk2(4) * ztc + devk3(4) * ztc * ztc )
+               zbuf2  = 0.5 * ( devk4(4) + devk5(4) * ztc )
+               akw3(ji,jj,jk) = zakw * 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(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 )
+
+
+               ! TOTAL BORATE CONCENTR. [MOLES/L]
+               borat(ji,jj,jk) = bor1 * zcl * bor2
+
+               ! Iron and SIO3 saturation concentration from ...
+               sio3eq(ji,jj,jk) = EXP(  LOG( 10.) * ( 6.44 - 968. / ztkel )  ) * 1.e-6
+               fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ( 273.15 + ztc ) )
+
+            END DO
+         END DO
+      END DO
+      !
+   END SUBROUTINE p4z_che
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No PISCES bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE p4z_che( kt )                   ! Empty routine
+      INTEGER, INTENT( in ) ::   kt
+      WRITE(*,*) 'p4z_che: You should not have seen this print! error?', kt
+   END SUBROUTINE p4z_che
+#endif 
+
+   !!======================================================================
+END MODULE  p4zche
Index: trunk/NEMO/TOP_SRC/PISCES/p4zday.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zday.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zday.F90	(revision 935)
@@ -0,0 +1,96 @@
+MODULE p4zday
+   !!======================================================================
+   !!                         ***  MODULE p4zday  ***
+   !! TOP :   PISCES compute the day length depending on latitude and the day
+   !!======================================================================
+   !! History :    -   !  1993     (E. Maier-Reimer) Original code GBC 1993
+   !!              -   !  1999     (C. Le Quere)
+   !!             1.0  !  2004     (O. Aumont) Original code
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'                                       PISCES bio-model
+   !!----------------------------------------------------------------------
+   !!   p4z_day       :   compute the day length depending on latitude and the day
+   !!----------------------------------------------------------------------
+   USE oce_trc         !
+   USE trp_trc         ! 
+   USE sms             ! 
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   p4z_day    ! called in p4zprod.F90
+
+   !! * Shared module variables
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::  &  !:
+      strn              !:  Day length in hours
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Header:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE p4z_day
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_day  ***
+      !!
+      !! ** Purpose :   compute the day length depending on latitude and the day
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      INTEGER  ::   ji, jj
+      INTEGER  ::   iyy
+      REAL(wp) ::   zrum, zdelta, zcodel, zphi, zargu
+      !!---------------------------------------------------------------------
+
+      ! Get year
+      ! --------
+
+      iyy = ndastp / 10000
+
+      IF(lwp) write(numout,*)
+      IF(lwp) write(numout,*) 'p4zday : - Julian day ', nday_year
+      IF(lwp) write(numout,*) '~~~~~~'
+
+
+      IF( nleapy == 1 .AND. MOD( iyy, 4 ) == 0 ) THEN
+         zrum = FLOAT( nday_year - 80 ) / 366.
+      ELSE
+         zrum = FLOAT( nday_year - 80 ) / 365.
+      ENDIF
+
+      zdelta = SIN( zrum * rpi * 2. ) * sin( rpi * 23.5 / 180. )
+      zcodel = ASIN( zdelta )
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zphi = gphit(ji,jj) * rpi / 180.
+            zargu = TAN( zcodel ) * TAN( zphi )
+            strn(ji,jj) = 0.e0
+            zargu = MIN(  1., zargu )
+            zargu = MAX( -1., zargu )
+            strn(ji,jj) = 24.- 2.* ACOS( zargu ) * 180./ rpi / 15.
+            strn(ji,jj) = MAX( strn(ji,jj), 0.e0 )
+         END DO
+      END DO
+      !
+   END SUBROUTINE p4z_day
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No PISCES bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE p4z_day                    ! Empty routine
+   END SUBROUTINE p4z_day
+#endif 
+
+   !!======================================================================
+END MODULE  p4zday
Index: trunk/NEMO/TOP_SRC/PISCES/p4zflx.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zflx.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zflx.F90	(revision 935)
@@ -0,0 +1,250 @@
+MODULE p4zflx
+   !!======================================================================
+   !!                         ***  MODULE p4zflx  ***
+   !! TOP :   PISCES CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
+   !!======================================================================
+   !! History :    -   !  1988-07  (E. MAIER-REIMER) Original code
+   !!              -   !  1998     (O. Aumont) additions
+   !!              -   !  1999     (C. Le Quere) modifications
+   !!             1.0  !  2004     (O. Aumont) modifications
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'                                       PISCES bio-model
+   !!----------------------------------------------------------------------
+   !!   p4z_flx       :   CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
+   !!   p4z_flx_init  :   Read the namelist
+   !!----------------------------------------------------------------------
+   USE trc
+   USE oce_trc         !
+   USE trp_trc
+   USE sms
+   USE prtctl_trc
+   USE p4zche
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   p4z_flx    ! called in p4zprg.F90
+
+   REAL(wp) :: &  ! pre-industrial atmospheric [co2] (ppm) 	
+     atcox  = 0.20946 ,   &  !:
+     atcco2 = 278.           !:
+
+   REAL(wp) ::   &
+     xrhoa  = 1.22   ,   &   !: Air density kg/m3
+     xcd    = 1.5e-3         !: drag coefficient
+
+   REAL(wp) :: &
+     tco2flx = 0.            !: Total flux of carbon per year
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Header:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE p4z_flx ( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_flx  ***
+      !!
+      !! ** Purpose :   CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kt
+      INTEGER  ::   ji, jj, jrorr, nspyr
+      REAL(wp) ::   zttc, ztx, zty, ztau, zws, zkgwan
+      REAL(wp) ::   zfld, zflu, zfld16, zflu16, zfact
+      REAL(wp) ::   zph, zah2, zbot, zdic, zalk, zschmitto2, zalka, zschmittco2
+      REAL(wp), DIMENSION(jpi,jpj) ::   zkgco2, zkgo2, zh2co3, ztco2flx
+      CHARACTER (len=25) :: charout
+
+      !!---------------------------------------------------------------------
+
+
+      IF( kt == nittrc000  )   CALL p4z_flx_init      ! Initialization (first time-step only)
+
+
+      ! -----------------------------------------------------
+      !     ASSIGNATION TO EXPONENTS IN THE LISS AND MERLIVAT
+      !     FORMULATION OF THE GAS EXCHANGE RATE
+      ! -----------------------------------------------------
+
+      nspyr = INT( raass / rdt )
+
+      ! SURFACE CHEMISTRY (PCO2 AND [H+] IN
+      !     SURFACE LAYER); THE RESULT OF THIS CALCULATION
+      !     IS USED TO COMPUTE AIR-SEA FLUX OF CO2
+
+      DO jrorr = 1, 10
+
+!CDIR NOVERRCHK
+         DO jj = 1, jpj
+!CDIR NOVERRCHK
+            DO ji = 1, jpi
+
+               ! DUMMY VARIABLES FOR DIC, H+, AND BORATE
+               zbot  = borat(ji,jj,1)
+               zfact = rhop(ji,jj,1) / 1000. + rtrn
+               zdic  = trn(ji,jj,1,jpdic) / zfact
+               zph   = MAX( hi(ji,jj,1), 1.e-10 ) / zfact
+               zalka = trn(ji,jj,1,jptal) / zfact
+
+               ! CALCULATE [ALK]([CO3--], [HCO3-])
+               zalk  = zalka - (  akw3(ji,jj,1) / zph - zph + zbot / ( 1.+ zph / akb3(ji,jj,1) )  )
+
+               ! CALCULATE [H+] AND [H2CO3]
+               zah2   = SQRT(  (zdic-zalk)**2 + 4.* ( zalk * ak23(ji,jj,1)   &
+                  &                                        / ak13(ji,jj,1) ) * ( 2.* zdic - zalk )  )
+               zah2   = 0.5 * ak13(ji,jj,1) / zalk * ( ( zdic - zalk ) + zah2 )
+               zh2co3(ji,jj) = ( 2.* zdic - zalk ) / ( 2.+ ak13(ji,jj,1) / zah2 ) * zfact
+               hi(ji,jj,1)   = zah2 * zfact
+            END DO
+         END DO
+      END DO
+
+
+      ! --------------
+      ! COMPUTE FLUXES
+      ! --------------
+
+      ! FIRST COMPUTE GAS EXCHANGE COEFFICIENTS
+      ! -------------------------------------------
+
+!CDIR NOVERRCHK
+      DO jj = 1, jpj
+!CDIR NOVERRCHK
+         DO ji = 1, jpi
+
+            zttc = MIN( 35., tn(ji,jj,1) )
+
+      ! Compute the schmidt Number both O2 and CO2
+      ! ------------------------------------------
+
+            zschmittco2 = 2073.1 - 125.62 * zttc + 3.6276 * zttc**2 - 0.043126 * zttc**3
+            zschmitto2  = 1953.4 - 128.0 * zttc + 3.9918 * zttc**2 - 0.050091 * zttc**3
+
+            ! Estimation of wind speed as a function of wind stress
+            ztx  = utau(ji,jj) * umask(ji,jj,1)
+            zty  = vtau(ji,jj) * vmask(ji,jj,1)
+            ztau = SQRT( ztx * ztx + zty * zty )
+            zws  = SQRT ( ztau / ( xrhoa * xcd ) )
+
+      ! Compute the piston velocity for O2 and CO2
+      ! ------------------------------------------
+
+            zkgwan = (  0.3 * zws * zws    &
+               &   + 2.5 * ( 0.5246 + zttc * ( 0.016256+zttc*0.00049946 ) ) ) &
+# if defined key_off_degrad
+               &         * facvol(ji,jj,1)      &
+# endif
+               &   / (100. * 3600.)* ( 1.- freeze(ji,jj) ) * tmask(ji,jj,1)  
+
+            ! COMPUTE GAS EXCHANGE FOR CO2
+            zkgco2(ji,jj) = zkgwan * SQRT( 660./ zschmittco2 )
+            zkgo2(ji,jj)  = zkgwan * SQRT( 660./ zschmitto2 )
+
+         END DO
+      END DO
+
+      ztco2flx(:,:) = 0.
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+
+            ! Compute CO2 flux for the sea and air
+            zfld = atcco2 * tmask(ji,jj,1) * chemc(ji,jj,1) * zkgco2(ji,jj)
+            zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj)
+            tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) / fse3t(ji,jj,1)
+
+            ! compute flux of carbon
+            ztco2flx(ji,jj) = ( zfld - zflu ) * rfact &
+               &             * e1t(ji,jj) * e2t(ji,jj) * tmask(ji,jj,1) * 1000.
+
+            ! Compute O2 flux 
+            zfld16 = atcox * chemc(ji,jj,2) *tmask(ji,jj,1) * zkgo2(ji,jj)
+            zflu16 = trn(ji,jj,1,jpoxy) * tmask(ji,jj,1) * zkgo2(ji,jj)
+            tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + ( zfld16 - zflu16 ) / fse3t(ji,jj,1)
+
+# if defined key_trc_diaadd
+            ! Save diagnostics
+            trc2d(ji,jj,1) = ( zfld - zflu ) * 1000.
+            trc2d(ji,jj,2) = ( zfld16 - zflu16 ) * 1000.
+            trc2d(ji,jj,3) = zkgco2(ji,jj)
+            trc2d(ji,jj,4) = atcco2 - zh2co3(ji,jj) / ( chemc(ji,jj,1) + rtrn )
+# endif
+         END DO
+      END DO
+      !
+
+      ! Total Flux of Carbon
+      DO jj = 1, jpj
+        DO ji = 1, jpi
+           tco2flx = tco2flx + ztco2flx(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,*) '-------------------- : ',tco2flx * 12./1E15
+        WRITE(numout,*) '(GtC/an)'
+        tco2flx = 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
+
+   SUBROUTINE p4z_flx_init
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE p4z_flx_init  ***
+      !!
+      !! ** Purpose :   Initialization of atmospheric conditions
+      !!
+      !! ** Method  :   Read the natext namelist and check the parameters
+      !!      called at the first timestep (nittrc000)
+      !! ** input   :   Namelist natext
+      !!
+      !!----------------------------------------------------------------------
+
+      NAMELIST/natext/ atcco2
+
+      REWIND( numnat )                     ! read numnat
+      READ  ( numnat, natext )
+
+      IF(lwp) THEN                         ! control print
+         WRITE(numout,*) ' '
+         WRITE(numout,*) ' Namelist parameters for air-sea exchange, natext'
+         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '    Atmospheric pCO2      atcco2      =', atcco2
+      ENDIF
+
+   END SUBROUTINE p4z_flx_init
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No PISCES bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE p4z_flx( kt )                   ! Empty routine
+      INTEGER, INTENT( in ) ::   kt
+      WRITE(*,*) 'p4z_flx: You should not have seen this print! error?', kt
+   END SUBROUTINE p4z_flx
+#endif 
+
+   !!======================================================================
+END MODULE  p4zflx
Index: trunk/NEMO/TOP_SRC/PISCES/p4zint.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zint.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zint.F90	(revision 935)
@@ -0,0 +1,89 @@
+MODULE p4zint
+   !!======================================================================
+   !!                         ***  MODULE p4zint  ***
+   !! TOP :   PISCES interpolation and computation of various accessory fields
+   !!======================================================================
+   !! History :   1.0  !  2004-03 (O. Aumont) Original code
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'                                       PISCES bio-model
+   !!----------------------------------------------------------------------
+   !!   p4z_int        :  interpolation and computation of various accessory fields
+   !!----------------------------------------------------------------------
+   USE oce_trc         !
+   USE trp_trc
+   USE sms
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   p4z_int    ! called in p4zprg.F90
+
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !:
+      Tgfunc,            &  !:  Temp. dependancy of various biological rates
+      Tgfunc2               !:  Temp. dependancy of mesozooplankton rates 
+
+   !! * Module variables
+   REAL(wp) :: &
+      xksilim = 16.5E-6        ! Half-saturation constant for the computation of the Si half-saturation constant
+
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Header:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE p4z_int
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_int  ***
+      !!
+      !! ** Purpose :   interpolation and computation of various accessory fields
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      !!
+      INTEGER  ::   ji, jj
+      REAL(wp) ::   zdum
+      !!---------------------------------------------------------------------
+
+      ! Computation of phyto and zoo metabolic rate
+      ! -------------------------------------------
+
+      Tgfunc (:,:,:) = EXP( 0.063913 * tn(:,:,:) )
+      Tgfunc2(:,:,:) = EXP( 0.07608  * tn(:,:,:) )
+
+      ! Computation of the silicon dependant half saturation
+      ! constant for silica uptake
+      ! ---------------------------------------------------
+
+      DO ji = 1, jpi
+         DO jj = 1, jpj
+            zdum = trn(ji,jj,1,jpsil) * trn(ji,jj,1,jpsil)
+            xksimax(ji,jj) = MAX( xksimax(ji,jj), ( 1.+ 7.* zdum / ( xksilim * xksilim + zdum ) ) * 1e-6 )
+         END DO
+      END DO
+
+      IF( nday_year == 365 ) THEN
+         xksi    = xksimax
+         xksimax = 0.e0
+      ENDIF
+      !
+   END SUBROUTINE p4z_int
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No PISCES bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE p4z_int                   ! Empty routine
+      WRITE(*,*) 'p4z_int: You should not have seen this print! error?'
+   END SUBROUTINE p4z_int
+#endif 
+
+   !!======================================================================
+END MODULE  p4zint
Index: trunk/NEMO/TOP_SRC/PISCES/p4zlim.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zlim.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zlim.F90	(revision 935)
@@ -0,0 +1,224 @@
+MODULE p4zlim
+   !!======================================================================
+   !!                         ***  MODULE p4zlim  ***
+   !! TOP :   PISCES 
+   !!======================================================================
+   !! History :   1.0  !  2004     (O. Aumont) Original code
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'                                       PISCES bio-model
+   !!----------------------------------------------------------------------
+   !!   p4z_lim        :   Compute the nutrients limitation terms 
+   !!   p4z_lim_init   :   Read the namelist 
+   !!----------------------------------------------------------------------
+   USE trc
+   USE oce_trc         !
+   USE trp_trc         ! 
+   USE sms             ! 
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   p4z_lim    ! called in p4zprg.F90
+
+   !! * Shared module variables
+   REAL(wp), PUBLIC ::   conc0  = 2.e-6_wp         ,  &  !:
+     &                   conc1  = 10.e-6_wp        ,  &  !:
+     &                   conc2  = 2.e-11_wp        ,  &  !:
+     &                   conc2m = 8.E-11_wp        ,  &  !:
+     &                   conc3  = 1.e-10_wp        ,  &  !:
+     &                   conc3m = 4.e-10_wp        ,  &  !:
+     &                   concnnh4  = 1.e-7_wp      ,  &  !:
+     &                   concdnh4  = 5.e-7_wp      ,  &   !:
+     &                   xksi1     = 2.E-6_wp      ,  &   !:
+     &                   xksi2     = 3.33E-6_wp    ,  &   !:
+     &                   xkdoc     = 417.E-6_wp    ,  &   !:
+     &                   caco3r    = 0.3_wp               !:
+
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Header:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE p4z_lim( kt, jnt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_lim  ***
+      !!
+      !! ** Purpose :   Compute the co-limitations by the various nutrients
+      !!              for the various phytoplankton species
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
+      INTEGER  ::   ji, jj, jk
+      REAL(wp) ::   zlim1, zlim2, zlim3, zlim4, zno3, zferlim
+      REAL(wp) ::   zconctemp, zconctemp2, zconctempn, zconctempn2
+      REAL(wp) ::   ztemp, zdenom
+      !!---------------------------------------------------------------------
+
+
+      IF( ( kt * jnt ) == nittrc000  )   CALL p4z_lim_init      ! Initialization (first time-step only)
+
+
+!  Tuning of the iron concentration to a minimum
+!  level that is set to the detection limit
+!  -------------------------------------
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zno3=trn(ji,jj,jk,jpno3)
+               zferlim = MAX( 1.5e-11*(zno3/40E-6)**2, 3e-12 )
+               zferlim = MIN( zferlim, 1.5e-11 )
+               trn(ji,jj,jk,jpfer) = MAX( trn(ji,jj,jk,jpfer), zferlim )
+            END DO
+         END DO
+      END DO
+
+!  Computation of a variable Ks for iron on diatoms
+!  taking into account that increasing biomass is
+!  made of generally bigger cells
+!  ------------------------------------------------
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zconctemp   = MAX( 0.e0 , trn(ji,jj,jk,jpdia)-5e-7 )
+               zconctemp2  = trn(ji,jj,jk,jpdia) - zconctemp
+               zconctempn  = MAX( 0.e0 , trn(ji,jj,jk,jpphy)-1e-6 )
+               zconctempn2 = trn(ji,jj,jk,jpphy) - zconctempn
+               concdfe(ji,jj,jk) = ( zconctemp2 * conc3 + conc3m * zconctemp)   &
+                  &              / ( trn(ji,jj,jk,jpdia) + rtrn )
+               concdfe(ji,jj,jk) = MAX( conc3, concdfe(ji,jj,jk) )
+               concnfe(ji,jj,jk) = ( zconctempn2 * conc2 + conc2m * zconctempn)   &
+                  &              / ( trn(ji,jj,jk,jpphy) + rtrn )
+               concnfe(ji,jj,jk) = MAX( conc2, concnfe(ji,jj,jk) )
+            END DO
+         END DO
+      END DO
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+    
+!      Michaelis-Menten Limitation term for nutrients
+!      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
+
+               zlim1 = xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk)
+               zlim2 = trn(ji,jj,jk,jppo4) / ( trn(ji,jj,jk,jppo4) + concnnh4          ) 
+               zlim3 = trn(ji,jj,jk,jpfer) / ( trn(ji,jj,jk,jpfer) + concnfe(ji,jj,jk) )
+               xlimphy(ji,jj,jk) = MIN( zlim1, zlim2, zlim3 )
+               zlim1 = trn(ji,jj,jk,jpnh4) / ( concnnh4 + trn(ji,jj,jk,jpnh4) )
+               zlim3 = trn(ji,jj,jk,jpfer) / ( conc2    + trn(ji,jj,jk,jpfer) )
+               zlim4 = trn(ji,jj,jk,jpdoc) / ( xkdoc   + trn(ji,jj,jk,jpdoc) )
+               xlimbac(ji,jj,jk) = MIN( zlim1, zlim2, zlim3 ) * zlim4
+
+            END DO
+         END DO
+      END DO
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+!   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 * zdenom
+               xdiatnh4(ji,jj,jk) = trn(ji,jj,jk,jpnh4) * conc1    * zdenom 
+
+               zlim1 = xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk)
+               zlim2 = trn(ji,jj,jk,jppo4) / ( trn(ji,jj,jk,jppo4) + concdnh4          )
+               zlim3 = trn(ji,jj,jk,jpsil) / ( trn(ji,jj,jk,jpsil) + xksi   (ji,jj)    )
+               zlim4 = trn(ji,jj,jk,jpfer) / ( trn(ji,jj,jk,jpfer) + concdfe(ji,jj,jk) )
+               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
+      END DO
+      !
+   END SUBROUTINE p4z_lim
+
+   SUBROUTINE p4z_lim_init
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE p4z_lim_init  ***
+      !!
+      !! ** Purpose :   Initialization of nutrient limitation parameters
+      !!
+      !! ** Method  :   Read the natlim namelist and check the parameters
+      !!      called at the first timestep (nittrc000)
+      !!
+      !! ** input   :   Namelist natlim
+      !!
+      !!----------------------------------------------------------------------
+
+      NAMELIST/natlim/ conc0, conc1, conc2, conc2m, conc3, conc3m,   &
+         &             concnnh4, concdnh4, xksi1, xksi2, xkdoc, caco3r
+
+      REWIND( numnat )                     ! read numnat
+      READ  ( numnat, natlim )
+
+      IF(lwp) THEN                         ! control print
+         WRITE(numout,*) ' '
+         WRITE(numout,*) ' Namelist parameters for nutrient limitations, natlim'
+         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '    mean rainratio                            caco3r    =', caco3r
+         WRITE(numout,*) '    NO3, PO4 half saturation                  conc0      =', conc0
+         WRITE(numout,*) '    half saturation constant for Si uptake    xksi1     =', xksi1
+         WRITE(numout,*) '    half saturation constant for Si/C         xksi2     =', xksi2
+         WRITE(numout,*) '    2nd half-sat. of DOC remineralization     xkdoc    =', xkdoc
+         WRITE(numout,*) '    Phosphate half saturation for diatoms     conc1     =', conc1
+         WRITE(numout,*) '    Iron half saturation for phyto            conc2     =', conc2
+         WRITE(numout,*) '    Max iron half saturation for phyto        conc2m    =', conc2m
+         WRITE(numout,*) '    Iron half saturation for diatoms          conc3     =', conc3
+         WRITE(numout,*) '    Maxi iron half saturation for diatoms     conc3m    =', conc3m
+         WRITE(numout,*) '    NH4 half saturation for phyto             concnnh4  =', concnnh4
+         WRITE(numout,*) '    NH4 half saturation for diatoms           concdnh4  =', concdnh4
+      ENDIF
+
+   END SUBROUTINE p4z_lim_init
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No PISCES bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE p4z_lim                   ! Empty routine
+   END SUBROUTINE p4z_lim
+#endif 
+
+   !!======================================================================
+END MODULE  p4zlim
Index: trunk/NEMO/TOP_SRC/PISCES/p4zlys.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zlys.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zlys.F90	(revision 935)
@@ -0,0 +1,206 @@
+MODULE p4zlys
+   !!======================================================================
+   !!                         ***  MODULE p4zlys  ***
+   !! TOP :   PISCES 
+   !!======================================================================
+   !! History :    -   !  1988-07  (E. MAIER-REIMER) Original code
+   !!              -   !  1998     (O. Aumont) additions
+   !!              -   !  1999     (C. Le Quere) modifications
+   !!             1.0  !  2004     (O. Aumont) modifications
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'                                       PISCES bio-model
+   !!----------------------------------------------------------------------
+   !!   p4z_lys        :   Compute the CaCO3 dissolution 
+   !!   p4z_lys_init   :   Read the namelist parameters
+   !!----------------------------------------------------------------------
+   USE trc
+   USE oce_trc         !
+   USE trp_trc
+   USE sms
+   USE prtctl_trc
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   p4z_lys    ! called in p4zprg.F90
+
+   !! * Shared module variables
+   REAL(wp), PUBLIC ::   kdca    = 0.327e3_wp   ,  &  !:
+     &                   nca  = 1.0_wp                !:
+
+   !! * Module variables
+   REAL(wp) :: &
+      calcon = 1.03E-2        ! mean calcite concentration [Ca2+] in sea water [mole/kg solution]
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Header:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE p4z_lys( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_lys  ***
+      !!
+      !! ** Purpose :   CALCULATES DEGREE OF CACO3 SATURATION IN THE WATER
+      !!                COLUMN, DISSOLUTION/PRECIPITATION OF CACO3 AND LOSS
+      !!                OF CACO3 TO THE CACO3 SEDIMENT POOL.
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt ! ocean time step
+      INTEGER  ::   ji, jj, jk, jn
+      REAL(wp) ::   zbot, zalk, zdic, zph, zremco3, zah2
+      REAL(wp) ::   zdispot, zfact, zalka
+      REAL(wp) ::   zomegaca, zexcess, zexcess0
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zco3
+      CHARACTER (len=25) :: charout
+      !!---------------------------------------------------------------------
+
+      IF( kt == nittrc000  )   CALL p4z_lys_init      ! Initialization (first time-step only)
+
+
+
+      !     -------------------------------------------
+      !     COMPUTE [CO3--] and [H+] CONCENTRATIONS
+      !     -------------------------------------------
+      
+      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
+
+                  ! SET DUMMY VARIABLE FOR TOTAL BORATE
+                  zbot  = borat(ji,jj,jk)
+                  zfact = rhop (ji,jj,jk) / 1000. + rtrn
+
+                  ! SET DUMMY VARIABLE FOR [H+]
+                  zph   = hi(ji,jj,jk) * tmask(ji,jj,jk) / zfact + ( 1.-tmask(ji,jj,jk) ) * 1.e-9
+
+                  ! SET DUMMY VARIABLE FOR [SUM(CO2)]GIVEN 
+                  zdic  = trn(ji,jj,jk,jpdic) / zfact
+                  zalka = trn(ji,jj,jk,jptal) / zfact
+
+                  ! CALCULATE [ALK]([CO3--], [HCO3-])
+                  zalk  = zalka - (  akw3(ji,jj,jk) / zph - zph   &
+                     &             + zbot / (1.+ zph / akb3(ji,jj,jk) )  )
+
+                  ! CALCULATE [H+] and [CO3--]
+                  zah2 = SQRT( (zdic-zalk)*(zdic-zalk)+   &
+                     &     4.*(zalk*ak23(ji,jj,jk)/ak13(ji,jj,jk))   &
+                     &     *(2*zdic-zalk))
+
+                  zah2=0.5*ak13(ji,jj,jk)/zalk*((zdic-zalk)+zah2)
+                  zco3(ji,jj,jk) = zalk/(2.+zah2/ak23(ji,jj,jk))*zfact
+
+                  hi(ji,jj,jk)  = zah2*zfact
+
+               END DO
+            END DO
+         END DO
+         !
+      END DO 
+
+      !     ---------------------------------------------------------
+      !        CALCULATE DEGREE OF CACO3 SATURATION AND CORRESPONDING
+      !        DISSOLOUTION AND PRECIPITATION OF CACO3 (BE AWARE OF
+      !        MGCO3)
+      !     ---------------------------------------------------------
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+               ! DEVIATION OF [CO3--] FROM SATURATION VALUE
+               zomegaca = ( calcon * zco3(ji,jj,jk) ) / aksp(ji,jj,jk)
+
+               ! SET DEGREE OF UNDER-/SUPERSATURATION
+               zexcess0 = MAX( 0., ( 1.- zomegaca ) )
+               zexcess  = zexcess0**nca
+
+               ! AMOUNT CACO3 (12C) THAT RE-ENTERS SOLUTION
+               !       (ACCORDING TO THIS FORMULATION ALSO SOME PARTICULATE
+               !       CACO3 GETS DISSOLVED EVEN IN THE CASE OF OVERSATURATION)
+# if defined key_off_degrad
+              zdispot = kdca * zexcess * trn(ji,jj,jk,jpcal) * facvol(ji,jj,jk)
+# else
+              zdispot = kdca * zexcess * trn(ji,jj,jk,jpcal)
+# endif
+
+              !  CHANGE OF [CO3--] , [ALK], PARTICULATE [CACO3],
+              !       AND [SUM(CO2)] DUE TO CACO3 DISSOLUTION/PRECIPITATION
+              zremco3 = zdispot / rmoss
+              zco3(ji,jj,jk) = zco3(ji,jj,jk) + zremco3 * rfact
+              tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + 2.*zremco3
+              tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) -    zremco3
+              tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) +    zremco3
+
+            END DO
+         END DO
+      END DO
+
+# if defined key_trc_dia3d
+         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
+
+   SUBROUTINE p4z_lys_init
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE p4z_lys_init  ***
+      !!
+      !! ** Purpose :   Initialization of CaCO3 dissolution parameters
+      !!
+      !! ** Method  :   Read the natcal namelist and check the parameters
+      !!      called at the first timestep (nittrc000)
+      !!
+      !! ** input   :   Namelist natcal
+      !!
+      !!----------------------------------------------------------------------
+
+      NAMELIST/natcal/ kdca, nca
+
+      REWIND( numnat )                     ! read numnat
+      READ  ( numnat, natcal )
+
+      IF(lwp) THEN                         ! control print
+         WRITE(numout,*) ' '
+         WRITE(numout,*) ' Namelist parameters for CaCO3 dissolution, natcal'
+         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '    diss. rate constant calcite (per month)   kdca      =', kdca
+         WRITE(numout,*) '    order of reaction for calcite dissolution nca       =', nca
+      ENDIF
+
+   END SUBROUTINE p4z_lys_init
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No PISCES bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE p4z_lys( kt )                   ! Empty routine
+      INTEGER, INTENT( in ) ::   kt
+      WRITE(*,*) 'p4z_lys: You should not have seen this print! error?', kt
+   END SUBROUTINE p4z_lys
+#endif 
+
+   !!======================================================================
+END MODULE  p4zlys
Index: trunk/NEMO/TOP_SRC/PISCES/p4zmeso.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zmeso.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zmeso.F90	(revision 935)
@@ -0,0 +1,346 @@
+MODULE p4zmeso
+   !!======================================================================
+   !!                         ***  MODULE p4zmeso  ***
+   !! TOP :   PISCES Compute the sources/sinks for mesozooplankton
+   !!======================================================================
+   !! History :   1.0  !  2002     (O. Aumont) Original code
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'                                       PISCES bio-model
+   !!----------------------------------------------------------------------
+   !!   p4z_meso       :   Compute the sources/sinks for mesozooplankton
+   !!   p4z_meso_init  :   Initialization of the parameters for mesozooplankton
+   !!----------------------------------------------------------------------
+   USE trc
+   USE oce_trc         !
+   USE trp_trc         ! 
+   USE sms             ! 
+   USE prtctl_trc
+   USE p4zint
+   USE p4zsink
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   p4z_meso         ! called in p4zbio.F90
+
+   !! * Shared module variables
+   REAL(wp), PUBLIC ::   &
+      xprefc = 1.0_wp       ,  &  !: 
+      xprefp = 0.2_wp       ,  &  !:
+      xprefz = 1.0_wp       ,  &  !:
+      xprefpoc = 0.0_wp     ,  &  !:
+      resrat2 = 0.005_wp    ,  &  !:
+      mzrat2 = 0.03_wp      ,  &  !:
+      grazrat2 = 0.7_wp     ,  &  !:
+      xkgraz2 = 20E-6_wp    ,  &  !:
+      unass2  = 0.3_wp      ,  & !:
+      sigma2 = 0.6_wp       ,  & !:
+      epsher2 = 0.33_wp     ,  & !:   
+      grazflux = 5.E3_wp 
+
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Header:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE p4z_meso( kt,jnt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_meso  ***
+      !!
+      !! ** Purpose :   Compute the sources/sinks for mesozooplankton
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
+      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
+      !!---------------------------------------------------------------------
+
+
+      IF( ( kt * jnt ) == nittrc000  )   CALL p4z_meso_init      ! Initialization (first time-step only)
+
+      zstep = rfact2 / rjjss      ! Time step duration for biology
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+               zcompam = MAX( ( trn(ji,jj,jk,jpmes) - 1.e-9 ), 0.e0 )
+# if defined key_off_degrad
+               zfact   = zstep * tgfunc(ji,jj,jk) * zcompam * facvol(ji,jj,jk)
+# else
+               zfact   = zstep * tgfunc(ji,jj,jk) * zcompam
+# endif
+
+!     Respiration rates of both zooplankton
+!     -------------------------------------
+               zrespz2(ji,jj,jk)  = resrat2 * zfact * ( 1. + 3. * nitrfac(ji,jj,jk) )        &
+                  &     * trn(ji,jj,jk,jpmes) / ( xkmort + trn(ji,jj,jk,jpmes) )
+
+!     Zooplankton mortality. A square function has been selected with
+!     no real reason except that it seems to be more stable and may
+!     mimic predation.
+!     ---------------------------------------------------------------
+               ztortz2(ji,jj,jk) = mzrat2 * 1.e6 * zfact * trn(ji,jj,jk,jpmes)
+               !
+            END DO
+         END DO
+      END DO
+
+
+      DO jk = 1,jpkm1
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               zcompadi  = MAX( ( trn(ji,jj,jk,jpdia) - 1.e-8 ), 0.e0 )
+               zcompaz   = MAX( ( trn(ji,jj,jk,jpzoo) - 1.e-8 ), 0.e0 )
+               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
+!     ------------------------
+               zdenom = 1. / (  xkgraz2 + xprefc   * trn(ji,jj,jk,jpdia)   &
+                  &                     + xprefz   * trn(ji,jj,jk,jpzoo)   &
+                  &                     + xprefp   * trn(ji,jj,jk,jpphy)   &
+                  &                     + xprefpoc * trn(ji,jj,jk,jppoc)  )
+
+               zgraze2 = grazrat2 * zstep * Tgfunc2(ji,jj,jk) * zdenom    &
+# if defined key_off_degrad
+                  &     * facvol(ji,jj,jk)          &
+# endif
+                  &     * trn(ji,jj,jk,jpmes)
+
+               zgrazd(ji,jj,jk)   = zgraze2 * xprefc   * zcompadi
+               zgrazz(ji,jj,jk)   = zgraze2 * xprefz   * zcompaz
+               zgrazn(ji,jj,jk)   = zgraze2 * xprefp   * zcompaph
+               zgrazpoc(ji,jj,jk) = zgraze2 * xprefpoc * zcompapoc
+
+               zgraznf(ji,jj,jk)  = zgrazn(ji,jj,jk)   * trn(ji,jj,jk,jpnfe) &
+                  &                                     / (trn(ji,jj,jk,jpphy) + rtrn)
+               zgrazf(ji,jj,jk)   = zgrazd(ji,jj,jk)   * trn(ji,jj,jk,jpdfe) &
+                  &                                    / (trn(ji,jj,jk,jpdia) + rtrn)
+               zgrazpof(ji,jj,jk) = zgrazpoc(ji,jj,jk) * trn(ji,jj,jk,jpsfe) &
+                  &                                   / (trn(ji,jj,jk,jppoc) + rtrn)
+            END DO
+         END DO
+      END DO
+      
+      
+      DO jk = 1,jpkm1
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               
+!    Mesozooplankton flux feeding on GOC
+!    ----------------------------------
+# if ! defined key_kriest
+               zgrazffe(ji,jj,jk) = grazflux * zstep * wsbio4(ji,jj,jk)          &
+#  if defined key_off_degrad
+                  &     * facvol(ji,jj,jk)          &
+#  endif
+                  &     * tgfunc2(ji,jj,jk) * trn(ji,jj,jk,jpgoc) * trn(ji,jj,jk,jpmes)
+
+               zgrazfff(ji,jj,jk) = zgrazffe(ji,jj,jk)       &
+                  &     * trn(ji,jj,jk,jpbfe) / (trn(ji,jj,jk,jpgoc) + rtrn)
+# else
+! KRIEST3
+!!               zgrazffe(ji,jj,jk) = 0.5 * 1.3e-2 / 5.5e-7 * 0.3 * zstep * wsbio3(ji,jj,jk)     &
+!!                  &     * tgfunc(ji,jj,jk) * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpmes)    &
+#  if defined key_off_degrad
+!!                  &     * facvol(ji,jj,jk)          &
+#  endif
+!!                  &     /  (trn(ji,jj,jk,jppoc) * 1.e7 + 0.1)
+
+              zgrazffe(ji,jj,jk) = grazflux * zstep * wsbio3(ji,jj,jk)     &
+#  if defined key_off_degrad
+                  &     * facvol(ji,jj,jk)          &
+#  endif
+                  &     * tgfunc2(ji,jj,jk) * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpmes)
+
+               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_dmeso
+#else
+               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zgrapoc2
+#endif
+            END DO
+         END DO
+      END DO
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !
+               !   Update the arrays TRA which contain the biological sources and sinks
+               !   --------------------------------------------------------------------
+               zmortz2 = ztortz2(ji,jj,jk) + zrespz2(ji,jj,jk)
+               tra(ji,jj,jk,jpmes) = tra(ji,jj,jk,jpmes) - zmortz2  &
+                  &    + epsher2 * ( zgrazd(ji,jj,jk) + zgrazz(ji,jj,jk) + zgrazn(ji,jj,jk) &
+                  &    + zgrazpoc(ji,jj,jk) + zgrazffe(ji,jj,jk) )
+               tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) - zgrazd(ji,jj,jk)
+               tra(ji,jj,jk,jpzoo) = tra(ji,jj,jk,jpzoo) - zgrazz(ji,jj,jk)
+               tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) - zgrazn(ji,jj,jk)
+               tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) - zgrazn(ji,jj,jk) * trn(ji,jj,jk,jpnch)  &
+                  &    / ( trn(ji,jj,jk,jpphy) + rtrn )
+               tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) - zgrazd(ji,jj,jk) * trn(ji,jj,jk,jpdch) &
+                  &    / ( trn(ji,jj,jk,jpdia) + rtrn )
+               tra(ji,jj,jk,jpbsi) = tra(ji,jj,jk,jpbsi) - zgrazd(ji,jj,jk) * trn(ji,jj,jk,jpbsi) &
+                  &    / ( trn(ji,jj,jk,jpdia) + rtrn )
+               tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) +  zgrazd(ji,jj,jk) * trn(ji,jj,jk,jpbsi) &
+                  &    / ( trn(ji,jj,jk,jpdia) + rtrn )
+               tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) -  zgraznf(ji,jj,jk)
+               tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) -  zgrazf(ji,jj,jk)
+
+               zprcaca = xfracal(ji,jj,jk) * 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_dmeso &
+                  &    - 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
+
+   SUBROUTINE p4z_meso_init
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE p4z_meso_init  ***
+      !!
+      !! ** Purpose :   Initialization of mesozooplankton parameters
+      !!
+      !! ** Method  :   Read the natmes namelist and check the parameters
+      !!      called at the first timestep (nittrc000)
+      !!
+      !! ** input   :   Namelist natmes
+      !!
+      !!----------------------------------------------------------------------
+
+      NAMELIST/natmes/ grazrat2,resrat2,mzrat2,xprefc, xprefp, &
+         &             xprefz, xprefpoc, xkgraz2, epsher2, sigma2, unass2, grazflux
+
+      REWIND( numnat )                     ! read numnat
+      READ  ( numnat, natmes )
+
+
+      IF(lwp) THEN                         ! control print
+         WRITE(numout,*) ' ' 
+         WRITE(numout,*) ' Namelist parameters for mesozooplankton, natmes'
+         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '    zoo preference for phyto                  xprefc    =', xprefc
+         WRITE(numout,*) '    zoo preference for POC                    xprefp    =', xprefp
+         WRITE(numout,*) '    zoo preference for zoo                    xprefz    =', xprefz
+         WRITE(numout,*) '    zoo preference for poc                    xprefpoc  =', xprefpoc
+         WRITE(numout,*) '    exsudation rate of mesozooplankton        resrat2   =', resrat2
+         WRITE(numout,*) '    mesozooplankton mortality rate            mzrat2    =', mzrat2
+         WRITE(numout,*) '    maximal mesozoo grazing rate              grazrat2  =', grazrat2
+         WRITE(numout,*) '    mesozoo flux feeding rate                 grazflux  =', grazflux
+         WRITE(numout,*) '    non assimilated fraction of P by mesozoo  unass2    =', unass2
+         WRITE(numout,*) '    Efficicency of Mesozoo growth             epsher2   =', epsher2
+         WRITE(numout,*) '    Fraction of mesozoo excretion as DOM      sigma2    =', sigma2
+         WRITE(numout,*) '    half sturation constant for grazing 2     xkgraz2   =', xkgraz2
+      ENDIF
+
+   END SUBROUTINE p4z_meso_init
+
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No PISCES bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE p4z_meso                    ! Empty routine
+   END SUBROUTINE p4z_meso
+#endif 
+
+   !!======================================================================
+END MODULE  p4zmeso
Index: trunk/NEMO/TOP_SRC/PISCES/p4zmicro.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zmicro.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zmicro.F90	(revision 935)
@@ -0,0 +1,272 @@
+MODULE p4zmicro
+   !!======================================================================
+   !!                         ***  MODULE p4zmicro  ***
+   !! TOP :   PISCES Compute the sources/sinks for microzooplankton
+   !!======================================================================
+   !! History :   1.0  !  2004     (O. Aumont) Original code
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'                                       PISCES bio-model
+   !!----------------------------------------------------------------------
+   !!   p4z_micro       :   Compute the sources/sinks for microzooplankton
+   !!   p4z_micro_init  :   Initialize and read the appropriate namelist
+   !!----------------------------------------------------------------------
+   USE trc
+   USE oce_trc         !
+   USE trp_trc         ! 
+   USE sms             ! 
+   USE prtctl_trc
+   USE p4zint
+   USE p4zsink
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   p4z_micro    ! called in p4zbio.F90
+
+   !! * Shared module variables
+   REAL(wp), PUBLIC ::   &
+      xpref2c = 0.0_wp       ,  &  !:
+      xpref2p = 0.5_wp       ,  &  !:
+      xpref2d = 0.5_wp       ,  &  !:
+      resrat = 0.03_wp      ,  &  !:
+      mzrat = 0.0_wp        ,  &  !:
+      grazrat = 4.0_wp      ,  &  !:
+      xkgraz = 20E-6_wp     ,  &  !:
+      unass  = 0.3_wp       ,  &  !:
+      sigma1 = 0.6_wp       ,  &  !:
+      epsher = 0.33_wp
+
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Header:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE p4z_micro( kt,jnt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_micro  ***
+      !!
+      !! ** Purpose :   Compute the sources/sinks for microzooplankton
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
+      INTEGER  :: ji, jj, jk
+      REAL(wp) :: zcompadi, zcompadi2, zcompaz , zcompaph, zcompapoc
+      REAL(wp) :: zgraze  , zdenom  , zdenom2
+      REAL(wp) :: zfact   , zstep   , zinano , zidiat, zipoc
+      REAL(wp) :: 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
+
+      !!---------------------------------------------------------------------
+
+      IF( ( kt * jnt ) == nittrc000  )   CALL p4z_micro_init      ! Initialization (first time-step only)
+
+
+        zstep = rfact2 / rjjss      ! Time step duration for biology
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+               zcompaz = MAX( ( trn(ji,jj,jk,jpzoo) - 1.e-9 ), 0.e0 )
+# if defined key_off_degrad
+               zfact   = zstep * tgfunc(ji,jj,jk) * zcompaz *facvol(ji,jj,jk)
+# else
+               zfact   = zstep * tgfunc(ji,jj,jk) * zcompaz
+# endif
+
+!     Respiration rates of both zooplankton
+!     -------------------------------------
+
+               zrespz(ji,jj,jk) = resrat * zfact  * ( 1.+ 3.* nitrfac(ji,jj,jk) )     &
+                  &            * trn(ji,jj,jk,jpzoo) / ( xkmort + trn(ji,jj,jk,jpzoo) )
+
+!     Zooplankton mortality. A square function has been selected with
+!     no real reason except that it seems to be more stable and may
+!     mimic predation.
+!     ---------------------------------------------------------------
+               ztortz(ji,jj,jk) = mzrat * 1.e6 * zfact * trn(ji,jj,jk,jpzoo)
+
+            END DO
+         END DO
+      END DO
+
+
+ 
+      DO jk = 1,jpkm1
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               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
+               !     ------------------------
+               zdenom2 = 1./ ( xpref2p * zcompaph + xpref2c * zcompapoc + xpref2d * zcompadi2 + rtrn )
+
+               zgraze = grazrat * zstep * tgfunc(ji,jj,jk)     &
+# if defined key_off_degrad
+                  &      * facvol(ji,jj,jk)         &
+# endif
+                  &      * trn(ji,jj,jk,jpzoo)
+
+               zinano = xpref2p * zcompaph  * zdenom2
+               zipoc  = xpref2c * zcompapoc * zdenom2
+               zidiat = xpref2d * zcompadi2 * zdenom2
+
+               zdenom = 1./ ( xkgraz + zinano * zcompaph + zipoc * zcompapoc + zidiat * zcompadi2 )
+
+               zgrazp(ji,jj,jk)  = zgraze * zinano * zcompaph * zdenom
+               zgrazm(ji,jj,jk)  = zgraze * zipoc  * zcompapoc * zdenom
+               zgrazsd(ji,jj,jk) = zgraze * zidiat * zcompadi2 * zdenom
+
+               zgrazpf (ji,jj,jk) = zgrazp(ji,jj,jk)  * trn(ji,jj,jk,jpnfe) / (trn(ji,jj,jk,jpphy) + rtrn)
+               zgrazmf(ji,jj,jk)  = zgrazm(ji,jj,jk)  * trn(ji,jj,jk,jpsfe) / (trn(ji,jj,jk,jppoc) + rtrn)
+               zgrazsf(ji,jj,jk)  = zgrazsd(ji,jj,jk) * trn(ji,jj,jk,jpdfe) / (trn(ji,jj,jk,jpdia) + rtrn)
+
+            END DO
+         END DO
+      END DO
+      
+
+      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_ddiat
+#endif
+            END DO
+         END DO
+      END DO
+
+!
+!   Update the arrays TRA which contain the biological sources and sinks
+!   --------------------------------------------------------------------
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+               zmortz = ztortz(ji,jj,jk) + zrespz(ji,jj,jk)
+               tra(ji,jj,jk,jpzoo) = tra(ji,jj,jk,jpzoo) - zmortz  &
+                 &     + epsher * ( zgrazp(ji,jj,jk) + zgrazm(ji,jj,jk) + zgrazsd(ji,jj,jk))
+               tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) - zgrazp(ji,jj,jk)
+               tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) - zgrazsd(ji,jj,jk)
+               tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) - zgrazp(ji,jj,jk)  &
+                 &     * trn(ji,jj,jk,jpnch)/(trn(ji,jj,jk,jpphy)+rtrn)
+               tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) - zgrazsd(ji,jj,jk) &
+                 &     * trn(ji,jj,jk,jpdch)/(trn(ji,jj,jk,jpdia)+rtrn)
+               tra(ji,jj,jk,jpbsi) = tra(ji,jj,jk,jpbsi) - zgrazsd(ji,jj,jk) &
+                 &     * trn(ji,jj,jk,jpbsi)/(trn(ji,jj,jk,jpdia)+rtrn)
+               tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) + zgrazsd(ji,jj,jk) &
+                 &     * trn(ji,jj,jk,jpbsi)/(trn(ji,jj,jk,jpdia)+rtrn)
+               tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) - zgrazpf(ji,jj,jk)
+               tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) - zgrazsf(ji,jj,jk)
+               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zmortz - zgrazm(ji,jj,jk)
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ferat3 * zmortz   &
+                 &     + unass * ( zgrazpf(ji,jj,jk) + zgrazsf (ji,jj,jk)) &
+                 &     - (1.-unass) * zgrazmf(ji,jj,jk)
+               zprcaca = xfracal(ji,jj,jk) * 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_ddiat
+#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
+
+
+   SUBROUTINE p4z_micro_init
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE p4z_micro_init  ***
+      !!
+      !! ** Purpose :   Initialization of microzooplankton parameters
+      !!
+      !! ** Method  :   Read the natzoo namelist and check the parameters
+      !!      called at the first timestep (nittrc000)
+      !!
+      !! ** input   :   Namelist natzoo
+      !!
+      !!----------------------------------------------------------------------
+
+      NAMELIST/natzoo/ grazrat,resrat,mzrat,xpref2c, xpref2p, &
+         &             xpref2d, xkgraz, epsher, sigma1, unass
+
+      REWIND( numnat )                     ! read numnat
+      READ  ( numnat, natzoo )
+
+      IF(lwp) THEN                         ! control print
+         WRITE(numout,*) ' '
+         WRITE(numout,*) ' Namelist parameters for microzooplankton, natzoo'
+         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '    zoo preference for POC                    xpref2c    =', xpref2c
+         WRITE(numout,*) '    zoo preference for nano                   xpref2p    =', xpref2p
+         WRITE(numout,*) '    zoo preference for diatoms                xpref2d    =', xpref2d
+         WRITE(numout,*) '    exsudation rate of microzooplankton       resrat    =', resrat
+         WRITE(numout,*) '    microzooplankton mortality rate           mzrat     =', mzrat
+         WRITE(numout,*) '    maximal microzoo grazing rate             grazrat   =', grazrat
+         WRITE(numout,*) '    non assimilated fraction of P by microzoo unass     =', unass
+         WRITE(numout,*) '    Efficicency of microzoo growth            epsher    =', epsher
+         WRITE(numout,*) '    Fraction of microzoo excretion as DOM     sigma1    =', sigma1
+         WRITE(numout,*) '    half sturation constant for grazing 1     xkgraz    =', xkgraz
+      ENDIF
+
+   END SUBROUTINE p4z_micro_init
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No PISCES bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE p4z_micro                    ! Empty routine
+   END SUBROUTINE p4z_micro
+#endif 
+
+   !!======================================================================
+END MODULE  p4zmicro
Index: trunk/NEMO/TOP_SRC/PISCES/p4zmort.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zmort.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zmort.F90	(revision 935)
@@ -0,0 +1,277 @@
+MODULE p4zmort
+   !!======================================================================
+   !!                         ***  MODULE p4zmort  ***
+   !! TOP :   PISCES Compute the mortality terms for phytoplankton
+   !!======================================================================
+   !! History :   1.0  !  2002     (O. Aumont)  Original code
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'                                       PISCES bio-model
+   !!----------------------------------------------------------------------
+   !!   p4z_mort       :   Compute the mortality terms for phytoplankton
+   !!   p4z_mort_init  :   Initialize the mortality params for phytoplankton
+   !!----------------------------------------------------------------------
+   USE trc
+   USE oce_trc         !
+   USE trp_trc         ! 
+   USE sms             ! 
+   USE p4zsink
+   USE prtctl_trc
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   p4z_mort         ! called in p4zbio.F90
+
+
+   !! * Shared module variables
+   REAL(wp), PUBLIC ::   wchl = 0.001_wp      ,  &  !:
+     &                   wchld = 0.02_wp      ,  &  !:
+     &                   mprat = 0.01_wp      ,  &  !:
+     &                   mprat2 = 0.01_wp     ,  &  !:
+     &                   mpratm = 0.01_wp           !:
+
+   !! * Module variables
+   REAL(wp) :: zstep
+
+
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Header:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE p4z_mort( kt, jnt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_mort  ***
+      !!
+      !! ** Purpose :   Calls the different subroutine to initialize and compute
+      !!                the different phytoplankton mortality terms
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
+      !!---------------------------------------------------------------------
+
+      IF( ( kt * jnt ) == nittrc000  )   CALL p4z_mort_init      ! Initialization (first time-step only)
+
+      zstep = rfact2 / rjjss      ! Time step duration for biology
+
+      CALL p4z_nano            ! nanophytoplankton
+
+      CALL p4z_diat            ! diatoms
+
+   END SUBROUTINE p4z_mort
+
+
+   SUBROUTINE p4z_nano
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_nano  ***
+      !!
+      !! ** Purpose :   Compute the mortality terms for nanophytoplankton
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      INTEGER  :: ji, jj, jk
+      REAL(wp) :: zcompaph
+      REAL(wp) :: zfactfe,zfactch,zprcaca,zfracal
+      REAL(wp) :: ztortp,zrespp,zmortp
+      CHARACTER (len=25) :: charout
+      !!---------------------------------------------------------------------
+
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+               zcompaph = MAX( ( trn(ji,jj,jk,jpphy) - 1e-8 ), 0.e0 )
+
+!     Squared mortality of Phyto similar to a sedimentation term during
+!     blooms (Doney et al. 1996)
+!     -----------------------------------------------------------------
+               zrespp = wchl * 1.e6 * zstep * xdiss(ji,jj,jk)   &
+# if defined key_off_degrad
+                  &        * facvol(ji,jj,jk)     &
+# endif
+                  &        * zcompaph * trn(ji,jj,jk,jpphy)
+
+!     Phytoplankton mortality. This mortality loss is slightly
+!     increased when nutrients are limiting phytoplankton growth
+!     as observed for instance in case of iron limitation.
+!     ----------------------------------------------------------
+               ztortp = mprat * zstep * trn(ji,jj,jk,jpphy)          &
+# if defined key_off_degrad
+                  &          * facvol(ji,jj,jk)     &
+# endif
+                  &   / ( xkmort + trn(ji,jj,jk,jpphy) ) * zcompaph
+
+
+               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_dnano + zrespp * xkr_ddiat
+               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
+
+   SUBROUTINE p4z_diat
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_diat  ***
+      !!
+      !! ** Purpose :   Compute the mortality terms for diatoms
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      INTEGER  ::  ji, jj, jk
+      REAL(wp) ::  zfactfe,zfactsi,zfactch, zcompadi
+      REAL(wp) ::  zrespp2, ztortp2, zmortp2
+      CHARACTER (len=25) :: charout
+ 
+      !!---------------------------------------------------------------------
+
+
+!    Aggregation term for diatoms is increased in case of nutrient
+!    stress as observed in reality. The stressed cells become more
+!    sticky and coagulate to sink quickly out of the euphotic zone
+!     ------------------------------------------------------------
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+               zcompadi = MAX( ( trn(ji,jj,jk,jpdia) - 1e-8), 0. )
+
+!    Aggregation term for diatoms is increased in case of nutrient
+!    stress as observed in reality. The stressed cells become more
+!    sticky and coagulate to sink quickly out of the euphotic zone
+!     ------------------------------------------------------------
+
+               zrespp2  = 1.e6 * zstep * (  wchl + wchld * ( 1.- xlimdia(ji,jj,jk) )  )    &
+# if defined key_off_degrad
+                  &       * facvol(ji,jj,jk)       &
+# endif
+                  &       * xdiss(ji,jj,jk) * zcompadi * trn(ji,jj,jk,jpdia)
+                                                                               
+
+!     Phytoplankton mortality. 
+!     ------------------------
+               ztortp2  = mprat2 * zstep * trn(ji,jj,jk,jpdia)     &
+# if defined key_off_degrad
+                  &        * facvol(ji,jj,jk)       &
+# endif
+                  &      / ( xkmort + trn(ji,jj,jk,jpdia) ) * zcompadi
+
+                zmortp2 = zrespp2 + ztortp2
+
+!   Update the arrays tra which contains the biological sources and sinks
+!   ---------------------------------------------------------------------
+               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 )
+
+               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_ddiat + zrespp2 * xkr_daggr
+               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
+
+   SUBROUTINE p4z_mort_init
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE p4z_mort_init  ***
+      !!
+      !! ** Purpose :   Initialization of phytoplankton parameters
+      !!
+      !! ** Method  :   Read the natmort namelist and check the parameters
+      !!      called at the first timestep (nittrc000)
+      !!
+      !! ** input   :   Namelist natmort
+      !!
+      !!----------------------------------------------------------------------
+
+      NAMELIST/natmort/ wchl, wchld, mprat, mprat2, mpratm
+
+      REWIND( numnat )                     ! read numnat
+      READ  ( numnat, natmort )
+
+      IF(lwp) THEN                         ! control print
+         WRITE(numout,*) ' '
+         WRITE(numout,*) ' Namelist parameters for phytoplankton mortality, natmort'
+         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '    quadratic mortality of phytoplankton      wchl      =', wchl
+         WRITE(numout,*) '    maximum quadratic mortality of diatoms    wchld     =', wchld
+         WRITE(numout,*) '    phytoplankton mortality rate              mprat     =', mprat
+         WRITE(numout,*) '    Diatoms mortality rate                    mprat2    =', mprat2
+         WRITE(numout,*) '    Phytoplankton minimum mortality rate      mpratm    =', mpratm
+      ENDIF
+
+   END SUBROUTINE p4z_mort_init
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No PISCES bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE p4z_mort                    ! Empty routine
+   END SUBROUTINE p4z_mort
+#endif 
+
+   !!======================================================================
+END MODULE  p4zmort
Index: trunk/NEMO/TOP_SRC/PISCES/p4zopt.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zopt.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zopt.F90	(revision 935)
@@ -0,0 +1,310 @@
+MODULE p4zopt
+   !!======================================================================
+   !!                         ***  MODULE p4zopt  ***
+   !! TOP :   PISCES Compute the light availability in the water column
+   !!======================================================================
+   !! History :   1.0  !  2004     (O. Aumont) Original code
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'                                       PISCES bio-model
+   !!----------------------------------------------------------------------
+   !!   p4z_opt        :   Compute the light availability in the water column
+   !!----------------------------------------------------------------------
+   USE trc
+   USE oce_trc         !
+   USE trp_trc
+   USE sms
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   p4z_opt   ! called in p4zprg.F90
+
+   !! * Shared module variables
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !:
+      etot, enano, ediat,       &  !: PAR for phyto, nano and diat 
+      emoy                         !: averaged PAR in the mixed layer
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !:
+      heup                         !: Depth of the euphotic zone
+
+   !! * Module variables
+   REAL(wp), DIMENSION(3,61)                ::   &   !:
+      xkrgb                 !: ???
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Header:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE p4z_opt(kt, jnt)
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_opt  ***
+      !!
+      !! ** Purpose :   Compute the light availability in the water column
+      !!              depending on the depth and the chlorophyll concentration
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
+      INTEGER  ::   ji, jj, jk
+      INTEGER  ::   irgb
+      REAL(wp) ::   zchl, zparlux
+      REAL(wp) ::   zrlight , zblight , zglight
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zdepmoy, zetmp
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zekg, zekr, zekb
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze1 , ze2 , ze3
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3lum, ze4lum
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze5lum, ze6lum
+      !!---------------------------------------------------------------------
+
+
+      IF( ( kt * jnt ) == nittrc000  )   CALL p4z_opt_init      ! Initialization (first time-step only)
+
+
+!     Initialisation of variables used to compute PAR
+!     -----------------------------------------------
+      ze1 (:,:,:) = 0.e0
+      ze2 (:,:,:) = 0.e0
+      ze3 (:,:,:) = 0.e0
+      etot(:,:,:) = 0.e0
+        
+      zparlux = 0.43 / 3.
+
+!    IF activated, computation of the qsr for the dynamics
+!    -----------------------------------------------------
+      IF( ln_qsr_sms ) THEN
+         ze3lum(:,:,:) = 0.e0
+         ze4lum(:,:,:) = 0.e0
+         ze5lum(:,:,:) = 0.e0
+         ze6lum(:,:,:) = 0.e0
+      ENDIF
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+!     Separation in three light bands: red, green, blue
+!     -------------------------------------------------
+               zchl = ( trn(ji,jj,jk,jpnch) + trn(ji,jj,jk,jpdch) + rtrn ) * 1.e6
+               zchl = MAX( 0.03, zchl )
+               zchl = MIN( 10. , zchl )
+                                                                                
+               irgb = INT( 41 + 20.* LOG10( zchl ) + rtrn )
+                                                                                
+               zekb(ji,jj,jk) = xkrgb(1,irgb)
+               zekg(ji,jj,jk) = xkrgb(2,irgb)
+               zekr(ji,jj,jk) = xkrgb(3,irgb)
+
+            END DO
+         END DO
+      END DO
+
+!CDIR NOVERRCHK
+      DO jj = 1,jpj
+!CDIR NOVERRCHK
+         DO ji = 1,jpi
+
+!     Separation in three light bands: red, green, blue
+!     -------------------------------------------------
+
+            zblight = 0.5 * zekb(ji,jj,1) * fse3t(ji,jj,1)
+            zglight = 0.5 * zekg(ji,jj,1) * fse3t(ji,jj,1)
+            zrlight = 0.5 * zekr(ji,jj,1) * fse3t(ji,jj,1)
+
+            ze1(ji,jj,1) = zparlux * qsr(ji,jj) * EXP(-zblight)
+            ze2(ji,jj,1) = zparlux * qsr(ji,jj) * EXP(-zglight)
+            ze3(ji,jj,1) = zparlux * qsr(ji,jj) * EXP(-zrlight)
+
+         END DO
+      END DO
+
+!CDIR NOVERRCHK
+      DO jk = 2, jpkm1
+!CDIR NOVERRCHK
+          DO jj = 1, jpj
+!CDIR NOVERRCHK
+            DO ji = 1, jpi
+
+!     Separation in three light bands: red, green, blue
+!     -------------------------------------------------
+
+               zblight = 0.5 * ( zekb(ji,jj,jk-1) * fse3t(ji,jj,jk-1)   &
+                  &            + zekb(ji,jj,jk  ) * fse3t(ji,jj,jk  ) )
+               zglight = 0.5 * ( zekg(ji,jj,jk-1) * fse3t(ji,jj,jk-1)   &
+                  &            + zekg(ji,jj,jk  ) * fse3t(ji,jj,jk  ) )
+               zrlight = 0.5 * ( zekr(ji,jj,jk-1) * fse3t(ji,jj,jk-1)   &
+                  &            + zekr(ji,jj,jk  ) * fse3t(ji,jj,jk  ) )
+
+               ze1(ji,jj,jk) = ze1(ji,jj,jk-1) * EXP(-zblight)
+               ze2(ji,jj,jk) = ze2(ji,jj,jk-1) * EXP(-zglight)
+               ze3(ji,jj,jk) = ze3(ji,jj,jk-1) * EXP(-zrlight)
+
+            END DO
+         END DO
+      END DO
+
+      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
+
+!   In the following, the vertical attenuation of qsr for the dynamics is computed
+!   ------------------------------------------------------------------------------
+
+!CDIR NOVERRCHK
+         DO jj = 1, jpj
+!CDIR NOVERRCHK
+            DO ji = 1, jpi
+
+!     Separation in three light bands: red, green, blue
+!     -------------------------------------------------
+
+               zblight = 0.5 * zekb(ji,jj,1) * fse3t(ji,jj,1)
+               zglight = 0.5 * zekg(ji,jj,1) * fse3t(ji,jj,1)
+               zrlight = 0.5 * zekr(ji,jj,1) * fse3t(ji,jj,1)
+
+               ze3lum(ji,jj,1) = zparlux * qsr(ji,jj)
+               ze4lum(ji,jj,1) = zparlux * qsr(ji,jj)
+               ze5lum(ji,jj,1) = zparlux * qsr(ji,jj)
+               ze6lum(ji,jj,1) = (1.-3. * zparlux) * qsr(ji,jj)
+
+            END DO
+         END DO
+
+!CDIR NOVERRCHK
+         DO jk = 2, jpkm1
+!CDIR NOVERRCHK
+            DO jj = 1, jpj
+!CDIR NOVERRCHK
+               DO ji = 1, jpi
+
+!     Separation in three light bands: red, green, blue
+!     -------------------------------------------------
+
+                  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 )
+                  ze4lum(ji,jj,jk) = ze4lum(ji,jj,jk-1) * EXP( -zglight )
+                  ze5lum(ji,jj,jk) = ze5lum(ji,jj,jk-1) * EXP( -zrlight )
+                  ze6lum(ji,jj,jk) = ze6lum(ji,jj,jk-1) * EXP( -fse3t(ji,jj,jk-1) / xsi1 )
+
+               END DO
+            END DO
+         END DO
+
+         etot3(:,:,:) = ze3lum(:,:,:) + ze4lum(:,:,:) + ze5lum(:,:,:) + ze6lum(:,:,:)
+
+      ENDIF
+
+!     Computation of the euphotic depth
+!     ---------------------------------
+    
+      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) )   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
+!    --------------------------------------------------------
+
+      zdepmoy(:,:)   = 0.e0
+      zetmp  (:,:)   = 0.e0
+      emoy   (:,:,:) = 0.e0
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
+                  zetmp  (ji,jj) = zetmp  (ji,jj) + etot(ji,jj,jk) * fse3t(ji,jj,jk)
+                  zdepmoy(ji,jj) = zdepmoy(ji,jj) + fse3t(ji,jj,jk)
+               ENDIF
+            END DO
+         END DO
+      END DO
+
+      emoy(:,:,:) = etot(:,:,:)
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
+                  emoy(ji,jj,jk) = zetmp(ji,jj) / ( zdepmoy(ji,jj) + rtrn )
+               ENDIF
+            END DO
+         END DO
+      END DO
+
+
+# if defined key_trc_diaadd
+      trc2d(:,:,11) = heup(:,:)
+# endif
+      !
+   END SUBROUTINE p4z_opt
+
+   SUBROUTINE p4z_opt_init
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE p4z_opt_init  ***
+      !!
+      !! ** Purpose :   Initialization of of the optical scheme
+      !!
+      !! ** Method  :   read the look up table for the optical coefficients
+      !!
+      !! ** input   :   xKRGB61
+      !!
+      !!----------------------------------------------------------------------
+
+      INTEGER :: ichl, iband
+      INTEGER :: numlight
+      REAL(wp) ::   ztoto
+
+      !  FROM THE NEW BIOOPTIC MODEL PROPOSED JM ANDRE, WE READ HERE
+      !  A PRECOMPUTED ARRAY CORRESPONDING TO THE ATTENUATION COEFFICIENT
+
+      CALL ctlopn( numlight, 'kRGB61.txt', 'OLD', 'FORMATTED', 'SEQUENTIAL',   &
+         &           1, numout, .TRUE., 1 )
+
+      DO ichl = 1,61
+         READ(numlight,*) ztoto, ( xkrgb(iband,ichl), iband = 1,3 )
+      END DO
+
+      CLOSE(numlight)
+
+      IF(lwp) THEN                         ! control print
+         WRITE(numout,*) ' '
+         WRITE(numout,*) ' Initialization of the optical look-up table done'
+         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+      ENDIF
+
+   END SUBROUTINE p4z_opt_init
+
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No PISCES bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE p4z_opt                   ! Empty routine
+   END SUBROUTINE p4z_opt
+#endif 
+
+   !!======================================================================
+END MODULE  p4zopt
Index: trunk/NEMO/TOP_SRC/PISCES/p4zprod.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zprod.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zprod.F90	(revision 935)
@@ -0,0 +1,421 @@
+MODULE p4zprod
+   !!======================================================================
+   !!                         ***  MODULE p4zprod  ***
+   !! TOP :   PISCES 
+   !!======================================================================
+   !! History :   1.0  !  2004     (O. Aumont) Original code
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'                                       PISCES bio-model
+   !!----------------------------------------------------------------------
+   !!   p4z_prod       :  
+   !!----------------------------------------------------------------------
+   USE trc
+   USE oce_trc         !
+   USE trp_trc         ! 
+   USE sms             ! 
+   USE p4zday          !
+   USE prtctl_trc
+   USE p4zint
+   USE p4zlim
+   USE p4zopt
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   p4z_prod    ! called in p4zbio.F90
+
+   !! * Shared module variables
+   REAL(wp), PUBLIC ::   pislope   = 3.0_wp          ,  &  !:
+     &                   pislope2  = 3.0_wp          ,  &  !:
+     &                   excret  = 10.e-5_wp         , &   !:
+     &                   excret2 = 0.05_wp           , &   !:
+     &                   chlcnm  = 0.033_wp          , &   !:
+     &                   chlcdm  = 0.05_wp           , &   !:
+     &                   fecnm  = 10.E-6_wp          , &   !:
+     &                   fecdm  = 15.E-6_wp          , &   !:
+     &                   grosip = 0.151_wp
+
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk)  ::        &
+     &                   prmax
+   
+   REAL(wp) ::   &
+     tpp = 0.                               !: Total primary production
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Header:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE p4z_prod( kt , jnt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_prod  ***
+      !!
+      !! ** Purpose :   Compute the phytoplankton production depending on
+      !!              light, temperature and nutrient availability
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kt, jnt
+      INTEGER  ::   ji, jj, jk, nspyr
+      REAL(wp) ::   zsilfac, zfact, zrfact2
+      REAL(wp) ::   zprdiachl, zprbiochl, zsilim, ztn, zadap, zadap2
+      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, zstrn
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpislopead , zpislopead2
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zprdia     , zprbio, zysopt
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zprorca    , zprorcad, zprofed
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zprofen   , zprochln, zprochld
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpronew    , zpronewd
+      CHARACTER (len=25) :: charout
+      !!---------------------------------------------------------------------
+
+
+      IF( ( kt * jnt ) == nittrc000  )   CALL p4z_prod_init      ! Initialization (first time-step only)
+
+
+      zprorca (:,:,:) = 0.0
+      zprorcad(:,:,:) = 0.0
+      zprofed(:,:,:) = 0.0
+      zprofen(:,:,:) = 0.0
+      zprochln(:,:,:) = 0.0
+      zprochld(:,:,:) = 0.0
+      zpronew (:,:,:) = 0.0
+      zpronewd(:,:,:) = 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
+!     -------------------------------------
+
+# if defined key_off_degrad
+      prmax(:,:,:) = 0.6 / rjjss * tgfunc(:,:,:) * facvol(:,:,:)
+# else
+      prmax(:,:,:) = 0.6 / rjjss * tgfunc(:,:,:)
+# endif
+
+      CALL p4z_day       ! Computation of the day length
+
+!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
+!      --------------------------------------------------
+                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 * 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
+!    ------------------------
+!    Si/C increases with iron stress and silicate availability
+!    Si/C is arbitrariliy increased for very high Si concentrations
+!    to mimic the very high ratios observed in the Southern Ocean (silpot2)
+
+                  zlim1  = trn(ji,jj,jk,jpsil) / ( trn(ji,jj,jk,jpsil) + xksi1 )
+                  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)
+                  zysopt(ji,jj,jk) = grosip * zlim1 * zsilfac
+
+              ENDIF
+            END DO
+         END DO
+      END DO
+
+!    Computation of the limitation term due to
+!    A mixed layer deeper than the euphotic depth
+!    --------------------------------------------
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zmxltst = MAX( 0.e0, hmld(ji,jj) - heup(ji,jj) )
+            zmxlday = zmxltst**2 / rjjss
+            zmixnano(ji,jj) = 1.- zmxlday / ( 1.+ zmxlday )
+            zmixdiat(ji,jj) = 1.- zmxlday / ( 3.+ zmxlday )
+         END DO
+      END DO
+                                                                                
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
+
+!     Mixed-layer effect on production
+!     --------------------------------
+                  zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * zmixnano(ji,jj)
+                  zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * zmixdiat(ji,jj)
+               ENDIF
+            END DO
+         END DO
+      END DO
+
+!     Computation of the fractionnal day length
+!     -----------------------------------------
+      zstrn(:,:) = strn(:,:)
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+
+!      Computation of the maximum light intensity
+!      ------------------------------------------
+            IF( zstrn(ji,jj) < 1.e0 )   zstrn(ji,jj) = 24.
+            zstrn(ji,jj) = 24./zstrn(ji,jj)
+         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.
+!     ----------------------------------------------------------
+                  zetot2 = enano(ji,jj,jk) * zstrn(ji,jj)
+                  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 )
+                  zprod = rjjss * zprorca(ji,jj,jk) * zprbiochl * trn(ji,jj,jk,jpphy) *zmax
+
+                  zprofen(ji,jj,jk) = (fecnm)**2 * zprod / chlcnm            &
+                  &              / (  zpislopead(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpnfe) + rtrn  )
+
+                  zprochln(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
+!       -------------------------------------------------------
+                  zetot2 = ediat(ji,jj,jk) * zstrn(ji,jj)
+                  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 )  )
+
+                  zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * trn(ji,jj,jk,jpdia) * rfact2
+
+                  zpronewd(ji,jj,jk) = zprorcad(ji,jj,jk) * xdiatno3(ji,jj,jk)     &
+                  &              / ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) + rtrn )
+
+                  zprod = rjjss * zprorcad(ji,jj,jk) * zprdiachl * trn(ji,jj,jk,jpdia) * zmax
+
+                  zprofed(ji,jj,jk) = (fecdm)**2 * zprod / chlcdm                   &
+                  &              / ( zpislopead2(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpdfe) + rtrn )
+
+                  zprochld(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 = zprorcad(ji,jj,jk) - zpronewd(ji,jj,jk)
+              tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - zprorca(ji,jj,jk) - zprorcad(ji,jj,jk)
+              tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) - zpronew(ji,jj,jk) - zpronewd(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) + zprochln(ji,jj,jk) * zexcret
+              tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) + zprofen(ji,jj,jk) * zexcret
+              tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) + zprorcad(ji,jj,jk) * zexcret2
+              tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) + zprochld(ji,jj,jk) * zexcret2
+              tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) + zprofed(ji,jj,jk) * zexcret2
+              tra(ji,jj,jk,jpbsi) = tra(ji,jj,jk,jpbsi) + zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) * zexcret2
+              tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + &
+              &                     excret2 * zprorcad(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) + zpronewd(ji,jj,jk) )
+              tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) &
+              &                     - zexcret * zprofen(ji,jj,jk) - zexcret2 * zprofed(ji,jj,jk)
+              tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) &
+              &                     - zexcret2 * zprorcad(ji,jj,jk) * zysopt(ji,jj,jk)
+              tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprorca(ji,jj,jk) - zprorcad(ji,jj,jk)
+              tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) &
+              &                    + rno3 * ( zpronew(ji,jj,jk) + zpronewd(ji,jj,jk) )
+          END DO
+        END DO
+     END DO
+
+     ! Total primary production per year
+     DO jk = 1, jpkm1
+        DO jj = 1, jpj
+          DO ji = 1, jpi
+             tpp  = tpp + ( zprorca(ji,jj,jk) + zprorcad(ji,jj,jk) )  &
+#if defined key_off_degrad
+             &              * facvol(ji,jj,jk)   &
+#endif
+             &              * 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,*) '-------------------- : ',tpp * 12. / 1.E12
+        WRITE(numout,*) '(GtC/an)'
+        tpp = 0.
+      ENDIF
+
+#if defined key_trc_dia3d
+      zrfact2 = 1.e3 * rfact2r
+!   Supplementary diagnostics
+!   -------------------------
+      trc3d(:,:,:,4)  = etot(:,:,:)
+      trc3d(:,:,:,5)  = zprorca(:,:,:)  * zrfact2
+      trc3d(:,:,:,6)  = zprorcad(:,:,:) * zrfact2
+      trc3d(:,:,:,7)  = zpronew(:,:,:)  * zrfact2
+      trc3d(:,:,:,8)  = zpronewd(:,:,:) * zrfact2
+      trc3d(:,:,:,9)  = zprorcad(:,:,:) * zysopt(:,:,:) * zrfact2
+      trc3d(:,:,:,10) = zprofed(:,:,:) * zrfact2
+#if ! defined key_kriest
+      trc3d(:,:,:,11) = zprofen(:,:,:) * 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
+
+   SUBROUTINE p4z_prod_init
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE p4z_prod_init  ***
+      !!
+      !! ** Purpose :   Initialization of phytoplankton production parameters
+      !!
+      !! ** Method  :   Read the natprod namelist and check the parameters
+      !!      called at the first timestep (nittrc000)
+      !!
+      !! ** input   :   Namelist natprod
+      !!
+      !!----------------------------------------------------------------------
+
+      NAMELIST/natprod/ pislope, pislope2, excret, excret2, chlcnm, chlcdm,   &
+         &              fecnm, fecdm, grosip
+
+      REWIND( numnat )                     ! read numnat
+      READ  ( numnat, natprod )
+
+      IF(lwp) THEN                         ! control print
+         WRITE(numout,*) ' '
+         WRITE(numout,*) ' Namelist parameters for phytoplankton growth, natprod'
+         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '    mean Si/C ratio                           grosip    =', grosip
+         WRITE(numout,*) '    P-I slope                                 pislope   =', pislope
+         WRITE(numout,*) '    excretion ratio of nanophytoplankton      excret    =', excret
+         WRITE(numout,*) '    excretion ratio of diatoms                excret2   =', excret2
+         WRITE(numout,*) '    P-I slope  for diatoms                    pislope2  =', pislope2
+         WRITE(numout,*) '    Minimum Chl/C in nanophytoplankton        chlcnm    =', chlcnm
+         WRITE(numout,*) '    Minimum Chl/C in diatoms                  chlcdm    =', chlcdm
+         WRITE(numout,*) '    Maximum Fe/C in nanophytoplankton         fecnm     =', fecnm
+         WRITE(numout,*) '    Minimum Fe/C in diatoms                   fecdm     =', fecdm
+      ENDIF
+
+   END SUBROUTINE p4z_prod_init
+
+
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No PISCES bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE p4z_prod                    ! Empty routine
+   END SUBROUTINE p4z_prod
+#endif 
+
+   !!======================================================================
+END MODULE  p4zprod
Index: trunk/NEMO/TOP_SRC/PISCES/p4zrem.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zrem.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zrem.F90	(revision 935)
@@ -0,0 +1,449 @@
+MODULE p4zrem
+   !!======================================================================
+   !!                         ***  MODULE p4zrem  ***
+   !! TOP :   PISCES Compute remineralization/scavenging of organic compounds
+   !!======================================================================
+   !! History :   1.0  !  2004     (O. Aumont) Original code
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_top'       and                                      TOP models
+   !!   'key_pisces'                                       PISCES bio-model
+   !!----------------------------------------------------------------------
+   !!   p4z_rem       :   Compute remineralization/scavenging of organic compounds
+   !!----------------------------------------------------------------------
+   USE trc
+   USE oce_trc         !
+   USE trp_trc         ! 
+   USE sms             ! 
+   USE prtctl_trc
+   USE p4zint
+   USE p4zopt
+   USE p4zmeso
+   USE p4zprod
+   USE p4zche
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   p4z_rem    ! called in p4zbio.F90
+
+   !! * Shared module variables
+   REAL(wp), PUBLIC ::   xremik  = 0.3_wp      ,  & !:
+     &                   xremip  = 0.025_wp    ,  & !:
+     &                   nitrif  = 0.05_wp     ,  & !:
+     &                   xsirem  = 0.015_wp    ,  & !:
+     &                   xlam1   = 0.005_wp    ,  & !:
+     &                   oxymin  = 1.e-6_wp         !:
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::    & !:
+     &                   denitr                     !: denitrification array
+
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Header:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE p4z_rem(kt, jnt)
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_rem  ***
+      !!
+      !! ** Purpose :   Compute remineralization/scavenging of organic compounds
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
+      INTEGER  ::   ji, jj, jk
+      REAL(wp) ::   zremip, zremik , zlam1b
+      REAL(wp) ::   zkeq  , zfeequi, zsiremin
+      REAL(wp) ::   zsatur, zsatur2, znusil
+      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, zolimi
+      CHARACTER (len=25) :: charout
+
+      !!---------------------------------------------------------------------
+
+
+      IF( ( kt * jnt ) == nittrc000  )   CALL p4z_rem_init      ! Initialization (first time-step only)
+
+       zstep = rfact2 / rjjss      ! Time step duration for the biology
+
+
+!      Computation of the mean phytoplankton concentration as
+!      a crude estimate of the bacterial biomass
+!      --------------------------------------------------
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( fsdept(ji,jj,jk) < 120. ) THEN
+                  zdepbac(ji,jj,jk) = MIN( 0.7 * ( trn(ji,jj,jk,jpzoo) + 2.* trn(ji,jj,jk,jpmes) ), 4.e-6 )
+                  ztempbac(ji,jj)   = zdepbac(ji,jj,jk)
+               ELSE
+                  zdepbac(ji,jj,jk) = MIN( 1., 120./ fsdept(ji,jj,jk) ) * ztempbac(ji,jj)
+               ENDIF
+            END DO
+         END DO
+      END DO
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+!    DENITRIFICATION FACTOR COMPUTED FROM O2 LEVELS
+!    ----------------------------------------------
+
+               nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - trn(ji,jj,jk,jpoxy) )    &
+                  &                                / ( oxymin + trn(ji,jj,jk,jpoxy) )  )
+            END DO
+         END DO
+      END DO
+
+      nitrfac(:,:,:) = MIN( 1., nitrfac(:,:,:) )
+
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+!     DOC ammonification. Depends on depth, phytoplankton biomass
+!     and a limitation term which is supposed to be a parameterization
+!     of the bacterial activity. 
+!     ----------------------------------------------------------------
+               zremik = xremik * zstep / 1.e-6 * xlimbac(ji,jj,jk)         &
+# if defined key_off_degrad
+                  &            * facvol(ji,jj,jk)              &
+# endif
+                  &            * zdepbac(ji,jj,jk)
+               zremik = MAX( zremik, 5.5e-4 * zstep )
+
+!     Ammonification in oxic waters with oxygen consumption
+!     -----------------------------------------------------
+               zolimi(ji,jj,jk) = MIN(  ( trn(ji,jj,jk,jpoxy) - rtrn ) / o2ut,  &
+                  &                    zremik * ( 1.- nitrfac(ji,jj,jk) ) * trn(ji,jj,jk,jpdoc)  ) 
+
+!     Ammonification in suboxic waters with denitrification
+!     -------------------------------------------------------
+               denitr(ji,jj,jk) = MIN(  ( trn(ji,jj,jk,jpno3) - rtrn ) / rdenit,   &
+                  &                     zremik * nitrfac(ji,jj,jk) * trn(ji,jj,jk,jpdoc)  )
+            END DO
+         END DO
+      END DO
+
+      zolimi (:,:,:) = MAX( 0.e0, zolimi (:,:,:) )
+      denitr (:,:,:) = MAX( 0.e0, denitr (:,:,:) )
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+!    NH4 nitrification to NO3. Ceased for oxygen concentrations
+!    below 2 umol/L. Inhibited at strong light 
+!    ----------------------------------------------------------
+               zonitr  = nitrif * 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
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+!    Bacterial uptake of iron. No iron is available in DOC. So
+!    Bacteries are obliged to take up iron from the water. Some
+!    studies (especially at Papa) have shown this uptake to be
+!    significant
+!    ----------------------------------------------------------
+               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 )  )
+
+               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
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+!    POC disaggregation by turbulence and bacterial activity. 
+!    -------------------------------------------------------------
+               zremip = xremip * zstep * tgfunc(ji,jj,jk)   &
+# if defined key_off_degrad
+                  &            * facvol(ji,jj,jk)              &
+# endif
+                  &            * ( 1.- 0.5 * nitrfac(ji,jj,jk) )
+
+!    POC disaggregation rate is reduced in anoxic zone as shown by
+!    sediment traps data. In oxic area, the exponent of the martin s
+!    law is around -0.87. In anoxic zone, it is around -0.35. This
+!    means a disaggregation constant about 0.5 the value in oxic zones
+!    -----------------------------------------------------------------
+               zorem  = zremip * trn(ji,jj,jk,jppoc)
+               zofer  = zremip * trn(ji,jj,jk,jpsfe)
+#if ! defined key_kriest
+               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
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+!     Remineralization rate of BSi depedant on T and saturation
+!     ---------------------------------------------------------
+               zsatur  = ( sio3eq(ji,jj,jk) - trn(ji,jj,jk,jpsil) ) / ( sio3eq(ji,jj,jk) + rtrn )
+               zsatur  = MAX( rtrn, zsatur )
+               zsatur2 = zsatur * ( 1. + tn(ji,jj,jk) / 400.)**4
+               znusil  = 0.225  * ( 1. + tn(ji,jj,jk) / 15.) * zsatur + 0.775 * zsatur2**9
+#    if defined key_off_degrad
+               zsiremin = xsirem * zstep * znusil * facvol(ji,jj,jk)
+# else
+               zsiremin = xsirem * zstep * znusil
+#    endif
+               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
+         END DO
+      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
+!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
+!     load in particles on which they are adsorbed. The
+!     parameterization has been taken from studies on Th
+!     ------------------------------------------------------------
+               zkeq = fekeq(ji,jj,jk)
+               zfeequi = ( -( 1. + zfesatur(ji,jj,jk) * zkeq - zkeq * trn(ji,jj,jk,jpfer) )               &
+                  &        + SQRT( ( 1. + zfesatur(ji,jj,jk) * zkeq - zkeq * trn(ji,jj,jk,jpfer) )**2       &
+                  &               + 4. * trn(ji,jj,jk,jpfer) * zkeq) ) / ( 2. * zkeq )
+
+#if defined key_kriest
+               zlam1b = 3.e-5 + xlam1 * (  trn(ji,jj,jk,jppoc)                   &
+                  &                      + trn(ji,jj,jk,jpcal) + trn(ji,jj,jk,jpdsi)  ) * 1.e6
+#else
+               zlam1b = 3.e-5 + xlam1 * (  trn(ji,jj,jk,jppoc) + trn(ji,jj,jk,jpgoc)   &
+                  &                      + trn(ji,jj,jk,jpcal) + trn(ji,jj,jk,jpdsi)  ) * 1.e6
+#endif
+
+# if defined key_off_degrad
+               zscave = zfeequi * zlam1b * zstep  * facvol(ji,jj,jk)
+# else
+               zscave = zfeequi * zlam1b * zstep
+# endif
+
+!  Increased scavenging for very high iron concentrations
+!  found near the coasts due to increased lithogenic particles
+!  and let s say it unknown processes (precipitation, ...)
+!  -----------------------------------------------------------
+               zlamfac = MAX( 0.e0, ( gphit(ji,jj) + 55.) / 30. )
+               zlamfac = MIN( 1.  , zlamfac )
+#if ! defined key_kriest
+               zlam1b = (  80.* ( trn(ji,jj,jk,jpdoc) + 35.e-6 )                           &
+                  &     + 698.*   trn(ji,jj,jk,jppoc) + 1.05e4 * trn(ji,jj,jk,jpgoc)  )                    &
+                  &   * 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)  )                    &
+                  &   * 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
+               zaggdfe = zlam1b * zstep * 0.5 * ( trn(ji,jj,jk,jpfer) - zfeequi ) * facvol(ji,jj,jk)
+# else
+               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
+
+   SUBROUTINE p4z_rem_init
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE p4z_rem_init  ***
+      !!
+      !! ** Purpose :   Initialization of remineralization parameters
+      !!
+      !! ** Method  :   Read the natrem namelist and check the parameters
+      !!      called at the first timestep (nittrc000)
+      !!
+      !! ** input   :   Namelist natrem
+      !!
+      !!----------------------------------------------------------------------
+
+      NAMELIST/natrem/ xremik, xremip, nitrif, xsirem, xlam1, oxymin
+
+      REWIND( numnat )                     ! read numnat
+      READ  ( numnat, natrem )
+
+      IF(lwp) THEN                         ! control print
+         WRITE(numout,*) ' '
+         WRITE(numout,*) ' Namelist parameters for remineralization, natrem'
+         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '    remineralisation rate of POC              xremip    =', xremip
+         WRITE(numout,*) '    remineralization rate of DOC              xremik    =', xremik
+         WRITE(numout,*) '    remineralization rate of Si               xsirem    =', xsirem
+         WRITE(numout,*) '    scavenging rate of Iron                   xlam1     =', xlam1
+         WRITE(numout,*) '    NH4 nitrification rate                    nitrif    =', nitrif
+         WRITE(numout,*) '    halk saturation constant for anoxia       oxymin    =', oxymin
+      ENDIF
+
+   END SUBROUTINE p4z_rem_init
+
+
+
+
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No PISCES bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE p4z_rem                    ! Empty routine
+   END SUBROUTINE p4z_rem
+#endif 
+
+   !!======================================================================
+END MODULE p4zrem
Index: trunk/NEMO/TOP_SRC/PISCES/p4zsed.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zsed.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zsed.F90	(revision 935)
@@ -0,0 +1,616 @@
+MODULE p4zsed
+   !!======================================================================
+   !!                         ***  MODULE p4sed  ***
+   !! TOP :   PISCES Compute loss of organic matter in the sediments
+   !!======================================================================
+   !! History :   1.0  !  2004-03 (O. Aumont) Original code
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'                                       PISCES bio-model
+   !!----------------------------------------------------------------------
+   !!   p4z_sed        :  Compute loss of organic matter in the sediments
+   !!   p4z_sbc        :  Read and interpolate time-varying nutrients fluxes
+   !!   p4z_sed_init   :  Initialization of p4z_sed
+   !!----------------------------------------------------------------------
+   USE trc
+   USE oce_trc         !
+   USE trp_trc
+   USE sms
+   USE lib_mpp
+   USE prtctl_trc
+   USE p4zbio
+   USE p4zint
+   USE p4zopt
+   USE p4zsink
+   USE p4zrem
+   USE p4zlim
+   USE lbclnk
+   USE iom
+
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   p4z_sed    ! called in p4zprg.F90
+
+   !! * Shared module variables
+   LOGICAL, PUBLIC ::    bdustfer = .FALSE.      ,  &  !:
+     &                   briver = .FALSE.        ,  &  !:
+     &                   bndepo = .FALSE.        ,  &  !:
+     &                   bsedinput = .FALSE.           !:
+
+   REAL(wp), PUBLIC ::   sedfeinput = 1.E-9_wp   ,  &  !:
+     &                   dustsolub  = 0.014_wp         !:
+
+   !! * Module variables
+   INTEGER ::          &
+      numdust,          &  ! logical unit for surface fluxes data
+      nflx1 , nflx2,   &  !  first and second record used
+      nflx11, nflx12      ! ???
+   REAL(wp), DIMENSION(jpi,jpj,2) ::    &  !:
+     dustmo                                !: 2 consecutive set of dust fields 
+   REAL(wp), DIMENSION(jpi,jpj)   ::    &
+     rivinp, cotdep, nitdep, dust
+   REAL(wp), DIMENSION(jpi,jpj,jpk)  ::   &
+     ironsed
+   REAL(wp) :: sumdepsi, rivalkinput, rivpo4input, nitdepinput
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Header:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE p4z_sed(kt, jnt)
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_sed  ***
+      !!
+      !! ** Purpose :   Compute loss of organic matter in the sediments. This
+      !!              is by no way a sediment model. The loss is simply 
+      !!              computed to balance the inout from rivers and dust
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
+      INTEGER  ::   ji, jj, jk
+      INTEGER  ::   ikt
+      REAL(wp) ::   zsumsedsi, zsumsedpo4, zsumsedcal
+      REAL(wp) ::   zconctmp , zdenitot  , znitrpottot
+      REAL(wp) ::   zlim, zconctmp2, zstep, zfact
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zsidep
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   znitrpot, zirondep
+      CHARACTER (len=25) :: charout
+      !!---------------------------------------------------------------------
+
+
+      IF( ( kt * jnt ) == nittrc000  )   CALL p4z_sed_init      ! Initialization (first time-step only)
+
+      IF( (jnt == 1) .and. (bdustfer) )  CALL p4z_sbc( kt )
+
+      zstep = rfact2 / rjjss      ! Time step duration for the biology
+
+      zirondep(:,:,:) = 0.e0          ! Initialisation of variables used to compute deposition
+      zsidep  (:,:)   = 0.e0
+
+      ! Iron and Si deposition at the surface
+      ! -------------------------------------
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zirondep(ji,jj,1) = ( dustsolub * dust(ji,jj) / ( 55.85 * rmoss ) + 3.e-10 / raass )   &
+               &             * rfact2 / fse3t(ji,jj,1)
+            zsidep  (ji,jj)   = 8.8 * 0.075 * dust(ji,jj) * rfact2 / ( fse3t(ji,jj,1) * 28.1 * rmoss )
+         END DO
+      END DO
+
+      ! Iron solubilization of particles in the water column
+      ! ----------------------------------------------------
+
+      DO jk = 2, jpkm1
+         zirondep(:,:,jk) = dust(:,:) / ( 10. * 55.85 * rmoss ) * rfact2 * 1.e-4
+      END DO
+
+      ! Add the external input of nutrients, carbon and alkalinity
+      ! ----------------------------------------------------------
+
+      trn(:,:,1,jppo4) = trn(:,:,1,jppo4) + rivinp(:,:) * rfact2 
+      trn(:,:,1,jpno3) = trn(:,:,1,jpno3) + (rivinp(:,:) + nitdep(:,:)) * rfact2
+      trn(:,:,1,jpfer) = trn(:,:,1,jpfer) + rivinp(:,:) * 3.e-5 * rfact2
+      trn(:,:,1,jpsil) = trn(:,:,1,jpsil) + zsidep (:,:) + cotdep(:,:)   * rfact2 / 6.
+      trn(:,:,1,jpdic) = trn(:,:,1,jpdic) + rivinp(:,:) * 2.631 * rfact2
+      trn(:,:,1,jptal) = trn(:,:,1,jptal) + (cotdep(:,:) - rno3*(rivinp(:,:) +  nitdep(:,:) ) ) * rfact2
+
+
+      ! Add the external input of iron which is 3D distributed
+      ! (dust, river and sediment mobilization)
+      ! ------------------------------------------------------
+
+      DO jk = 1, jpkm1
+         trn(:,:,jk,jpfer) = trn(:,:,jk,jpfer)   &
+         &       + zirondep(:,:,jk) + ironsed(:,:,jk) * rfact2
+      END DO
+
+      ! Initialisation of variables used to compute Sinking Speed
+      ! ---------------------------------------------------------
+
+      zsumsedsi  = 0.e0
+      zsumsedpo4 = 0.e0
+      zsumsedcal = 0.e0
+
+      ! Loss of biogenic silicon, Caco3 organic carbon in the sediments. 
+      ! First, the total loss is computed.
+      ! The factor for calcite comes from the alkalinity effect
+      ! -------------------------------------------------------------
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ikt = MAX( mbathy(ji,jj)-1, 1 )
+            zfact = e1t(ji,jj) * e2t(ji,jj) / rjjss * tmask_i(ji,jj)
+
+# if defined key_kriest
+            zsumsedsi  = zsumsedsi  + zfact *  trn(ji,jj,ikt,jpdsi) * wscal (ji,jj,ikt)
+            zsumsedpo4 = zsumsedpo4 + zfact *  trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt)
+# else
+            zsumsedsi  = zsumsedsi  + zfact *  trn(ji,jj,ikt,jpdsi) * wsbio4(ji,jj,ikt)
+            zsumsedpo4 = zsumsedpo4 + zfact *( trn(ji,jj,ikt,jpgoc) * wsbio4(ji,jj,ikt)   &
+               &       + trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt) )
+# endif
+
+            zsumsedcal = zsumsedcal + zfact *  trn(ji,jj,ikt,jpcal) * wscal (ji,jj,ikt) * 2.e0
+
+         END DO
+      END DO
+
+      IF( lk_mpp ) THEN
+         CALL mpp_sum( zsumsedsi  )   ! sums over the global domain
+         CALL mpp_sum( zsumsedcal )   ! sums over the global domain
+         CALL mpp_sum( zsumsedpo4 )   ! sums over the global domain
+      ENDIF
+
+      ! Then this loss is scaled at each bottom grid cell for
+      ! equilibrating the total budget of silica in the ocean.
+      ! Thus, the amount of silica lost in the sediments equal
+      ! the supply at the surface (dust+rivers)
+      ! ------------------------------------------------------
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ikt = MAX( mbathy(ji,jj) - 1, 1 )
+            zconctmp = trn(ji,jj,ikt,jpdsi) * zstep / fse3t(ji,jj,ikt)   &
+# if ! defined key_kriest
+     &             * wscal (ji,jj,ikt)
+# else
+     &             * wsbio4(ji,jj,ikt)
+# endif
+
+            trn(ji,jj,ikt,jpdsi) = trn(ji,jj,ikt,jpdsi) - zconctmp
+            trn(ji,jj,ikt,jpsil) = trn(ji,jj,ikt,jpsil) + zconctmp   &
+            &      * ( 1.- ( sumdepsi + rivalkinput / raass / 6. ) / zsumsedsi )
+         END DO
+      END DO
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ikt = MAX( mbathy(ji,jj) - 1, 1 )
+            zconctmp = trn(ji,jj,ikt,jpcal) * wscal(ji,jj,ikt) * zstep / fse3t(ji,jj,ikt)
+            trn(ji,jj,ikt,jpcal) = trn(ji,jj,ikt,jpcal) - zconctmp
+            trn(ji,jj,ikt,jptal) = trn(ji,jj,ikt,jptal) + zconctmp   &
+               &   * ( 1.- ( rivalkinput / raass ) / zsumsedcal ) * 2.e0
+            trn(ji,jj,ikt,jpdic) = trn(ji,jj,ikt,jpdic) + zconctmp   &
+               &   * ( 1.- ( rivalkinput / raass ) / zsumsedcal )
+         END DO
+      END DO
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ikt = MAX( mbathy(ji,jj) - 1, 1 )
+            zfact = zstep / fse3t(ji,jj,ikt)
+
+# if ! defined key_kriest
+            zconctmp  = trn(ji,jj,ikt,jpgoc)
+            zconctmp2 = trn(ji,jj,ikt,jppoc)
+            trn(ji,jj,ikt,jpgoc) = trn(ji,jj,ikt,jpgoc) - zconctmp  * wsbio4(ji,jj,ikt) * zfact
+            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) - zconctmp2 * wsbio3(ji,jj,ikt) * zfact
+            trn(ji,jj,ikt,jpdoc) = trn(ji,jj,ikt,jpdoc)    &
+            &      + ( zconctmp  * wsbio4(ji,jj,ikt) + zconctmp2 * wsbio3(ji,jj,ikt) ) * zfact   &
+            &      * ( 1.- rivpo4input / (raass * zsumsedpo4 ) )
+            trn(ji,jj,ikt,jpbfe) = trn(ji,jj,ikt,jpbfe) - trn(ji,jj,ikt,jpbfe) * wsbio4(ji,jj,ikt) * zfact
+            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) - trn(ji,jj,ikt,jpsfe) * wsbio3(ji,jj,ikt) * zfact
+
+# else
+
+            zconctmp  = trn(ji,jj,ikt,jpnum)
+            zconctmp2 = trn(ji,jj,ikt,jppoc)
+            trn(ji,jj,ikt,jpnum) = trn(ji,jj,ikt,jpnum)   &
+            &      - zconctmp  * wsbio4(ji,jj,ikt) * zfact
+            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc)   &
+            &      - zconctmp2 * wsbio3(ji,jj,ikt) * zfact
+            trn(ji,jj,ikt,jpdoc) = trn(ji,jj,ikt,jpdoc)    &
+            &      + ( zconctmp2 * wsbio3(ji,jj,ikt) )   &
+            &      * zfact * ( 1.- rivpo4input / ( raass * zsumsedpo4 ) )
+            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe)   &
+            &      - trn(ji,jj,ikt,jpsfe) * wsbio3(ji,jj,ikt) * zfact
+
+# endif
+         END DO
+      END DO
+
+      ! Nitrogen fixation (simple parameterization). The total gain
+      ! from nitrogen fixation is scaled to balance the loss by 
+      ! denitrification
+      ! -------------------------------------------------------------
+
+      zdenitot = 0.e0
+      DO jk = 1, jpkm1
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               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) * xnegtr(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      IF( lk_mpp )   CALL mpp_sum( zdenitot )      ! sum over the global domain
+
+      ! Potential nitrogen fication dependant on temperature and iron
+      ! -------------------------------------------------------------
+
+!CDIR NOVERRCHK
+      DO jk = 1, jpk
+!CDIR NOVERRCHK
+         DO jj = 1, jpj
+!CDIR NOVERRCHK
+            DO ji = 1, jpi
+               zlim = ( 1.- xnanono3(ji,jj,jk) - xnanonh4(ji,jj,jk) )
+               IF( zlim <= 0.2 )   zlim = 0.01
+               znitrpot(ji,jj,jk) = MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) / rjjss )   &
+# if defined key_off_degrad
+               &                  * facvol(ji,jj,jk)   &
+# endif
+               &                  * zlim * rfact2 * trn(ji,jj,jk,jpfer)   &
+               &                  / ( conc3 + trn(ji,jj,jk,jpfer) ) * ( 1.- EXP( -etot(ji,jj,jk) / 50.) )
+            END DO
+         END DO 
+      END DO
+
+      znitrpottot = 0.e0
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               znitrpottot = znitrpottot + znitrpot(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)   &
+                  &                                           * tmask(ji,jj,jk) * tmask_i(ji,jj) 
+            END DO
+         END DO
+      END DO
+
+      IF( lk_mpp )   CALL mpp_sum( znitrpottot )  ! sum over the global domain
+
+      ! Nitrogen change due to nitrogen fixation
+      ! ----------------------------------------
+
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+# if ! defined key_c1d && ( defined key_orca_r4 || defined key_orca_r2 || defined key_orca_r05 || defined key_orca_r025 )
+!!             zfact = znitrpot(ji,jj,jk) * zdenitot / znitrpottot
+               zfact = znitrpot(ji,jj,jk) * 1.e-7
+# else
+               zfact = znitrpot(ji,jj,jk) * 1.e-7
+# endif
+               trn(ji,jj,jk,jpnh4) = trn(ji,jj,jk,jpnh4) + zfact
+               trn(ji,jj,jk,jpoxy) = trn(ji,jj,jk,jpoxy) + zfact   * o2nit
+               trn(ji,jj,jk,jppo4) = trn(ji,jj,jk,jppo4) + 30./ 46.* zfact
+            END DO
+         END DO
+      END DO
+
+# if defined key_trc_diaadd
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            trc2d(ji,jj,13) = znitrpot(ji,jj,1) * 1.e-7 * fse3t(ji,jj,1) * 1.e+3 / rfact2
+            trc2d(ji,jj,12) = zirondep(ji,jj,1) * 1.e+3 * rfact2r * fse3t(ji,jj,1)
+         END DO
+      END DO
+# 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
+
+   SUBROUTINE p4z_sbc(kt)
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE p4z_sbc  ***
+      !!
+      !! ** Purpose :   Read and interpolate the external sources of 
+      !!                nutrients
+      !!
+      !! ** Method  :   Read the files and interpolate the appropriate variables
+      !!
+      !! ** input   :   external netcdf files
+      !!
+      !!----------------------------------------------------------------------
+      !! * arguments
+      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
+
+      !! * Local declarations
+      INTEGER ::   &
+         imois, imois2,       &  ! temporary integers
+         i15  , iman             !    "          "
+      REAL(wp) ::   &
+         zxy                     !    "         "
+
+
+      !!---------------------------------------------------------------------
+
+      ! Initialization
+      ! --------------
+
+      i15 = INT( 2 * FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) )
+      iman  = INT( raamo )
+      imois = nmonth + i15 - 1
+      IF( imois == 0 ) imois = iman
+      imois2 = nmonth
+
+      ! 1. first call kt=nit000
+      ! -----------------------
+
+      IF( kt == nit000 ) THEN
+         ! initializations
+         nflx1  = 0
+         nflx11 = 0
+         ! open the file
+         IF(lwp) THEN
+            WRITE(numout,*) ' '
+            WRITE(numout,*) ' **** Routine p4z_sbc'
+         ENDIF
+         CALL iom_open ( 'dust.orca.nc', numdust )
+      ENDIF
+
+
+     ! Read monthly file
+      ! ----------------
+
+      IF( kt == nit000 .OR. imois /= nflx1 ) THEN
+
+         ! Calendar computation
+
+         ! nflx1 number of the first file record used in the simulation
+         ! nflx2 number of the last  file record
+
+         nflx1 = imois
+         nflx2 = nflx1+1
+         nflx1 = MOD( nflx1, iman )
+         nflx2 = MOD( nflx2, iman )
+         IF( nflx1 == 0 )   nflx1 = iman
+         IF( nflx2 == 0 )   nflx2 = iman
+         IF(lwp) WRITE(numout,*) 'first record file used nflx1 ',nflx1
+         IF(lwp) WRITE(numout,*) 'last  record file used nflx2 ',nflx2
+
+         ! Read monthly fluxes data
+
+         ! humidity
+         CALL iom_get ( numdust, jpdom_data, 'dust', dustmo(:,:,1), nflx1 )
+         CALL iom_get ( numdust, jpdom_data, 'dust', dustmo(:,:,2), nflx2 )
+
+         IF(lwp .AND. nitend-nit000 <= 100 ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) ' read clio flx ok'
+            WRITE(numout,*)
+               WRITE(numout,*)
+               WRITE(numout,*) 'Clio month: ',nflx1,'  field: dust'
+               CALL prihre( dustmo(:,:,1),jpi,jpj,1,jpi,20,1,jpj,10,1e9,numout )
+         ENDIF
+
+      ENDIF
+
+     ! 3. at every time step interpolation of fluxes
+      ! ---------------------------------------------
+
+      zxy = FLOAT( nday ) / FLOAT( nobis(nflx1) ) + 0.5 - i15
+
+      dust(:,:) = ( (1.-zxy) * dustmo(:,:,1) + zxy * dustmo(:,:,2) )
+
+      IF( kt == nitend ) CALL iom_close (numdust)
+
+   END SUBROUTINE p4z_sbc
+
+
+   SUBROUTINE p4z_sed_init
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE p4z_sed_init  ***
+      !!
+      !! ** Purpose :   Initialization of the external sources of nutrients
+      !!
+      !! ** Method  :   Read the files and compute the budget
+      !!      called at the first timestep (nittrc000)
+      !!
+      !! ** input   :   external netcdf files
+      !!
+      !!----------------------------------------------------------------------
+
+      INTEGER ::   ji, jj, jk, jm
+      INTEGER , PARAMETER ::   jpmois = 12, jpan = 1
+      INTEGER :: numriv, numbath, numdep
+
+
+      REAL(wp) ::   zcoef
+      REAL(wp) ::   expide, denitide,zmaskt
+      REAL(wp) , DIMENSION (jpi,jpj)     ::   riverdoc, river, ndepo
+      REAL(wp) , DIMENSION (jpi,jpj,jpk) ::   cmask
+      REAL(wp), DIMENSION(jpi,jpj,12)    ::   zdustmo
+
+      NAMELIST/natsms/ bdustfer, briver, bndepo, bsedinput, sedfeinput, dustsolub
+
+
+      REWIND( numnat )                     ! read numnat
+      READ  ( numnat, natsms )
+
+      IF(lwp) THEN
+         WRITE(numout,*) ' '
+         WRITE(numout,*) ' Namelist : natsms '
+         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '
+         WRITE(numout,*) '    Dust input from the atmosphere           bdustfer  = ', bdustfer
+         WRITE(numout,*) '    River input of nutrients                 briver    = ', briver
+         WRITE(numout,*) '    Atmospheric deposition of N              bndepo    = ', bndepo
+         WRITE(numout,*) '    Fe input from sediments                  bsedinput = ', bsedinput
+         WRITE(numout,*) '    Coastal release of Iron                  sedfeinput =', sedfeinput
+         WRITE(numout,*) '    Solubility of the dust                   dustsolub  =', dustsolub
+      ENDIF
+
+      ! Dust input from the atmosphere
+      ! ------------------------------
+      IF( bdustfer ) THEN 
+         IF(lwp) WRITE(numout,*) '    Initialize dust input from atmosphere '
+         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
+         CALL iom_open ( 'dust.orca.nc', numdust )
+         DO jm = 1, jpmois
+            CALL iom_get( numdust, jpdom_data, 'dust', zdustmo(:,:,jm), jm )
+         END DO
+         CALL iom_close( numdust )
+      ELSE
+         zdustmo(:,:,:) = 0.e0
+         dust(:,:) = 0.0
+      ENDIF
+
+      ! Nutrient input from rivers
+      ! --------------------------
+      IF( briver ) THEN
+         IF(lwp) WRITE(numout,*) '    Initialize the nutrient input by rivers from river.orca.nc file'
+         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+         CALL iom_open ( 'river.orca.nc', numriv )
+         CALL iom_get  ( numriv, jpdom_data, 'riverdic', river   (:,:), jpan )
+         CALL iom_get  ( numriv, jpdom_data, 'riverdoc', riverdoc(:,:), jpan )
+         CALL iom_close( numriv )
+      ELSE
+         river   (:,:) = 0.e0
+         riverdoc(:,:) = 0.e0
+      endif
+
+      ! Nutrient input from dust
+      ! ------------------------
+      IF( bndepo ) THEN
+         IF(lwp) WRITE(numout,*) '    Initialize the nutrient input by dust from ndeposition.orca.nc'
+         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+         CALL iom_open ( 'ndeposition.orca.nc', numdep )
+         CALL iom_get  ( numdep, jpdom_data, 'ndep', ndepo(:,:), jpan )
+         CALL iom_close( numdep )
+      ELSE
+         ndepo(:,:) = 0.e0
+      ENDIF
+
+      ! Coastal and island masks
+      ! ------------------------
+      IF( bsedinput ) THEN     
+         IF(lwp) WRITE(numout,*) '    Computation of an island mask to enhance coastal supply of iron'
+         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+         IF(lwp) WRITE(numout,*) '       from bathy.orca.nc file '
+         CALL iom_open ( 'bathy.orca.nc', numbath )
+         CALL iom_get  ( numbath, jpdom_data, 'bathy', cmask(:,:,:), jpan )
+         CALL iom_close( numbath )
+         !
+         DO jk = 1, 5
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  IF( tmask(ji,jj,jk) /= 0. ) THEN
+                     zmaskt = tmask(ji+1,jj,jk) * tmask(ji-1,jj,jk) * tmask(ji,jj+1,jk)    &
+                        &                       * tmask(ji,jj-1,jk) * tmask(ji,jj,jk+1)
+                     IF( zmaskt == 0. )   cmask(ji,jj,jk ) = 0.1
+                  ENDIF
+               END DO
+            END DO
+         END DO
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  expide   = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) )
+                  denitide = -0.9543 + 0.7662 * LOG( expide ) - 0.235 * LOG( expide )**2
+                  cmask(ji,jj,jk) = cmask(ji,jj,jk) * MIN( 1., EXP( denitide ) / 0.5 )
+               END DO
+            END DO
+         END DO
+      ELSE
+         cmask(:,:,:) = 0.e0
+      ENDIF
+
+      CALL lbc_lnk( cmask , 'T', 1. )      ! Lateral boundary conditions on cmask   (sign unchanged)
+
+
+      ! total atmospheric supply of Si
+      ! ------------------------------
+      sumdepsi = 0.e0
+      DO jm = 1, jpmois
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               sumdepsi = sumdepsi + zdustmo(ji,jj,jm) / (12.*rmoss) * 8.8        &
+                  &     * 0.075/28.1 * e1t(ji,jj) * e2t(ji,jj) * tmask(ji,jj,1) * tmask_i(ji,jj)
+            END DO
+         END DO
+      END DO
+      IF( lk_mpp )   CALL mpp_sum( sumdepsi )  ! sum over the global domain
+
+      ! N/P and Si releases due to coastal rivers
+      ! -----------------------------------------
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zcoef = raass * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)
+            cotdep(ji,jj) =  river(ji,jj)                  *1E9 / ( 12. * zcoef + rtrn ) * tmask(ji,jj,1)
+            rivinp(ji,jj) = (river(ji,jj)+riverdoc(ji,jj)) *1E9 / ( 31.6* zcoef + rtrn ) * tmask(ji,jj,1)
+            nitdep(ji,jj) = 7.6 * ndepo(ji,jj)                  / ( 14E6*raass*fse3t(ji,jj,1) + rtrn ) * tmask(ji,jj,1)
+         END DO
+      END DO
+      ! Lateral boundary conditions on ( cotdep, rivinp, nitdep )   (sign unchanged)
+      CALL lbc_lnk( cotdep , 'T', 1. )  ;  CALL lbc_lnk( rivinp , 'T', 1. )  ;  CALL lbc_lnk( nitdep , 'T', 1. )
+
+      rivpo4input=0.e0
+      rivalkinput=0.e0
+      nitdepinput=0.e0
+      DO jj = 2 , jpjm1
+         DO ji = 2, jpim1
+            zcoef = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1) * tmask(ji,jj,1) * tmask_i(ji,jj) * raass
+            rivpo4input = rivpo4input + rivinp(ji,jj) * zcoef
+            rivalkinput = rivalkinput + cotdep(ji,jj) * zcoef
+            nitdepinput = nitdepinput + nitdep(ji,jj) * zcoef
+         END DO
+     END DO
+      IF( lk_mpp ) THEN
+         CALL mpp_sum( rivpo4input )  ! sum over the global domain
+         CALL mpp_sum( rivalkinput )  ! sum over the global domain
+         CALL mpp_sum( nitdepinput )  ! sum over the global domain
+      ENDIF
+
+
+      ! Coastal supply of iron
+      ! -------------------------
+      DO jk = 1, jpkm1
+         ironsed(:,:,jk) = sedfeinput * cmask(:,:,jk) / ( fse3t(:,:,jk) * rjjss )
+      END DO
+      CALL lbc_lnk( ironsed , 'T', 1. )      ! Lateral boundary conditions on ( ironsed )   (sign unchanged)
+
+
+   END SUBROUTINE p4z_sed_init
+
+
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No PISCES bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE p4z_sed                         ! Empty routine
+   END SUBROUTINE p4z_sed
+#endif 
+
+   !!======================================================================
+END MODULE  p4zsed
Index: trunk/NEMO/TOP_SRC/PISCES/p4zsink.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zsink.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zsink.F90	(revision 935)
@@ -0,0 +1,627 @@
+MODULE p4zsink
+   !!======================================================================
+   !!                         ***  MODULE p4zsink  ***
+   !! TOP :   PISCES Compute vertical flux of particulate matter due to gravitational sinking
+   !!======================================================================
+   !! History :   1.0  !  2004     (O. Aumont) Original code
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   p4z_sink       :  Compute vertical flux of particulate matter due to gravitational sinking
+   !!----------------------------------------------------------------------
+   USE trc
+   USE oce_trc         !
+   USE trp_trc
+   USE sms
+   USE p4zsink2        !
+   USE prtctl_trc
+
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   p4z_sink    ! called in p4zbio.F90
+
+   !! * Shared module variables
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !:
+     wsbio3, wsbio4,      &    !: POC and GOC sinking speeds
+     wscal                     !: Calcite and BSi sinking speeds
+
+   !! * Module variables
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &   !:
+     sinking, sinking2,   &    !: POC sinking fluxes (different meanings depending on the parameterization
+     sinkcal, sinksil,    &    !: CaCO3 and BSi sinking fluxes
+     sinkfer                   !: Small BFe sinking flux
+
+#if  defined key_kriest
+   REAL(wp)          ::       &   
+      xkr_sfact    = 250.  ,  &   !: Sinking factor
+      xkr_stick    = 0.2   ,  &   !: Stickiness
+      xkr_nnano    = 2.337 ,  &   !: Nbr of cell in nano size class
+      xkr_ndiat    = 3.718 ,  &   !: Nbr of cell in diatoms size class
+      xkr_nmeso    = 7.147 ,  &   !: Nbr of cell in mesozoo  size class
+      xkr_naggr    = 9.877        !: Nbr of cell in aggregates  size class
+
+   REAL(wp)          ::       &   
+      xkr_frac
+
+   REAL(wp), PUBLIC ::        &
+      xkr_dnano            ,  &   !: Size of particles in nano pool
+      xkr_ddiat            ,  &   !: Size of particles in diatoms pool
+      xkr_dmeso            ,  &   !: Size of particles in mesozoo pool
+      xkr_daggr            ,  &   !: Size of particles in aggregates pool
+      xkr_wsbio_min        ,  &   !: min vertical particle speed
+      xkr_wsbio_max               !: max vertical particle speed
+
+   REAL(wp), PUBLIC, DIMENSION(jpk) ::   &   !:
+      xnumm                       !:     maximum number of particles in aggregates
+
+#endif
+
+#if ! defined key_kriest
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &   !:
+     sinkfer2                  !: Big Fe sinking flux
+#endif 
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Header:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+#if defined key_kriest
+
+   SUBROUTINE p4z_sink ( kt, jnt )
+      !!---------------------------------------------------------------------
+      !!                ***  ROUTINE p4z_sink  ***
+      !!
+      !! ** Purpose :   Compute vertical flux of particulate matter due to
+      !!              gravitational sinking - Kriest parameterization
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+
+      INTEGER, INTENT(in) :: kt, jnt
+      INTEGER  :: ji, jj, jk
+      INTEGER  :: iksed
+      REAL(wp) :: zagg1, zagg2, zagg3, zagg4, zagg5, zaggsi, zaggsh
+      REAL(wp) :: zagg , zaggdoc, znumdoc
+      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
+#endif
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   znum3d
+      CHARACTER (len=25) :: charout
+
+      !!---------------------------------------------------------------------
+
+      IF( ( kt * jnt ) == nittrc000  )   CALL p4z_sink_init   ! Initialization (first time-step only)
+
+       zstep = rfact2 / rjjss      ! Time step duration for biology
+
+
+!     Initialisation of variables used to compute Sinking Speed
+!     ---------------------------------------------------------
+
+       znum3d(:,:,:) = 0.e0
+       iksed = 10
+       zval1 = 1. + xkr_zeta
+       zval2 = 1. + xkr_zeta + xkr_eta
+       zval3 = 1. + xkr_eta
+
+!     Computation of the vertical sinking speed : Kriest et Evans, 2000
+!     -----------------------------------------------------------------
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
+                  znum = trn(ji,jj,jk,jppoc) / ( trn(ji,jj,jk,jpnum) + rtrn ) / xkr_massp
+! -------------- To avoid sinking speed over 50 m/day -------
+                  znum  = MIN( xnumm(jk), znum )
+                  znum  = MAX( 1.1      , znum )
+                  znum3d(ji,jj,jk) = znum
+!------------------------------------------------------------
+                  zeps  = ( zval1 * znum - 1. )/ ( znum - 1. )
+                  zfm   = xkr_frac**( 1. - zeps )
+                  zgm   = xkr_frac**( zval1 - zeps )
+                  zdiv  = MAX( 1.e-4, ABS( zeps - zval2 ) ) * SIGN( 1., ( zeps - zval2 ) )
+                  zdiv1 = zeps - zval3
+                  wsbio3(ji,jj,jk) = xkr_wsbio_min * ( zeps - zval1 ) / zdiv    &
+     &                             - xkr_wsbio_max *   zgm * xkr_eta  / zdiv
+                  wsbio4(ji,jj,jk) = xkr_wsbio_min *   ( zeps-1. )    / zdiv1   &
+     &                             - xkr_wsbio_max *   zfm * xkr_eta  / zdiv1
+                  IF( znum == 1.1)   wsbio3(ji,jj,jk) = wsbio4(ji,jj,jk)
+               ENDIF
+            END DO
+         END DO
+      END DO
+
+      wscal(:,:,:) = MAX( wsbio3(:,:,:), 50. )
+
+
+!   INITIALIZE TO ZERO ALL THE SINKING ARRAYS
+!   -----------------------------------------
+
+      sinking (:,:,:) = 0.e0
+      sinking2(:,:,:) = 0.e0
+      sinkcal (:,:,:) = 0.e0
+      sinkfer (:,:,:) = 0.e0
+      sinksil (:,:,:) = 0.e0
+
+!   Compute the sedimentation term using p4zsink2 for all
+!   the sinking particles
+!   -----------------------------------------------------
+
+      CALL p4z_sink2( wsbio3, sinking , jppoc )
+      CALL p4z_sink2( wsbio4, sinking2, jpnum )
+      CALL p4z_sink2( wsbio3, sinkfer , jpsfe )
+      CALL p4z_sink2( wscal , sinksil , jpdsi )
+      CALL p4z_sink2( wscal , sinkcal , jpcal )
+
+!  Exchange between organic matter compartments due to
+!  coagulation/disaggregation
+!  ---------------------------------------------------
+
+      zval1 = 1. + xkr_zeta
+      zval2 = 1. + xkr_eta
+      zval3 = 3. + xkr_eta
+      zval4 = 4. + xkr_eta
+
+      DO jk = 1,jpkm1
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
+
+                  znum = trn(ji,jj,jk,jppoc)/(trn(ji,jj,jk,jpnum)+rtrn) / xkr_massp
+! -------------- To avoid sinking speed over 50 m/day -------
+                  znum  = min(xnumm(jk),znum)
+                  znum  = MAX( 1.1,znum)
+!------------------------------------------------------------
+                  zeps  = ( zval1 * znum - 1.) / ( znum - 1.)
+                  zdiv  = MAX( 1.e-4, ABS( zeps - zval3) ) * SIGN( 1., zeps - zval3 )
+                  zdiv1 = MAX( 1.e-4, ABS( zeps - 4.   ) ) * SIGN( 1., zeps - 4.    )
+                  zdiv2 = zeps - 2.
+                  zdiv3 = zeps - 3.
+                  zdiv4 = zeps - zval2
+                  zdiv5 = 2.* zeps - zval4
+                  zfm   = xkr_frac**( 1.- zeps )
+                  zsm   = xkr_frac**xkr_eta
+
+!    Part I : Coagulation dependant on turbulence
+!    ----------------------------------------------
+
+                  zagg1 = ( 0.163 * trn(ji,jj,jk,jpnum)**2               &
+                     &            * 2.*( (zfm-1.)*(zfm*xkr_mass_max**3-xkr_mass_min**3)    &
+                     &            * (zeps-1)/zdiv1 + 3.*(zfm*xkr_mass_max-xkr_mass_min)    &
+                     &            * (zfm*xkr_mass_max**2-xkr_mass_min**2)                  &
+                     &            * (zeps-1.)**2/(zdiv2*zdiv3))            &
+# if defined key_off_degrad
+                     &                 *facvol(ji,jj,jk)       &
+# endif
+                     &    )
+
+                  zagg2 = (  2*0.163*trn(ji,jj,jk,jpnum)**2*zfm*                       &
+                     &                   ((xkr_mass_max**3+3.*(xkr_mass_max**2          &
+                     &                    *xkr_mass_min*(zeps-1.)/zdiv2                 &
+                     &                    +xkr_mass_max*xkr_mass_min**2*(zeps-1.)/zdiv3)    &
+                     &                    +xkr_mass_min**3*(zeps-1)/zdiv1)                  &
+                     &                    -zfm*xkr_mass_max**3*(1.+3.*((zeps-1.)/           &
+                     &                    (zeps-2.)+(zeps-1.)/zdiv3)+(zeps-1.)/zdiv1))      &
+#    if defined key_off_degrad
+                     &                 *facvol(ji,jj,jk)             &
+#    endif
+                     &    )
+
+                  zagg3 = (  0.163*trn(ji,jj,jk,jpnum)**2*zfm**2*8. * xkr_mass_max**3   &
+#    if defined key_off_degrad
+                     &                 *facvol(ji,jj,jk)             &
+#    endif
+                     &    )
+
+                  zaggsh = ( zagg1 + zagg2 + zagg3 ) * rfact2 * xdiss(ji,jj,jk) / 1000.
+
+!    Aggregation of small into large particles
+!    Part II : Differential settling
+!    ----------------------------------------------
+
+                  zagg4 = (  2.*3.141*0.125*trn(ji,jj,jk,jpnum)**2*                       &
+                     &                 xkr_wsbio_min*(zeps-1.)**2                         &
+                     &                 *(xkr_mass_min**2*((1.-zsm*zfm)/(zdiv3*zdiv4)      &
+                     &                 -(1.-zfm)/(zdiv*(zeps-1.)))-                       &
+                     &                 ((zfm*zfm*xkr_mass_max**2*zsm-xkr_mass_min**2)     &
+                     &                 *xkr_eta)/(zdiv*zdiv3*zdiv5) )                     &
+# if defined key_off_degrad
+                     &                 *facvol(ji,jj,jk)        &
+# endif
+                     &    )
+
+                  zagg5 = (  2.*3.141*0.125*trn(ji,jj,jk,jpnum)**2                         &
+                     &                 *(zeps-1.)*zfm*xkr_wsbio_min                        &
+                     &                 *(zsm*(xkr_mass_min**2-zfm*xkr_mass_max**2)         &
+                     &                 /zdiv3-(xkr_mass_min**2-zfm*zsm*xkr_mass_max**2)    &
+                     &                 /zdiv)                   &
+# if defined key_off_degrad
+                     &                 *facvol(ji,jj,jk)        &
+# endif
+                     &    )
+
+                  zaggsi = ( zagg4 + zagg5 ) * zstep / 10.
+
+                  zagg = 0.5 * xkr_stick * ( zaggsh + zaggsi )
+
+!     Aggregation of DOC to small particles
+!     --------------------------------------
+
+                  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
+            END DO
+         END DO
+      END DO
+
+#    if defined key_trc_dia3d
+      zrfact2 = 1.e3 * rfact2r
+      trc2d(:,:, 5)   = sinking (:,:,iksed+1) * zrfact2
+      trc2d(:,:, 6)   = sinking2(:,:,iksed+1) * zrfact2
+      trc2d(:,:, 7)   = sinkfer (:,:,iksed+1) * zrfact2
+      trc2d(:,:, 9)   = sinksil (:,:,iksed+1) * zrfact2
+      trc2d(:,:,10)   = sinkcal (:,:,iksed+1) * zrfact2
+      trc3d(:,:,:,12) = sinking (:,:,:)       * zrfact2
+      trc3d(:,:,:,13) = sinking2(:,:,:)       * zrfact2
+      trc3d(:,:,:,14) = sinksil (:,:,:)       * zrfact2
+      trc3d(:,:,:,15) = sinkcal (:,:,:)       * zrfact2
+      trc3d(:,:,:,16) = znum3d  (:,:,:)
+      trc3d(:,:,:,17) = wsbio3  (:,:,:)
+      trc3d(:,:,:,18) = wsbio4  (:,:,:)
+#    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
+
+   SUBROUTINE p4z_sink_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE p4z_sink_init  ***
+      !!
+      !! ** Purpose :   Initialization of sinking parameters
+      !!                Kriest parameterization only
+      !!
+      !! ** Method  :   Read the natkriest namelist and check the parameters
+      !!      called at the first timestep (nittrc000)
+      !!
+      !! ** input   :   Namelist natkriest
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER  ::   jk, jn, kiter
+      REAL(wp) ::   znum, zdiv
+      REAL(wp) ::   zws, zwr, zwl,wmax, znummax
+      REAL(wp) ::   zmin, zmax, zl, zr, xacc
+
+      NAMELIST/natkrsize/ xkr_sfact, xkr_stick ,  &
+         &                xkr_nnano, xkr_ndiat, xkr_nmeso, xkr_naggr
+
+      !!----------------------------------------------------------------------
+      !                               ! natkriest : kriest parameters
+      !                               ! -----------------------------
+      REWIND( numnat )                     ! read natkriest
+      READ  ( numnat, natkrsize )
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) ' Namelist : natkrsize'
+         WRITE(numout,*) '    Sinking factor                           xkr_sfact    = ', xkr_sfact
+         WRITE(numout,*) '    Stickiness                               xkr_stick    = ', xkr_stick
+         WRITE(numout,*) '    Nbr of cell in nano size class           xkr_nnano    = ', xkr_nnano
+         WRITE(numout,*) '    Nbr of cell in diatoms size class        xkr_ndiat    = ', xkr_ndiat
+         WRITE(numout,*) '    Nbr of cell in mesozoo size class        xkr_nmeso    = ', xkr_nmeso
+         WRITE(numout,*) '    Nbr of cell in aggregates size class     xkr_naggr    = ', xkr_naggr
+     ENDIF
+
+
+     ! max and min vertical particle speed
+     xkr_wsbio_min = xkr_sfact * xkr_mass_min**xkr_eta
+     xkr_wsbio_max = xkr_sfact * xkr_mass_max**xkr_eta
+     WRITE(numout,*) ' max and min vertical particle speed ', xkr_wsbio_min, xkr_wsbio_max
+
+     !
+     !    effect of the sizes of the different living pools on particle numbers
+     !    nano = 2um-20um -> mean size=6.32 um -> ws=2.596 -> xnum=xnnano=2.337
+     !    diat and microzoo = 10um-200um -> 44.7 -> 8.732 -> xnum=xndiat=3.718
+     !    mesozoo = 200um-2mm -> 632.45 -> 45.14 -> xnum=xnmeso=7.147
+     !    aggregates = 200um-10mm -> 1414 -> 74.34 -> xnum=xnaggr=9.877
+     !    doc aggregates = 1um
+     ! ----------------------------------------------------------
+
+     xkr_dnano = 1. / ( xkr_massp * xkr_nnano )
+     xkr_ddiat = 1. / ( xkr_massp * xkr_ndiat )
+     xkr_dmeso = 1. / ( xkr_massp * xkr_nmeso )
+     xkr_daggr = 1. / ( xkr_massp * xkr_naggr )
+
+      !!---------------------------------------------------------------------
+      !!    'key_kriest'                                                  ???
+      !!---------------------------------------------------------------------
+      !  COMPUTATION OF THE VERTICAL PROFILE OF MAXIMUM SINKING SPEED
+      !  Search of the maximum number of particles in aggregates for each k-level.
+      !  Bissection Method
+      !--------------------------------------------------------------------
+      WRITE(numout,*)
+      WRITE(numout,*)'    kriest : Compute maximum number of particles in aggregates'
+
+      xacc     =  0.001
+      kiter    = 50
+      zmin     =  1.10
+      zmax     = xkr_mass_max / xkr_mass_min
+      xkr_frac = zmax
+
+      DO jk = 1,jpk
+         zl = zmin
+         zr = zmax
+         wmax = 0.5 * fse3t(1,1,jk) * rjjss / rfact2
+         zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
+         znum = zl - 1.
+         zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
+            & - ( xkr_wsbio_max * xkr_eta * znum * &
+            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
+            & - wmax
+
+         zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
+         znum = zr - 1.
+         zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
+            & - ( xkr_wsbio_max * xkr_eta * znum * &
+            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
+            & - wmax
+iflag:  DO jn = 1, kiter
+           IF( zwl == 0.e0 ) THEN
+              znummax = zl
+           ELSE IF ( zwr == 0.e0 ) THEN
+              znummax = zr
+           ELSE
+              znummax = ( zr + zl ) / 2.
+              zdiv = xkr_zeta + xkr_eta - xkr_eta * znummax
+              znum = znummax - 1.
+              zws =  xkr_wsbio_min * xkr_zeta / zdiv &
+                 & - ( xkr_wsbio_max * xkr_eta * znum * &
+                 &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
+                 & - wmax
+              IF( zws * zwl < 0. ) THEN
+                 zr = znummax
+              ELSE
+                 zl = znummax
+              ENDIF
+              zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
+              znum = zl - 1.
+              zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
+                 & - ( xkr_wsbio_max * xkr_eta * znum * &
+                 &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
+                 & - wmax
+
+              zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
+              znum = zr - 1.
+              zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
+                 & - ( xkr_wsbio_max * xkr_eta * znum * &
+                 &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
+                 & - wmax
+
+              IF ( ABS ( zws )  <= xacc ) EXIT iflag
+
+           ENDIF
+
+        END DO iflag
+
+        xnumm(jk) = znummax
+        WRITE(numout,*) '       jk = ', jk, ' wmax = ', wmax,' xnum max = ', xnumm(jk)
+
+     END DO
+
+  END SUBROUTINE p4z_sink_init
+
+#else
+
+   SUBROUTINE p4z_sink ( kt, jnt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_sink  ***
+      !!
+      !! ** Purpose :   Compute vertical flux of particulate matter due to 
+      !!                gravitational sinking
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kt, jnt
+      INTEGER  ::   ji, jj, jk
+      INTEGER  ::   iksed
+      REAL(wp) ::   zagg1, zagg2, zagg3, zagg4
+      REAL(wp) ::   zagg , zaggfe, zaggdoc, zaggdoc2
+      REAL(wp) ::   zfact, zstep, zwsmax
+#if defined key_trc_dia3d
+      REAL(wp) ::   zrfact2
+#endif
+      CHARACTER (len=25) :: charout
+      !!---------------------------------------------------------------------
+
+       zstep = rfact2 / rjjss      ! Timestep duration for biology
+
+
+!    Sinking speeds of detritus is increased with depth as shown
+!    by data and from the coagulation theory
+!    -----------------------------------------------------------
+
+      iksed = 10
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji=1,jpi
+               zfact = MAX( 0., fsdepw(ji,jj,jk+1)-hmld(ji,jj) ) / 4000.
+               wsbio4(ji,jj,jk) = wsbio2 + ( 200.- wsbio2 ) * zfact
+            END DO
+         END DO
+      END DO
+
+!      LIMIT THE VALUES OF THE SINKING SPEEDS 
+!      TO AVOID NUMERICAL INSTABILITIES
+
+      wsbio3(:,:,:) = wsbio
+!
+! OA Below, this is garbage. the ideal would be to find a time-splitting
+! OA algorithm that does not increase the computing cost by too much
+! OA In ROMS, I have included a time-splitting procedure. But it is 
+! OA too expensive as the loop is computed globally. Thus, a small e3t
+! OA at one place determines the number of subtimesteps globally
+! OA AWFULLY EXPENSIVE !! Not able to find a better approach. Damned !!
+
+      DO jk = 1,jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zwsmax = 0.8 * fse3t(ji,jj,jk) / zstep
+               wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax )
+               wsbio3(ji,jj,jk) = MIN( wsbio3(ji,jj,jk), zwsmax )
+            END DO
+         END DO
+      END DO
+
+      wscal(:,:,:) = wsbio4(:,:,:)
+
+!   INITIALIZE TO ZERO ALL THE SINKING ARRAYS
+!   -----------------------------------------
+
+      sinking (:,:,:) = 0.e0
+      sinking2(:,:,:) = 0.e0
+      sinkcal (:,:,:) = 0.e0
+      sinkfer (:,:,:) = 0.e0
+      sinksil (:,:,:) = 0.e0
+      sinkfer2(:,:,:) = 0.e0
+
+!   Compute the sedimentation term using p4zsink2 for all
+!   the sinking particles
+!   -----------------------------------------------------
+
+      CALL p4z_sink2( wsbio3, sinking , jppoc )
+      CALL p4z_sink2( wsbio3, sinkfer , jpsfe )
+      CALL p4z_sink2( wsbio4, sinking2, jpgoc )
+      CALL p4z_sink2( wsbio4, sinkfer2, jpbfe )
+      CALL p4z_sink2( wsbio4, sinksil , jpdsi )
+      CALL p4z_sink2( wscal , sinkcal , jpcal )
+
+!  Exchange between organic matter compartments due to
+!  coagulation/disaggregation
+!  ---------------------------------------------------
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zfact = zstep * xdiss(ji,jj,jk)
+
+!    Part I : Coagulation dependent on turbulence
+!    ----------------------------------------------
+
+# if defined key_off_degrad
+               zagg1 = 940.* zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc) * facvol(ji,jj,jk)
+# else
+               zagg1 = 940.* zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
+# endif
+
+# if defined key_off_degrad
+               zagg2 = 1.054e4 * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc) * facvol(ji,jj,jk)
+# else
+               zagg2 = 1.054e4 * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
+# endif
+
+!    Aggregation of small into large particles
+!    Part II : Differential settling
+!    ----------------------------------------------
+
+# if defined key_off_degrad
+               zagg3 = 0.66 * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc) * facvol(ji,jj,jk)
+# else
+               zagg3 = 0.66 * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
+# endif
+
+# if defined key_off_degrad
+               zagg4 = 0.e0 * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc) * facvol(ji,jj,jk)
+# else
+               zagg4 = 0.e0 * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
+# endif
+
+               zagg   = zagg1 + zagg2 + zagg3 + zagg4
+               zaggfe = zagg * trn(ji,jj,jk,jpsfe) / ( trn(ji,jj,jk,jppoc) + rtrn )
+
+!     Aggregation of DOC to small particles
+!     --------------------------------------
+
+               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)
+!
+!  Update the trends
+!  -----------------
+!
+               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zagg + zaggdoc
+               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zagg + zaggdoc2
+               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
+         END DO
+      END DO
+
+# if defined key_trc_dia3d
+      zrfact2 = 1.e3 * rfact2r
+      trc2d(:,:, 5) = sinking (:,:,iksed+1) * zrfact2
+      trc2d(:,:, 6) = sinking2(:,:,iksed+1) * zrfact2
+      trc2d(:,:, 7) = sinkfer (:,:,iksed+1) * zrfact2
+      trc2d(:,:, 8) = sinkfer2(:,:,iksed+1) * zrfact2
+      trc2d(:,:, 9) = sinksil (:,:,iksed+1) * zrfact2
+      trc2d(:,:,10) = sinkcal (:,:,iksed+1) * zrfact2
+# 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
+
+#endif
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No PISCES bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE p4z_sink                    ! Empty routine
+   END SUBROUTINE p4z_sink
+#endif 
+
+   !!======================================================================
+END MODULE  p4zsink
Index: trunk/NEMO/TOP_SRC/PISCES/p4zsink2.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/p4zsink2.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/p4zsink2.F90	(revision 935)
@@ -0,0 +1,167 @@
+MODULE p4zsink2
+   !!======================================================================
+   !!                         ***  MODULE p4zsink2  ***
+   !! TOP :   PISCES Compute vertical flux of particulate matter due to gravitational sinking
+   !!======================================================================
+   !! History :   1.0  !  2004     (O. Aumont) Original code
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'                                       PISCES bio-model
+   !!----------------------------------------------------------------------
+   !!   p4z_sink2       :  Compute vertical flux of particulate matter due to gravitational sinking
+   !!----------------------------------------------------------------------
+   USE oce_trc         !
+   USE trp_trc
+   USE sms
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   p4z_sink2    ! called in p4zbio.F90
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Header:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE p4z_sink2( wstmp, sinktemp, jn )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_sink2  ***
+      !!
+      !! ** Purpose :   Compute the sedimentation terms for the various sinking
+      !!     particles. The scheme used to compute the trends is based
+      !!     on MUSCL.
+      !!
+      !! ** Method  : - this ROUTINE compute not exactly the advection but the
+      !!      transport term, i.e.  div(u*tra).
+      !!---------------------------------------------------------------------
+      INTEGER , INTENT(in   )                         ::   jn         ! tracer index index      
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   wstmp      ! ???
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   sinktemp   ! ???
+      !!
+      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  , zwstmp2
+      !!---------------------------------------------------------------------
+
+      zstep  = rfact2 / 2.
+
+      ztraz(:,:,:) = 0.e0
+      zkz  (:,:,:) = 0.e0
+      zakz (:,:,:) = 0.e0
+
+      DO jk = 1, jpkm1
+# if defined key_off_degrad
+         zwstmp2(:,:,jk+1)=-wstmp(:,:,jk)/rjjss*tmask(:,:,jk+1)*facvol(:,:,jk)
+# else
+         zwstmp2(:,:,jk+1)=-wstmp(:,:,jk)/rjjss*tmask(:,:,jk+1)
+
+# endif
+      END DO
+ 
+      zwstmp2(:,:,1) = 0.e0
+!
+! Vertical advective flux
+!-------------------------------
+
+      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.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
+         !
+         ! 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
+         
+         ! 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
+               ! 
+               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(:,:,:)
+
+      !
+   END SUBROUTINE p4z_sink2
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No PISCES bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE p4z_sink2( wstmp, sinktemp, jn )         ! Empty routine
+      INTEGER, INTENT( in ) ::   jn
+      REAL   , INTENT( in ) ::   wstmp,sinktemp    
+      WRITE(*,*) 'p4z_sink2: You should not have seen this print! error?', jn, wstmp, sinktemp
+   END SUBROUTINE p4z_sink2
+#endif 
+
+   !!======================================================================
+END MODULE  p4zsink2
Index: trunk/NEMO/TOP_SRC/PISCES/par_pisces.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/par_pisces.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/par_pisces.F90	(revision 935)
@@ -0,0 +1,113 @@
+MODULE par_pisces
+   !!======================================================================
+   !!                        ***  par_pisces  ***
+   !! TOP :   set the PISCES parameters
+   !!======================================================================
+   !! History :   2.0  !  2007-12  (C. Ethe, G. Madec)  revised architecture
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Id:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+   USE par_lobster, ONLY : jp_lobster   !: number of tracers in LOBSTER
+
+   IMPLICIT NONE
+   PUBLIC
+
+   INTEGER, PUBLIC, PARAMETER ::   jp_l  = jp_lobster      !: cumulative number of already defined TRC
+
+#if defined key_pisces  &&  defined key_kriest
+   !!---------------------------------------------------------------------
+   !!   'key_pisces' & 'key_kriest'                 PISCES bio-model + ???
+   !!---------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_pisces     = .TRUE.  !: PISCES flag 
+   LOGICAL, PUBLIC, PARAMETER ::   lk_kriest     = .TRUE.  !: Kriest flag 
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces     =  23     !: number of passive tracers
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_2d  =  13     !: additional 2d output ('key_trc_diaadd')
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_3d  =  18     !: additional 3d output ('key_trc_diaadd')
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_trd =   1     !: number of sms trends for PISCES
+
+   ! assign an index in trc arrays for each LOBSTER prognostic variables
+   !    WARNING: be carefull about the order when reading the restart
+        !   !!gm  this warning should be obsolet with IOM
+   INTEGER, PUBLIC, PARAMETER ::   jpdic = jp_l +  1    !: dissolved inoganic carbon concentration 
+   INTEGER, PUBLIC, PARAMETER ::   jptal = jp_l +  2    !: total alkalinity 
+   INTEGER, PUBLIC, PARAMETER ::   jpoxy = jp_l +  3    !: oxygen carbon concentration 
+   INTEGER, PUBLIC, PARAMETER ::   jpcal = jp_l +  4    !: calcite  concentration 
+   INTEGER, PUBLIC, PARAMETER ::   jppo4 = jp_l +  5    !: phosphate concentration 
+   INTEGER, PUBLIC, PARAMETER ::   jppoc = jp_l +  6    !: small particulate organic phosphate concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpsil = jp_l +  7    !: silicate concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpphy = jp_l +  8    !: phytoplancton concentration 
+   INTEGER, PUBLIC, PARAMETER ::   jpzoo = jp_l +  9    !: zooplancton concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpdoc = jp_l + 10    !: dissolved organic carbon concentration 
+   INTEGER, PUBLIC, PARAMETER ::   jpdia = jp_l + 11    !: Diatoms Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpmes = jp_l + 12    !: Mesozooplankton Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpbsi = jp_l + 13    !: (big) Silicate Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpfer = jp_l + 14    !: Iron Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpnum = jp_l + 15    !: Big iron particles Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpsfe = jp_l + 16    !: number of particulate organic phosphate concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpdfe = jp_l + 17    !: Diatoms iron Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpdsi = jp_l + 18    !: Diatoms Silicate Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpnfe = jp_l + 19    !: Nano iron Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpnch = jp_l + 20    !: Nano Chlorophyll Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpdch = jp_l + 21    !: Diatoms Chlorophyll Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpno3 = jp_l + 22    !: Nitrates Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpnh4 = jp_l + 23    !: Ammonium Concentration
+
+#elif defined key_pisces
+   !!---------------------------------------------------------------------
+   !!   'key_pisces'   :                         standard PISCES bio-model
+   !!---------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_pisces     = .TRUE.  !: PISCES flag 
+   LOGICAL, PUBLIC, PARAMETER ::   lk_kriest     = .FALSE. !: Kriest flag 
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces     = 24      !: number of PISCES passive tracers
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_2d  = 13      !: additional 2d output ('key_trc_diaadd')
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_3d  = 11      !: additional 3d output ('key_trc_diaadd')
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_trd =  1      !: number of sms trends for PISCES
+
+   ! assign an index in trc arrays for each LOBSTER prognostic variables
+   !    WARNING: be carefull about the order when reading the restart
+        !   !!gm  this warning should be obsolet with IOM
+   INTEGER, PUBLIC, PARAMETER ::   jpdic = jp_l +  1    !: dissolved inoganic carbon concentration 
+   INTEGER, PUBLIC, PARAMETER ::   jptal = jp_l +  2    !: total alkalinity 
+   INTEGER, PUBLIC, PARAMETER ::   jpoxy = jp_l +  3    !: oxygen carbon concentration 
+   INTEGER, PUBLIC, PARAMETER ::   jpcal = jp_l +  4    !: calcite  concentration 
+   INTEGER, PUBLIC, PARAMETER ::   jppo4 = jp_l +  5    !: phosphate concentration 
+   INTEGER, PUBLIC, PARAMETER ::   jppoc = jp_l +  6    !: small particulate organic phosphate concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpsil = jp_l +  7    !: silicate concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpphy = jp_l +  8    !: phytoplancton concentration 
+   INTEGER, PUBLIC, PARAMETER ::   jpzoo = jp_l +  9    !: zooplancton concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpdoc = jp_l + 10    !: dissolved organic carbon concentration 
+   INTEGER, PUBLIC, PARAMETER ::   jpdia = jp_l + 11    !: Diatoms Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpmes = jp_l + 12    !: Mesozooplankton Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpbsi = jp_l + 13    !: (big) Silicate Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpfer = jp_l + 14    !: Iron Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpbfe = jp_l + 15    !: Big iron particles Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpgoc = jp_l + 16    !: big particulate organic phosphate concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpsfe = jp_l + 17    !: Small iron particles Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpdfe = jp_l + 18    !: Diatoms iron Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpdsi = jp_l + 19    !: Diatoms Silicate Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpnfe = jp_l + 20    !: Nano iron Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpnch = jp_l + 21    !: Nano Chlorophyll Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpdch = jp_l + 22    !: Diatoms Chlorophyll Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpno3 = jp_l + 23    !: Nitrates Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpnh4 = jp_l + 24    !: Ammonium Concentration
+
+#else
+   !!---------------------------------------------------------------------
+   !!   Default                                   No CFC geochemical model
+   !!---------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_pisces     = .FALSE.  !: CFC flag 
+   LOGICAL, PUBLIC, PARAMETER ::   lk_kriest     = .FALSE.  !: Kriest flag 
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces     =  0       !: No CFC tracers
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_2d  =  0       !: No CFC additional 2d output arrays 
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_3d  =  0       !: No CFC additional 3d output arrays 
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_trd =  0       !: number of sms trends for PISCES
+#endif
+
+   ! Starting/ending PISCES do-loop indices (N.B. no PISCES : jpl_pcs < jpf_pcs the do-loop are never done)
+   INTEGER, PUBLIC, PARAMETER ::   jp_pcs0 = jp_l + 1           !: First index of CFC passive tracers
+   INTEGER, PUBLIC, PARAMETER ::   jp_pcs1 = jp_l + jp_pisces   !: Last  index of CFC passive tracers
+
+   !!======================================================================
+END MODULE par_pisces
Index: trunk/NEMO/TOP_SRC/PISCES/sms_pisces.h90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/sms_pisces.h90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/sms_pisces.h90	(revision 935)
@@ -0,0 +1,60 @@
+   !!----------------------------------------------------------------------
+   !!                     ***  sms_pisces.h90  ***  
+   !! TOP :   PISCES Source Minus Sink variables
+   !!----------------------------------------------------------------------
+   !! History :   1.0  !  2000-02 (O. Aumont) original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 1.0 , LOCEAN-IPSL (2005) 
+   !! $Header$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !! Variable for chemistry of the CO2 cycle
+   !! ---------------------------------------------------------------------
+   !
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   akb3, ak13, ak23, aksp, akw3             !: ???
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   hi, borat                                !: ???
+
+   !!----------------------------------------------------------------------
+   !!  Time variables
+   !! ---------------------------------------------------------------------
+   INTEGER  ::   nrdttrc, ndayflxtr       !: ???
+   REAL(wp) ::   rfact, rfactr, rfact2, rfact2r   !: ???
+
+   !!---------------------------------------
+   !!  Biological parameters 
+   !! --------------------------------------
+   !
+   REAL(wp) ::   part, rno3, o2ut, po4r, rdenit, o2nit     !: ???
+   REAL(wp) ::   wsbio, wsbio2, xkmort, ferat3                     !: ???
+
+   !!---------------------------------------------
+   !!  Biological fluxes for light
+   !!---------------------------------------------
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   etot3                 !: ???
+
+!!----------------------------------------------------------
+!!  Biological fluxes for primary production
+!!----------------------------------------------------------
+   REAL(wp), DIMENSION(jpi,jpj) :: xksi, xksimax
+   !
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   xnanono3, xdiatno3, xnanonh4, xdiatnh4       !: ???
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   xlimphy, xlimdia, concdfe, concnfe !: ???
+
+   !!---------------------------------------------
+   !!  SMS for the organic matter
+   !!---------------------------------------------
+   !
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   xfracal, nitrfac, xlimbac, xdiss                  !: ??
+
+#if defined key_kriest
+   !!---------------------------------------------------------
+   !!  Kriest parameter for aggregation
+   !!---------------------------------------------------------
+   REAL(wp) ::  xkr_eta, xkr_zeta
+   REAL(wp) ::  xkr_mass_min, xkr_mass_max, xkr_massp
+#endif
+
Index: trunk/NEMO/TOP_SRC/PISCES/trcctl.pisces.h90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/trcctl.pisces.h90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/trcctl.pisces.h90	(revision 935)
@@ -0,0 +1,33 @@
+   !!======================================================================
+   !!                         ***  trcctl.cfc.h90  ***
+   !! TOP :   Control of CFC chemical model
+   !!======================================================================
+   !!----------------------------------------------------------------------
+   !! History :   1.0  !  2005-10 (C. Ethe) assign a parameter to name individual tracers
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 1.0 , LOCEAN-IPSL (2005) 
+   !! $Id: trcctl.pisces.h90 772 2007-12-17 11:59:33Z gm $ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' use PISCES biological model '
+
+   ! Check number of tracers
+   ! -----------------------
+#if  defined key_kriest
+      IF( jp_pisces /= 23) THEN 
+#else
+      IF( jp_pisces /= 24) THEN
+#endif
+          IF (lwp) THEN 
+              WRITE (numout,*) ' ===>>>> : w a r n i n g '
+              WRITE (numout,*) ' =======   ============= '
+              WRITE (numout,*)                               &
+              &   ' STOP, change jp_pisces',               & 
+              &   ' in parameter.passivetrc.pisces.h '  
+          END IF 
+          STOP 'TRC_CTL'
+      END IF 
+      
Index: trunk/NEMO/TOP_SRC/PISCES/trcini_pisces.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 935)
@@ -0,0 +1,96 @@
+MODULE trcini_pisces
+   !!======================================================================
+   !!                         ***  MODULE trcini_pisces  ***
+   !! TOP :   initialisation of the PISCES biochemical model
+   !!======================================================================
+   !! History :    -   !  1988-07  (E. Maier-Reiner) Original code
+   !!              -   !  1999-10  (O. Aumont, C. Le Quere)
+   !!              -   !  2002     (O. Aumont)  PISCES
+   !!             1.0  !  2005-03  (O. Aumont, A. El Moussaoui) F90
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec) from trcini.pisces.h90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'                                       PISCES bio-model
+   !!----------------------------------------------------------------------
+   !! trc_ini_pisces   : PISCES biochemical model initialisation
+   !!----------------------------------------------------------------------
+   USE par_trc         ! TOP parameters
+   USE sms             ! Source Minus Sink variables
+   USE oce_trc         ! ocean variables
+   USE trp_trc         !
+   USE p4zche 
+
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   trc_ini_pisces   ! called by trcini.F90 module
+
+
+#  include "domzgr_substitute.h90"
+#  include "top_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Id: trcini_pisces.F90 776 2007-12-19 14:10:14Z gm $ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE trc_ini_pisces
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE trc_ini_pisces ***
+      !!
+      !! ** Purpose :   Initialisation of the PISCES biochemical model
+      !!----------------------------------------------------------------------
+      INTEGER ::   jk
+
+      !!----------------------------------------------------------------------
+
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' trc_ini_pisces :   PISCES biochemical model initialisation'
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+
+
+      !                                            ! Time-step
+      rfact   = rdttra(1) * float(ndttrc)          ! ---------
+      rfactr  = 1. / rfact
+      rfact2  = rfact / float(nrdttrc)
+      rfact2r = 1. / rfact2
+
+      IF(lwp) WRITE(numout,*) '    Tracer  time step    rfact  = ', rfact, ' rdt = ', rdt
+      IF(lwp) write(numout,*) '    Biology time step    rfact2 = ', rfact2
+
+
+      !----------------------------------------------------------------------
+      ! Initialize biological variables 
+      !----------------------------------------------------------------------
+      ! Set biological ratios
+      ! ---------------------
+      rno3   = (16.+2.) / 122.
+      po4r   =   1.e0   / 122.
+      o2nit  =  32.     / 122.
+      rdenit =  97.6    /  16.
+      o2ut   = 140.     / 122.
+
+      CALL p4z_che        ! initialize the chemical constants
+
+      ndayflxtr = 0      !  Initialize a counter for the computation of chemistry
+
+      IF(lwp) WRITE(numout,*) ' Initialisation of PISCES done'
+      !
+   END SUBROUTINE trc_ini_pisces
+   
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module                            No PISCES biochemical model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_ini_pisces             ! Empty routine
+   END SUBROUTINE trc_ini_pisces
+#endif
+
+   !!======================================================================
+END MODULE trcini_pisces
Index: trunk/NEMO/TOP_SRC/PISCES/trclsm_pisces.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/trclsm_pisces.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/trclsm_pisces.F90	(revision 935)
@@ -0,0 +1,122 @@
+MODULE trclsm_pisces
+   !!======================================================================
+   !!                      ***  MODULE trclsm_lobster  ***
+   !! TOP :   initialisation of some run parameters for PISCES bio-model
+   !!======================================================================
+   !! History :    -   !  1999-10 (M.A. Foujols, M. Levy) original code
+   !!              -   !  2000-01 (L. Bopp) hamocc3, p3zd
+   !!             1.0  !  2003-08 (C. Ethe)  module F90
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec) from trclsm.pisces.h90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'   :                                   PISCES bio-model
+   !!----------------------------------------------------------------------
+   !! trc_lsm_pisces       : PISCES model namelist read
+   !!----------------------------------------------------------------------
+   USE oce_trc         ! Ocean variables
+   USE par_trc         ! TOP parameters
+   USE trc             ! TOP variables
+   USE sms             ! sms trends
+
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   trc_lsm_pisces   ! called by trclsm.F90 module
+
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Id: trclsm_pisces.F90 773 2007-12-17 16:55:32Z gm $ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE trc_lsm_pisces
+      !!----------------------------------------------------------------------
+      !!                     ***  trc_lsm_pisces  ***  
+      !!
+      !! ** Purpose :   read PISCES namelist
+      !!
+      !! ** input   :   file 'namelist.trc.sms' containing the following
+      !!             namelist: natext, natbio, natsms
+      !!                       natkriest ("key_kriest")
+      !!----------------------------------------------------------------------
+      CHARACTER (len=32) ::   clname
+      !!
+      NAMELIST/natbio/ part, nrdttrc, wsbio, xkmort, ferat3, wsbio2
+
+#if defined key_kriest
+      NAMELIST/natkrpar/xkr_eta, xkr_zeta, xkr_mass_min, xkr_mass_max
+#endif
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' trc_lsm_pisces : read PISCES namelists'
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+
+
+      !                               ! Open the namelist file
+      !                               ! ----------------------
+      clname ='namelist.trc.sms'
+      CALL ctlopn( numnat, clname, 'OLD', 'FORMATTED', 'SEQUENTIAL',   &
+         &           1, numout, .FALSE., 1 )
+
+      !                               ! natbio : biological parameters
+      !                               ! ------------------------------
+      REWIND( numnat )                     ! read natbio
+      READ  ( numnat, natbio )
+
+      IF(lwp) THEN                         ! control print
+         WRITE(numout,*) ' Namelist : natbio'
+         WRITE(numout,*) '    part of calcite not dissolved in guts     part      =', part
+         WRITE(numout,*) '    frequence pour la biologie                nrdttrc   =', nrdttrc
+         WRITE(numout,*) '    POC sinking speed                         wsbio     =', wsbio
+         WRITE(numout,*) '    half saturation constant for mortality    xkmort    =', xkmort
+         WRITE(numout,*) '    Fe/C in zooplankton                       ferat3    =', ferat3
+         WRITE(numout,*) '    Big particles sinking speed               wsbio2    =', wsbio2
+      ENDIF
+
+#if defined key_kriest
+
+      !                               ! natkrpar : kriest parameters
+      !                               ! -----------------------------
+      xkr_eta      = 0.62        
+      xkr_zeta     = 1.62        
+      xkr_mass_min = 0.0002     
+      xkr_mass_max = 1.      
+
+      REWIND( numnat )                     ! read natkriest
+      READ  ( numnat, natkrpar )
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) ' Namelist : natkrpar'
+         WRITE(numout,*) '    Sinking  exponent                        xkr_eta      = ', xkr_eta
+         WRITE(numout,*) '    N content exponent                       xkr_zeta     = ', xkr_zeta
+         WRITE(numout,*) '    Minimum mass for Aggregates              xkr_mass_min = ', xkr_mass_min
+         WRITE(numout,*) '    Maximum mass for Aggregates              xkr_mass_max = ', xkr_mass_max
+         WRITE(numout,*)
+     ENDIF
+
+
+     ! Computation of some variables
+     xkr_massp = 5.7E-6 * 7.6 * xkr_mass_min**xkr_zeta
+
+#endif
+      !
+   END SUBROUTINE trc_lsm_pisces
+
+#else
+   !!----------------------------------------------------------------------
+   !!  Dummy module :                                   No PISCES bio-model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_lsm_pisces                      ! Empty routine
+   END  SUBROUTINE  trc_lsm_pisces
+#endif  
+
+   !!======================================================================
+END MODULE trclsm_pisces
Index: trunk/NEMO/TOP_SRC/PISCES/trcsms_pisces.F90
===================================================================
--- trunk/NEMO/TOP_SRC/PISCES/trcsms_pisces.F90	(revision 935)
+++ trunk/NEMO/TOP_SRC/PISCES/trcsms_pisces.F90	(revision 935)
@@ -0,0 +1,113 @@
+MODULE trcsms_pisces
+   !!======================================================================
+   !!                         ***  MODULE trcsms_pisces  ***
+   !! TOP :   PISCES Source Minus Sink manager
+   !!======================================================================
+   !! History :   1.0  !  2004-03 (O. Aumont) Original code
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'                                       PISCES bio-model
+   !!----------------------------------------------------------------------
+   !!   trcsms_pisces        :  Time loop of passive tracers sms
+   !!----------------------------------------------------------------------
+   USE oce_trc         !
+   USE trp_trc
+   USE sms
+   USE lbclnk
+   USE lib_mpp
+   
+   USE p4zint          ! 
+   USE p4zche          ! 
+   USE p4zbio          ! 
+   USE p4zsed          ! 
+   USE p4zlys          ! 
+   USE p4zflx          ! 
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   trc_sms_pisces    ! called in trcsms.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Header:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE trc_sms_pisces( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_sms_pisces  ***
+      !!
+      !! ** Purpose :   Managment of the call to Biological sources and sinks 
+      !!              routines of PISCES bio-model
+      !!
+      !! ** Method  : - at each new day ...
+      !!              - several calls of bio and sed ???
+      !!              - ...
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index      
+      !!
+      INTEGER ::   jnt, jn
+      INTEGER ::   iyy, imm, idd
+
+
+      !!---------------------------------------------------------------------
+
+
+      iyy = ndastp/10000
+      imm = (ndastp - iyy*10000)/100
+      idd = (ndastp - iyy*10000 - imm*100)
+
+      IF( ndayflxtr /= idd) THEN      ! New days
+         !
+         ndayflxtr = idd
+
+         CALL p4z_che          ! computation of chemical constants
+
+         CALL p4z_int          ! computation of various rates for biogeochemistry
+         !
+      ENDIF
+
+
+      DO jnt = 1, nrdttrc          ! Potential time splitting if requested
+         !
+         CALL p4z_bio (kt, jnt)    ! Compute soft tissue production (POC)
+
+
+         CALL p4z_sed (kt, jnt)    ! compute soft tissue remineralisation
+
+         !
+         trb(:,:,:,:) = trn(:,:,:,:)
+         !
+      END DO
+
+      CALL p4z_lys( kt )             ! Compute CaCO3 saturation
+
+      CALL p4z_flx( kt )             ! Compute surface fluxes
+
+
+      DO jn = jp_pcs0, jp_pcs1
+        CALL lbc_lnk( trn(:,:,:,jn), 'T', 1. )
+        CALL lbc_lnk( trb(:,:,:,jn), 'T', 1. )
+        CALL lbc_lnk( tra(:,:,:,jn), 'T', 1. )
+      END DO
+      !
+   END SUBROUTINE trc_sms_pisces
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No PISCES bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE trc_sms_pisces( kt )                   ! Empty routine
+      INTEGER, INTENT( in ) ::   kt
+      WRITE(*,*) 'trc_sms_pisces: You should not have seen this print! error?', kt
+   END SUBROUTINE trc_sms_pisces
+#endif 
+
+   !!======================================================================
+END MODULE trcsms_pisces 
