CCC$Header$
CCC  TOP 1.0 , LOCEAN-IPSL (2005)
C This software is governed by CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
C ---------------------------------------------------------------------------
      SUBROUTINE p4zbio
CDIR$ LIST
#if defined key_passivetrc && defined key_trc_pisces
CCC   ------------------------------------------------------------------
CCC   
CCC   ROUTINE p4zbio : PISCES MODEL
CCC   *****************************
CCC   
CC
CC     PURPOSE.
CC     --------
CC          *P4ZBIO* ECOSYSTEM MODEL IN THE WHOLE OCEAN
CC                   THIS ROUTINE COMPUTES THE DIFFERENT INTERACTIONS
CC                   BETWEEN THE DIFFERENT COMPARTMENTS OF THE MODEL
CC     EXTERNAL :
CC     ----------
CC          p4zopt, p4zprod, p4znano, p4zdiat, p4zmicro, p4zmeso
CC          p4zsink, p4zrem
CC
CC   MODIFICATIONS:
CC   --------------
CC      original      : 2004    O. Aumont
CC ----------------------------------------------------------------
CC parameters and commons
CC ======================
CDIR$ NOLIST
      USE oce_trc
      USE trp_trc
      USE sms
      USE lib_mpp
      USE lbclnk
      IMPLICIT NONE
CDIR$ LIST
CC-----------------------------------------------------------------
CC local declarations
CC ==================
C     
      INTEGER ji, jj, jk

      REAL xcond,zdenom,zdenom1(jpi,jpj,jpk),zdenom2(jpi,jpj,jpk)
      REAL zneg, prodca
C     
      REAL irondep(jpi,jpj,jpk),sidep(jpi,jpj,jpk),prodt
      INTEGER jn
C
CC----------------------------------------------------------------------
CC statement functions
CC ===================
CDIR$ NOLIST
#include "domzgr_substitute.h90"
CDIR$ LIST
C     
C     SET HALF PRECISION CONSTANTS
C-----------------------------
C     
C     Initialisation of variables used to compute deposition
C     ------------------------------------------------------
C     
      irondep     = 0.
      sidep       = 0.
C
C
C     Iron and Si deposition at the surface
C     -------------------------------------
C
       do jj=1,jpj
         do ji=1,jpi
         irondep(ji,jj,1)=(0.01*dust(ji,jj)/(55.85*rmoss)
     &      +3E-10/raass)*rfact2/fse3t(ji,jj,1)
         sidep(ji,jj,1)=8.8*0.075*dust(ji,jj)*rfact2/
     &      (fse3t(ji,jj,1)*28.01*rmoss)
         end do
       end do
C
C     ASSIGN THE SHEAR RATE THAT IS USED FOR AGGREGATION
C     OF PHYTOPLANKTON AND DETRITUS
C
      zdiss=0.01
C
       DO jk=1,jpkm1
        DO jj=1,jpj
          DO ji=1,jpi
       if (fsdepw(ji,jj,jk+1).le.hmld(ji,jj)) zdiss(ji,jj,jk)=1.
          END DO 
        END DO
       END DO
C
C      Compute de different ratios for scavenging of iron
C      --------------------------------------------------
C
       DO jk=1,jpk
         DO jj=1,jpj
           DO ji=1,jpi
         zdenom=1./(trn(ji,jj,jk,jppoc)+trn(ji,jj,jk,jpgoc)
     $     +trn(ji,jj,jk,jpdsi)+trn(ji,jj,jk,jpcal)+rtrn)
C
         zdenom1(ji,jj,jk)=trn(ji,jj,jk,jppoc)*zdenom
         zdenom2(ji,jj,jk)=trn(ji,jj,jk,jpgoc)*zdenom
           END DO
         END DO
       END DO


C
C  Call optical routine to compute the PAR in the water column
C  -----------------------------------------------------------
C
      CALL p4zopt


C
C  Call production routine to compute phytoplankton growth rate
C  over the global ocean. Growth rates for each element is 
C  computed (C, Si, Fe, Chl)
C  ------------------------------------------------------------
C

      CALL p4zprod


C
C  Call phytoplankton mortality routines. Mortality losses for 
C  Each elements are computed (C, Fe, Si, Chl)
C  -----------------------------------------------------------
C
      CALL p4znano

      CALL p4zdiat

C
C  Call zooplankton sources/sinks routines. 
C  Each elements are computed (C, Fe, Si, Chl)
C  -----------------------------------------------------------
C
      CALL p4zmicro

      CALL p4zmeso

C
C     Call subroutine for computation of the vertical flux 
C     of particulate organic matter
C     ----------------------------------------------------
      CALL p4zsink

C
C     Call subroutine for computation of remineralization
C     terms of organic matter+scavenging of Fe
C     ----------------------------------------------------
      CALL p4zrem

C
C     Vertical loop to pre-compute concentration changes of the rapid
C     varying tracers for preventing them to fall below 0
C     ---------------------------------------------------------------
C
      DO jk = 1,jpkm1
        DO jj = 1,jpj
          DO ji = 1,jpi
C     
C     Evolution of PO4
C     ----------------
C     
         zneg = trn(ji,jj,jk,jppo4)
     &     -prorca(ji,jj,jk)-prorca2(ji,jj,jk)+denitr(ji,jj,jk)
     &     +grarem(ji,jj,jk)*sigma1+grarem2(ji,jj,jk)*sigma2
     &     +olimi(ji,jj,jk)+po4dep(ji,jj,jk)*rfact2
C     
C     Nullity test for PO4
C     --------------------
C     
         xcond=(0.5+sign(0.5,zneg)) 
         prorca(ji,jj,jk)=prorca(ji,jj,jk)*xcond
         prorca2(ji,jj,jk)=prorca2(ji,jj,jk)*xcond
         proreg(ji,jj,jk)=proreg(ji,jj,jk)*xcond
         proreg2(ji,jj,jk)=proreg2(ji,jj,jk)*xcond
         pronew(ji,jj,jk)=pronew(ji,jj,jk)*xcond
         pronew2(ji,jj,jk)=pronew2(ji,jj,jk)*xcond
C
C     Evolution of NO3
C     ----------------
C
         zneg = trn(ji,jj,jk,jpno3)
     &     -pronew(ji,jj,jk)-pronew2(ji,jj,jk)
     &     +po4dep(ji,jj,jk)*rfact2+onitr(ji,jj,jk)
     &     -denitr(ji,jj,jk)*rdenit+nitdep(ji,jj,jk)*rfact2
C
C     Nullity test for NO3
C     --------------------
C
         xcond=(0.5+sign(0.5,zneg))
         prorca(ji,jj,jk)=prorca(ji,jj,jk)*xcond
         prorca2(ji,jj,jk)=prorca2(ji,jj,jk)*xcond
         proreg(ji,jj,jk)=proreg(ji,jj,jk)*xcond
         proreg2(ji,jj,jk)=proreg2(ji,jj,jk)*xcond
         pronew(ji,jj,jk)=pronew(ji,jj,jk)*xcond
         pronew2(ji,jj,jk)=pronew2(ji,jj,jk)*xcond
         denitr(ji,jj,jk)=denitr(ji,jj,jk)*xcond
C
C     Evolution of NH4
C     ----------------
C
         zneg = trn(ji,jj,jk,jpnh4)
     &     -proreg(ji,jj,jk)-proreg2(ji,jj,jk)-onitr(ji,jj,jk)
     &     +grarem(ji,jj,jk)*sigma1+grarem2(ji,jj,jk)*sigma2
     &     +olimi(ji,jj,jk)+denitr(ji,jj,jk)
C
C     Nullity test for NH4
C     --------------------
C
         xcond=(0.5+sign(0.5,zneg))
         prorca(ji,jj,jk)=prorca(ji,jj,jk)*xcond
         prorca2(ji,jj,jk)=prorca2(ji,jj,jk)*xcond
         proreg(ji,jj,jk)=proreg(ji,jj,jk)*xcond
         proreg2(ji,jj,jk)=proreg2(ji,jj,jk)*xcond
         pronew(ji,jj,jk)=pronew(ji,jj,jk)*xcond
         pronew2(ji,jj,jk)=pronew2(ji,jj,jk)*xcond
         onitr(ji,jj,jk)=onitr(ji,jj,jk)*xcond
C
C     Evolution of IRON
C     -----------------
C
          zneg = trn(ji,jj,jk,jpfer)
     &      +(excret-1.)*prorca5(ji,jj,jk)-xaggdfe(ji,jj,jk)
     &      +(excret2-1.)*prorca4(ji,jj,jk)-xbactfer(ji,jj,jk)
     &      +grafer(ji,jj,jk)+grafer2(ji,jj,jk)
     &      +ofer(ji,jj,jk)-xscave(ji,jj,jk)+irondep(ji,jj,jk)
     &      +(ironsed(ji,jj,jk)+po4dep(ji,jj,jk)*9.E-5)*rfact2
C
C     Nullity test for iron
C     ---------------------
C
         xcond=(0.5+sign(0.5,zneg))
         prorca4(ji,jj,jk)=prorca4(ji,jj,jk)*xcond
         prorca5(ji,jj,jk)=prorca5(ji,jj,jk)*xcond
C
C     Evolution of O2
C     ---------------
C
         xcond=(0.5+sign(0.5,(trn(ji,jj,jk,jpoxy)-oxymin)))
         zneg = trn(ji,jj,jk,jpoxy)
     &     +o2ut*(proreg(ji,jj,jk)+proreg2(ji,jj,jk))
     &     +(o2ut+o2nit)*(pronew(ji,jj,jk)+pronew2(ji,jj,jk))
     &     -o2ut*olimi(ji,jj,jk)-o2ut*xcond*(grarem(ji,jj,jk)
     &     *sigma1+grarem2(ji,jj,jk)*sigma2)-o2nit*onitr(ji,jj,jk)
C
C     Nullity test for oxygen
C     -----------------------
C
         xcond=(0.5+sign(0.5,zneg))
         olimi(ji,jj,jk)=olimi(ji,jj,jk)*xcond
         onitr(ji,jj,jk)=onitr(ji,jj,jk)*xcond
C
          END DO
        END DO
      END DO

      DO jk = 1,jpkm1
        DO jj = 1,jpj
          DO ji = 1,jpi
C    
C     Evolution of nanophytoplankton
C     ------------------------------
C     
         zneg = trn(ji,jj,jk,jpphy)
     &     +prorca(ji,jj,jk)*(1.-excret)-tortp(ji,jj,jk)
     &     -grazp(ji,jj,jk)-grazn(ji,jj,jk)-respp(ji,jj,jk)
C     
C     Nullity test for Phyto
C     ----------------------
C     
         xcond=(0.5+sign(0.5,zneg))
         tortp(ji,jj,jk)=tortp(ji,jj,jk)*xcond
         grazp(ji,jj,jk)=grazp(ji,jj,jk)*xcond
         grazn(ji,jj,jk)=grazn(ji,jj,jk)*xcond
         respp(ji,jj,jk)=respp(ji,jj,jk)*xcond
C
C     Evolution of nanophytoplankton chlorophyll
C     ------------------------------
C
         zneg = trn(ji,jj,jk,jpnch)
     &     +prorca6(ji,jj,jk)*(1.-excret)-tortnch(ji,jj,jk)
     &     -grazpch(ji,jj,jk)-graznch(ji,jj,jk)
     &     -respnch(ji,jj,jk)
C
C     Nullity test for Phyto
C     ----------------------
C
         xcond=(0.5+sign(0.5,zneg))
         tortnch(ji,jj,jk)=tortnch(ji,jj,jk)*xcond
         graznch(ji,jj,jk)=graznch(ji,jj,jk)*xcond
         grazpch(ji,jj,jk)=grazpch(ji,jj,jk)*xcond
         respnch(ji,jj,jk)=respnch(ji,jj,jk)*xcond
C
C     Evolution of biogenic Iron in Nanophytoplankton
C     -----------------------------------------------
C
         zneg = trn(ji,jj,jk,jpnfe)
     &     +prorca5(ji,jj,jk)*(1.-excret)-tortnf(ji,jj,jk)
     &     -respnf(ji,jj,jk)-grazpf(ji,jj,jk)-graznf(ji,jj,jk)
C
C     Nullity test for Biogenic Iron in Nanophytoplankton
C     ---------------------------------------------------
C
          xcond=(0.5+sign(0.5,zneg))
          tortnf(ji,jj,jk)=tortnf(ji,jj,jk)*xcond
          respnf(ji,jj,jk)=respnf(ji,jj,jk)*xcond
          grazpf(ji,jj,jk)=grazpf(ji,jj,jk)*xcond
          graznf(ji,jj,jk)=graznf(ji,jj,jk)*xcond
C    
C     Evolution of Diatoms
C     ------------------
C
         zneg = trn(ji,jj,jk,jpdia)
     &     +prorca2(ji,jj,jk)*(1.-excret2)-tortp2(ji,jj,jk)
     &     -respp2(ji,jj,jk)-grazd(ji,jj,jk)-grazsd(ji,jj,jk)
C    
C     Nullity test for diatoms
C     ----------------------
C
         xcond=(0.5+sign(0.5,zneg))
         tortp2(ji,jj,jk)=tortp2(ji,jj,jk)*xcond
         respp2(ji,jj,jk)=respp2(ji,jj,jk)*xcond
         grazd(ji,jj,jk)=grazd(ji,jj,jk)*xcond
         grazsd(ji,jj,jk)=grazsd(ji,jj,jk)*xcond
C   
C     Evolution of Diatoms Chlorophyll
C     ------------------
C
         zneg = trn(ji,jj,jk,jpdch)
     &     +prorca7(ji,jj,jk)*(1.-excret2)-tortdch(ji,jj,jk)
     &     -respdch(ji,jj,jk)-grazdch(ji,jj,jk)-grazsch(ji,jj,jk)
C   
C     Nullity test for diatoms
C     ----------------------
C
         xcond=(0.5+sign(0.5,zneg))
         tortdch(ji,jj,jk)=tortdch(ji,jj,jk)*xcond
         respdch(ji,jj,jk)=respdch(ji,jj,jk)*xcond
         grazdch(ji,jj,jk)=grazdch(ji,jj,jk)*xcond
         grazsch(ji,jj,jk)=grazsch(ji,jj,jk)*xcond
C
C     Evolution of biogenic Iron in diatoms
C     -------------------------------------
C
          zneg = trn(ji,jj,jk,jpdfe)
     &     +prorca4(ji,jj,jk)*(1.-excret2)-grazsf(ji,jj,jk)
     &     -tortdf(ji,jj,jk)-respdf(ji,jj,jk)-grazf(ji,jj,jk)
C
C     Nullity test for Biogenic Iron in diatoms
C     -----------------------------------------
C
          xcond=(0.5+sign(0.5,zneg))
          tortdf(ji,jj,jk)=tortdf(ji,jj,jk)*xcond
          respdf(ji,jj,jk)=respdf(ji,jj,jk)*xcond
          grazf(ji,jj,jk)=grazf(ji,jj,jk)*xcond
          grazsf(ji,jj,jk)=grazsf(ji,jj,jk)*xcond
C
C     Evolution of biogenic Silica in diatoms
C     ---------------------------------------
C
         zneg = trn(ji,jj,jk,jpbsi)
     &     +prorca3(ji,jj,jk)*(1.-excret2)-tortds(ji,jj,jk)
     &     -respds(ji,jj,jk)-grazs(ji,jj,jk)-grazss(ji,jj,jk)
C
C     Nullity test for Biogenic Silica in Diatoms
C     -------------------------------------------
C
          xcond=(0.5+sign(0.5,zneg))
          tortds(ji,jj,jk)=tortds(ji,jj,jk)*xcond
          respds(ji,jj,jk)=respds(ji,jj,jk)*xcond
          grazs(ji,jj,jk)=grazs(ji,jj,jk)*xcond
          grazss(ji,jj,jk)=grazss(ji,jj,jk)*xcond
          END DO
        END DO
      END DO

      DO jk = 1,jpkm1
        DO jj = 1,jpj
          DO ji = 1,jpi
C    
C     Evolution of Zooplankton
C     ------------------------
C    
         zneg = trn(ji,jj,jk,jpzoo)+epsher*
     &     (grazp(ji,jj,jk)+grazm(ji,jj,jk)+grazsd(ji,jj,jk))
     &     -grazz(ji,jj,jk)-tortz(ji,jj,jk)-respz(ji,jj,jk)
C    
C     Nullity test for Zooplankton
C     ----------------------------
C    
         xcond=(0.5+sign(0.5,zneg))
         tortz(ji,jj,jk)=tortz(ji,jj,jk)*xcond
         respz(ji,jj,jk)=respz(ji,jj,jk)*xcond
         grazz(ji,jj,jk)=grazz(ji,jj,jk)*xcond
C
C     Evolution of Mesozooplankton
C     ------------------------
C
         zneg = trn(ji,jj,jk,jpmes)
     &     +epsher2*(grazd(ji,jj,jk)+grazn(ji,jj,jk)+grazz(ji,jj,jk)
     &     +grazpoc(ji,jj,jk)+grazffe(ji,jj,jk))-tortz2(ji,jj,jk)
     &     -respz2(ji,jj,jk)
C
C     Nullity test for Zooplankton
C     ----------------------------
C
         xcond=(0.5+sign(0.5,zneg))
         tortz2(ji,jj,jk)=tortz2(ji,jj,jk)*xcond
         respz2(ji,jj,jk)=respz2(ji,jj,jk)*xcond
          END DO
        END DO
      END DO

      DO jk = 1,jpkm1
        DO jj = 1,jpj
          DO ji = 1,jpi
C     
C     Evolution of detritus
C     ---------------------
C     
         zneg = trn(ji,jj,jk,jppoc)
     &     -grazpoc(ji,jj,jk)+grapoc(ji,jj,jk)-grazm(ji,jj,jk)
     &     +respz(ji,jj,jk)-xagg(ji,jj,jk)+xaggdoc(ji,jj,jk)
     &     +respp(ji,jj,jk)+tortp2(ji,jj,jk)+orem2(ji,jj,jk)
     &     +tortz(ji,jj,jk)+tortp(ji,jj,jk)-orem(ji,jj,jk)
     &     +(sinking(ji,jj,jk)-sinking(ji,jj,jk+1))
     &     /fse3t(ji,jj,jk)
C     
C     Nullity test for POC
C     --------------------
C     
         xcond=(0.5+sign(0.5,zneg))
         grazm(ji,jj,jk)=grazm(ji,jj,jk)*xcond
         sinking(ji,jj,jk+1)=sinking(ji,jj,jk+1)*xcond
         orem(ji,jj,jk)=orem(ji,jj,jk)*xcond
         xagg(ji,jj,jk)=xagg(ji,jj,jk)*xcond
         grazpoc(ji,jj,jk)=grazpoc(ji,jj,jk)*xcond
C   
C     Evolution of detritus
C     ---------------------
C    
         zneg = trn(ji,jj,jk,jpgoc)
     &     +grapoc2(ji,jj,jk)+respp2(ji,jj,jk)+xagg(ji,jj,jk)
     &     +tortz2(ji,jj,jk)+respz2(ji,jj,jk)-orem2(ji,jj,jk)
     &     +xaggdoc2(ji,jj,jk)-grazffe(ji,jj,jk)
     &     +(sinking2(ji,jj,jk)-sinking2(ji,jj,jk+1))
     &     /fse3t(ji,jj,jk)
C
C     Nullity test on goc212
C     ----------------------
C
         xcond=(0.5+sign(0.5,zneg))
         sinking2(ji,jj,jk+1)=sinking2(ji,jj,jk+1)*xcond
         orem2(ji,jj,jk)=orem2(ji,jj,jk)*xcond
C
C     Evolution of small biogenic Iron
C     --------------------------
C
         zdenom=1./(trn(ji,jj,jk,jppoc)+trn(ji,jj,jk,jpgoc)+rtrn)
C
         zneg = trn(ji,jj,jk,jpsfe)
     &     +unass*(grazpf(ji,jj,jk)+grazsf(ji,jj,jk))
     &     -grazpof(ji,jj,jk)-(1.-unass)*grazmf(ji,jj,jk)
     &     +tortdf(ji,jj,jk)+respnf(ji,jj,jk)+tortnf(ji,jj,jk)
     &     +ferat3*(tortz(ji,jj,jk)+respz(ji,jj,jk))-ofer(ji,jj,jk)
     &     +ofer2(ji,jj,jk)-xaggfe(ji,jj,jk)
     &     +xscave(ji,jj,jk)*zdenom1(ji,jj,jk)
     &     +(sinkfer(ji,jj,jk)-sinkfer(ji,jj,jk+1))
     &     /fse3t(ji,jj,jk)
C
C     Nullity test for biogenic iron
C     --------------------
C
         xcond=(0.5+sign(0.5,zneg))
         sinkfer(ji,jj,jk+1)=sinkfer(ji,jj,jk+1)*xcond
         ofer(ji,jj,jk)=ofer(ji,jj,jk)*xcond
         xaggfe(ji,jj,jk)=xaggfe(ji,jj,jk)*xcond
         grazmf(ji,jj,jk)=grazmf(ji,jj,jk)*xcond
C
C     Evolution of big biogenic Iron
C     --------------------------
C
         zneg = trn(ji,jj,jk,jpbfe)
     &     +unass2*(graznf(ji,jj,jk)+grazf(ji,jj,jk)+grazfff(ji,jj,jk)
     &     +grazpof(ji,jj,jk)+ferat3*grazz(ji,jj,jk))+ferat3*
     &     (tortz2(ji,jj,jk)+respz2(ji,jj,jk))-ofer2(ji,jj,jk)
     &     +respdf(ji,jj,jk)+xaggfe(ji,jj,jk)+xbactfer(ji,jj,jk)
     &     -grazfff(ji,jj,jk)+xscave(ji,jj,jk)*zdenom2(ji,jj,jk)
     &     +(sinkfer2(ji,jj,jk)-sinkfer2(ji,jj,jk+1))
     &     /fse3t(ji,jj,jk)
C
C     Nullity test for biogenic iron
C     --------------------
C
         xcond=(0.5+sign(0.5,zneg))
         sinkfer2(ji,jj,jk+1)=sinkfer2(ji,jj,jk+1)*xcond
         ofer2(ji,jj,jk)=ofer2(ji,jj,jk)*xcond
         grazfff(ji,jj,jk)=grazfff(ji,jj,jk)*xcond
C
C     Evolution of sinking biogenic silica
C     --------------------------
C
         zneg = trn(ji,jj,jk,jpdsi)
     &     +tortds(ji,jj,jk)+grazss(ji,jj,jk)
     &     +respds(ji,jj,jk)+grazs(ji,jj,jk)-osil(ji,jj,jk)
     &     +(sinksil(ji,jj,jk)-sinksil(ji,jj,jk+1))
     &     /fse3t(ji,jj,jk)
C
C     Nullity test for Biogenic Silica
C     --------------------------------
C
          xcond=(0.5+sign(0.5,zneg))
          sinksil(ji,jj,jk+1)=sinksil(ji,jj,jk+1)*xcond
          osil(ji,jj,jk)=osil(ji,jj,jk)*xcond
C     
          END DO
        END DO
      END DO
C
C  Recompute the SMS related to zooplankton grazing
C  ------------------------------------------------
C
      DO jk = 1,jpkm1
        DO jj = 1,jpj
          DO ji = 1,jpi
         grarem(ji,jj,jk)=(grazp(ji,jj,jk)+grazm(ji,jj,jk)
     &      +grazsd(ji,jj,jk))*(1.-epsher-unass)

        grafer(ji,jj,jk)=(grazpf(ji,jj,jk)+grazsf(ji,jj,jk)
     &      +grazmf(ji,jj,jk))*(1.-epsher-unass)
     &      +(grazm(ji,jj,jk)*max((trn(ji,jj,jk,jpsfe)/
     &      (trn(ji,jj,jk,jppoc)+rtrn)-ferat3),0.)
     &      +grazp(ji,jj,jk)*max((trn(ji,jj,jk,jpnfe)/
     &      (trn(ji,jj,jk,jpphy)+rtrn)-ferat3),0.)
     &      +grazsd(ji,jj,jk)*max((trn(ji,jj,jk,jpdfe)/
     &      (trn(ji,jj,jk,jpdia)+rtrn)-ferat3),0.))*epsher

        grarem2(ji,jj,jk)=(grazd(ji,jj,jk)+grazz(ji,jj,jk)
     &      +grazn(ji,jj,jk)+grazpoc(ji,jj,jk)+grazffe(ji,jj,jk))
     &      *(1.-epsher2-unass2)

        grafer2(ji,jj,jk)=(grazf(ji,jj,jk)+graznf(ji,jj,jk)
     &    +grazz(ji,jj,jk)*ferat3+grazpof(ji,jj,jk)
     &    +grazfff(ji,jj,jk))*(1.-epsher2-unass2)
     &    +epsher2*(grazd(ji,jj,jk)*max(
     &    (trn(ji,jj,jk,jpdfe)/(trn(ji,jj,jk,jpdia)+rtrn)
     &    -ferat3),0.)+grazn(ji,jj,jk)*max(
     &    (trn(ji,jj,jk,jpnfe)/(trn(ji,jj,jk,jpphy)+rtrn)
     &    -ferat3),0.)+grazpoc(ji,jj,jk)*max(
     &    (trn(ji,jj,jk,jpsfe)/(trn(ji,jj,jk,jppoc)+rtrn)
     &    -ferat3),0.)+grazffe(ji,jj,jk)*max(
     &    (trn(ji,jj,jk,jpbfe)/(trn(ji,jj,jk,jpgoc)+rtrn)
     &    -ferat3),0.))

        grapoc2(ji,jj,jk)=(grazd(ji,jj,jk)+grazz(ji,jj,jk)
     &    +grazn(ji,jj,jk)+grazpoc(ji,jj,jk)+grazffe(ji,jj,jk))*unass2

        grapoc(ji,jj,jk)=(grazp(ji,jj,jk)+grazm(ji,jj,jk)
     &      +grazsd(ji,jj,jk))*unass
          END DO
        END DO
      END DO
C     
C     Determination of tracers concentration as a function of 
C     biological sources and sinks
C     --------------------------------------------------------
C     
      DO jk = 1,jpkm1
        DO jj = 1,jpj
          DO ji = 1,jpi
C     
C     Evolution of PO4
C     ----------------
C     
          trn(ji,jj,jk,jppo4) = trn(ji,jj,jk,jppo4)
     &      -prorca(ji,jj,jk)-prorca2(ji,jj,jk)
     &      +olimi(ji,jj,jk)+grarem(ji,jj,jk)*sigma1+denitr(ji,jj,jk)
     &      +grarem2(ji,jj,jk)*sigma2+po4dep(ji,jj,jk)*rfact2
C
C     Evolution of NO3 and NH4
C     ------------------------
C
          trn(ji,jj,jk,jpno3) = trn(ji,jj,jk,jpno3)
     &      -pronew(ji,jj,jk)-pronew2(ji,jj,jk)+onitr(ji,jj,jk)
     &      -denitr(ji,jj,jk)*rdenit+po4dep(ji,jj,jk)*rfact2
     &      +nitdep(ji,jj,jk)*rfact2

          trn(ji,jj,jk,jpnh4) = trn(ji,jj,jk,jpnh4)
     &      -proreg(ji,jj,jk)-proreg2(ji,jj,jk)+olimi(ji,jj,jk)
     &      +grarem(ji,jj,jk)*sigma1+grarem2(ji,jj,jk)*sigma2
     &      -onitr(ji,jj,jk)+denitr(ji,jj,jk)

          END DO
        END DO
      END DO

      DO jk = 1,jpkm1
        DO jj = 1,jpj
          DO ji = 1,jpi

C    
C     Evolution of Phytoplankton
C     --------------------------
C     
          trn(ji,jj,jk,jpphy) = trn(ji,jj,jk,jpphy)
     &      +prorca(ji,jj,jk)*(1.-excret)-tortp(ji,jj,jk)
     &      -grazp(ji,jj,jk)-grazn(ji,jj,jk)-respp(ji,jj,jk)

          trn(ji,jj,jk,jpnch) = trn(ji,jj,jk,jpnch)
     &      +prorca6(ji,jj,jk)*(1.-excret)-tortnch(ji,jj,jk)
     &      -grazpch(ji,jj,jk)-graznch(ji,jj,jk)-respnch(ji,jj,jk)
C
C     Evolution of Diatoms
C     ------------------
C
          trn(ji,jj,jk,jpdia) = trn(ji,jj,jk,jpdia)
     &      +prorca2(ji,jj,jk)*(1.-excret2)-tortp2(ji,jj,jk)
     &      -respp2(ji,jj,jk)-grazd(ji,jj,jk)-grazsd(ji,jj,jk)

          trn(ji,jj,jk,jpdch) = trn(ji,jj,jk,jpdch)
     &      +prorca7(ji,jj,jk)*(1.-excret2)-tortdch(ji,jj,jk)
     &      -respdch(ji,jj,jk)-grazdch(ji,jj,jk)-grazsch(ji,jj,jk)
          END DO
        END DO
      END DO

      DO jk = 1,jpkm1
        DO jj = 1,jpj
          DO ji = 1,jpi
C    
C     Evolution of Zooplankton
C     ------------------------
C     
          trn(ji,jj,jk,jpzoo) = trn(ji,jj,jk,jpzoo)
     &      +epsher*(grazp(ji,jj,jk)+grazm(ji,jj,jk)+grazsd(ji,jj,jk))
     &      -grazz(ji,jj,jk)-tortz(ji,jj,jk)-respz(ji,jj,jk)
C
C     Evolution of Mesozooplankton
C     ------------------------
C
          trn(ji,jj,jk,jpmes) = trn(ji,jj,jk,jpmes)
     &      +epsher2*(grazd(ji,jj,jk)+grazz(ji,jj,jk)+grazn(ji,jj,jk)
     &      +grazpoc(ji,jj,jk)+grazffe(ji,jj,jk))-tortz2(ji,jj,jk)
     &      -respz2(ji,jj,jk)
          END DO
        END DO
      END DO

      DO jk = 1,jpkm1
        DO jj = 1,jpj
          DO ji = 1,jpi
C    
C     Evolution of DOC
C     ----------------
C     
          trn(ji,jj,jk,jpdoc) = trn(ji,jj,jk,jpdoc)
     &      +orem(ji,jj,jk)+excret2*prorca2(ji,jj,jk)
     &      +excret*prorca(ji,jj,jk)-olimi(ji,jj,jk)-denitr(ji,jj,jk)
     &      +grarem(ji,jj,jk)*(1.-sigma1)+grarem2(ji,jj,jk)
     &      *(1.-sigma2)-xaggdoc(ji,jj,jk)-xaggdoc2(ji,jj,jk)
          END DO
        END DO
      END DO

      DO jk = 1,jpkm1
        DO jj = 1,jpj
          DO ji = 1,jpi
C     
C     Evolution of Detritus
C     ---------------------
C     
          trn(ji,jj,jk,jppoc) = trn(ji,jj,jk,jppoc)
     &      -grazpoc(ji,jj,jk)+grapoc(ji,jj,jk)+tortp2(ji,jj,jk)
     &      -grazm(ji,jj,jk)+respp(ji,jj,jk)+tortz(ji,jj,jk)
     &      +tortp(ji,jj,jk)+respz(ji,jj,jk)-orem(ji,jj,jk)
     &      +orem2(ji,jj,jk)-xagg(ji,jj,jk)+xaggdoc(ji,jj,jk)
     &      +(sinking(ji,jj,jk)-sinking(ji,jj,jk+1))
     &      /fse3t(ji,jj,jk)
C    
C     Evolution of rapid Detritus
C     ---------------------
C    
          trn(ji,jj,jk,jpgoc) = trn(ji,jj,jk,jpgoc)
     &      +grapoc2(ji,jj,jk)+respp2(ji,jj,jk)+xagg(ji,jj,jk)
     &      +tortz2(ji,jj,jk)+respz2(ji,jj,jk)-orem2(ji,jj,jk)
     &      -grazffe(ji,jj,jk)+xaggdoc2(ji,jj,jk)
     &      +(sinking2(ji,jj,jk)-sinking2(ji,jj,jk+1))
     &      /fse3t(ji,jj,jk)
          END DO
        END DO
      END DO

      DO jk = 1,jpkm1
        DO jj = 1,jpj
          DO ji = 1,jpi
C    
C     Evolution of O2
C     ---------------
C     
         xcond=(0.5+sign(0.5,(trn(ji,jj,jk,jpoxy)-oxymin)))
         trn(ji,jj,jk,jpoxy) = trn(ji,jj,jk,jpoxy)
     &     +o2ut*(proreg(ji,jj,jk)+proreg2(ji,jj,jk)-olimi(ji,jj,jk)
     &     -xcond*(grarem(ji,jj,jk)*sigma1+grarem2(ji,jj,jk)*sigma2))
     &     +(o2ut+o2nit)*( pronew(ji,jj,jk)+pronew2(ji,jj,jk))
     &     -o2nit*onitr(ji,jj,jk)
C
          END DO
        END DO
      END DO

      DO jk = 1,jpkm1
        DO jj = 1,jpj
          DO ji = 1,jpi
C
C     Evolution of IRON
C     -----------------
C
          trn(ji,jj,jk,jpfer) = trn(ji,jj,jk,jpfer)
     &      +(excret-1.)*prorca5(ji,jj,jk)-xaggdfe(ji,jj,jk)
     &      +(excret2-1.)*prorca4(ji,jj,jk)-xbactfer(ji,jj,jk)
     &      +grafer(ji,jj,jk)+grafer2(ji,jj,jk)
     &      +ofer(ji,jj,jk)-xscave(ji,jj,jk)+irondep(ji,jj,jk)
     &      +(ironsed(ji,jj,jk)+po4dep(ji,jj,jk)*9E-5)*rfact2
          END DO
        END DO
      END DO

      DO jk = 1,jpkm1
        DO jj = 1,jpj
          DO ji = 1,jpi
C
C     Evolution of small biogenic Iron
C     --------------------------
C
          zdenom=1./(trn(ji,jj,jk,jppoc)+trn(ji,jj,jk,jpgoc)+rtrn)
C
          trn(ji,jj,jk,jpsfe) = trn(ji,jj,jk,jpsfe)
     &     +unass*(grazpf(ji,jj,jk)+grazsf(ji,jj,jk))
     &     -grazpof(ji,jj,jk)-(1.-unass)*grazmf(ji,jj,jk)
     &     +tortdf(ji,jj,jk)+respnf(ji,jj,jk)+tortnf(ji,jj,jk)
     &     +ferat3*(tortz(ji,jj,jk)+respz(ji,jj,jk))-ofer(ji,jj,jk)
     &     +ofer2(ji,jj,jk)-xaggfe(ji,jj,jk)
     &     +xscave(ji,jj,jk)*zdenom1(ji,jj,jk)
     &     +(sinkfer(ji,jj,jk)-sinkfer(ji,jj,jk+1))
     &     /fse3t(ji,jj,jk)
C
C     Evolution of big biogenic Iron
C     --------------------------
C
          trn(ji,jj,jk,jpbfe) = trn(ji,jj,jk,jpbfe)
     &     +unass2*(graznf(ji,jj,jk)+grazf(ji,jj,jk)+grazfff(ji,jj,jk)
     &     +grazpof(ji,jj,jk)+grazz(ji,jj,jk)*ferat3)+ferat3*
     &     (tortz2(ji,jj,jk)+respz2(ji,jj,jk))-ofer2(ji,jj,jk)
     &     +respdf(ji,jj,jk)+xaggfe(ji,jj,jk)+xbactfer(ji,jj,jk)
     &     -grazfff(ji,jj,jk)+xscave(ji,jj,jk)*zdenom2(ji,jj,jk)
     &     +(sinkfer2(ji,jj,jk)-sinkfer2(ji,jj,jk+1))
     &     /fse3t(ji,jj,jk)
          END DO
        END DO
      END DO

      DO jk = 1,jpkm1
        DO jj = 1,jpj
          DO ji = 1,jpi
C
C     Evolution of biogenic Silica
C     ----------------------------
C
          trn(ji,jj,jk,jpbsi) = trn(ji,jj,jk,jpbsi)
     &      +prorca3(ji,jj,jk)*(1.-excret2)-grazss(ji,jj,jk)
     &      -tortds(ji,jj,jk)-respds(ji,jj,jk)-grazs(ji,jj,jk)
C
          silpro(ji,jj,jk)=
     &      tortds(ji,jj,jk)+respds(ji,jj,jk)+grazs(ji,jj,jk)
     &      +grazss(ji,jj,jk)
C
          END DO
        END DO
      END DO

      DO jk = 1,jpkm1
        DO jj = 1,jpj
          DO ji = 1,jpi
C
C     Evolution of sinking biogenic silica
C     ------------------------------------
C
          trn(ji,jj,jk,jpdsi)=trn(ji,jj,jk,jpdsi)
     &      +tortds(ji,jj,jk)+respds(ji,jj,jk)+grazs(ji,jj,jk)
     &      -osil(ji,jj,jk)+grazss(ji,jj,jk)
     &      +(sinksil(ji,jj,jk)-sinksil(ji,jj,jk+1))
     &      /fse3t(ji,jj,jk)
C
          END DO
        END DO
      END DO

      DO jk = 1,jpkm1
        DO jj = 1,jpj
          DO ji = 1,jpi
C
C     Evolution of biogenic diatom Iron
C     ---------------------------------
C
          trn(ji,jj,jk,jpdfe) = trn(ji,jj,jk,jpdfe)
     &      +prorca4(ji,jj,jk)*(1.-excret2)-grazsf(ji,jj,jk)
     &      -tortdf(ji,jj,jk)-respdf(ji,jj,jk)-grazf(ji,jj,jk)
C
C     Evolution of biogenic nanophytoplankton Iron
C     --------------------------------------------
C
          trn(ji,jj,jk,jpnfe) = trn(ji,jj,jk,jpnfe)
     &      +prorca5(ji,jj,jk)*(1.-excret)-graznf(ji,jj,jk)
     &      -tortnf(ji,jj,jk)-respnf(ji,jj,jk)-grazpf(ji,jj,jk)
C
C     Evolution of dissolved Silica
C     -----------------------------
C
          trn(ji,jj,jk,jpsil) = trn(ji,jj,jk,jpsil)
     &      -(1.-excret2)*prorca3(ji,jj,jk)+osil(ji,jj,jk)
     &      +sidep(ji,jj,jk)+cotdep(ji,jj,jk)*rfact2/6.
C
          END DO
        END DO
      END DO
C     
C     Evolution of calcite and silicates as a function of the two tracers
C     -------------------------------------------------------------------
C     
      DO  jk = 1,jpkm1
        DO  jj = 1,jpj
          DO  ji = 1,jpi
C
          prodt = prorca(ji,jj,jk)+prorca2(ji,jj,jk)
     &      -olimi(ji,jj,jk)-grarem(ji,jj,jk)*sigma1
     &      -grarem2(ji,jj,jk)*sigma2-denitr(ji,jj,jk)

          prodca = pronew(ji,jj,jk)+pronew2(ji,jj,jk)
     &      -onitr(ji,jj,jk)+rdenit*denitr(ji,jj,jk)
     &      -po4dep(ji,jj,jk)*rfact2-nitdep(ji,jj,jk)*rfact2
C     
C     potential production of calcite and biogenic silicate
C     ------------------------------------------------------
C     
          prcaca(ji,jj,jk)=caco3r*(0.5*(unass*grazp(ji,jj,jk)+
     &      unass2*grazn(ji,jj,jk))+tortp(ji,jj,jk)+respp(ji,jj,jk))
     &      *xlimphy(ji,jj,jk)*xlimphy(ji,jj,jk)
C     
C     Consumption of Total (12C)O2
C     ----------------------------
C     
          trn(ji,jj,jk,jpdic) = trn(ji,jj,jk,jpdic)
     &      -prodt-prcaca(ji,jj,jk)+po4dep(ji,jj,jk)*rfact2*2.633
C     
C     Consumption of alkalinity due to ca++ uptake and increase 
C     of alkalinity due to nitrate consumption during organic 
C     soft tissue production
C     ---------------------------------------------------------
C     
          trn(ji,jj,jk,jptal) = trn(ji,jj,jk,jptal)
     &      +rno3*prodca-2.*prcaca(ji,jj,jk)
     &      +cotdep(ji,jj,jk)*rfact2
          END DO
        END DO
      END DO
C
      DO  jk = 1,jpkm1
        DO  jj = 1,jpj
          DO  ji = 1,jpi
C
C     Production of calcite due to biological production
C     --------------------------------------------------
C     
           trn(ji,jj,jk,jpcal) = trn(ji,jj,jk,jpcal)
     &        +prcaca(ji,jj,jk)+(sinkcal(ji,jj,jk)-
     &         sinkcal(ji,jj,jk+1))/fse3t(ji,jj,jk)
          END DO
        END DO
      ENDDO
C
      DO jn=1 , jptra
        CALL lbc_lnk(trn(:,:,:,jn), 'T', 1. )
      END DO

#    if defined key_trc_diaadd
       DO jj=1,jpj
         DO ji=1,jpi
          trc2d(ji,jj,12) = irondep(ji,jj,1)*1e3*rfact2r
     &       *fse3t(ji,jj,1)
         END DO
       END DO
#    endif
C
#    if defined key_trc_dia3d
          trc3d(:,:,:,4)=etot(:,:,:)
          trc3d(:,:,:,5)=prorca(:,:,:)*1e3*rfact2r
          trc3d(:,:,:,6)=prorca2(:,:,:)*1e3*rfact2r
          trc3d(:,:,:,7)=pronew(:,:,:)*1e3*rfact2r
          trc3d(:,:,:,8)=pronew2(:,:,:)*1e3*rfact2r
          trc3d(:,:,:,9)=prorca3(:,:,:)*1e3*rfact2r
          trc3d(:,:,:,10)=prorca4(:,:,:)*1e3*rfact2r
          trc3d(:,:,:,11)=prorca5(:,:,:)*1e3*rfact2r
#    endif
C     
#endif
C    
 
      RETURN
      END