CCC$Header$
CCC  TOP 1.0 , LOCEAN-IPSL (2005)
C This software is governed by CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
C ---------------------------------------------------------------------------
CCC
CCC      trcbio.lobster1.h
CCC      *****************
CCC
CC   defined key : key_trc_lobster1
CC   ===========
CC
CC   INPUT :
CC   -----
CC      argument
CC              ktask           : task identificator
CC              kt              : time step
CC      COMMON
CC            /comcoo/          : orthogonal curvilinear coordinates
CC                                and scale factors
CC                                depths
CC            /cottrp/          : present and next fields for passive
CC                              : tracers
CC            /comtsk/          : multitasking
CC            /comtke/          : emin, en()
CC            /cotbio/          : biological parameters
CC
CC   OUTPUT :
CC   ------
CC      COMMON
CC            /cottrp/ tra      : general tracer trend increased by the
CC                                now horizontal tracer advection trend
CC            /cottbd/ trbio    : now horizontal tracer advection trend
CC                                (IF 'key_trc_diabio' is activated)
CC
CC   WORKSPACE :
CC   ---------
CC      local
CC               zdet,zzoo,zphy,znh4,zno3,zdom    : now concentrations
CC               zlt,zlno3,zlnh4,zle              : limitation terms for phyto
CC               zfno3phy and so on..             : fluxes between bio boxes
CC               zphya,zzooa,zdeta, ...           : after bio trends
CC               zppz, zpdz, zpppz, zppdz, zfood  : preferences terms
CC               zfilpz, zfilpd                   : filtration terms
CC      COMMON
CC
CC   EXTERNAL :                   no
CC   --------
CC
CC   REFERENCES :                 no
CC   ----------
CC
CC   MODIFICATIONS:
CC   --------------
CC       original : 99-07 (M. Levy)
CC                  00-12 (E. Kestenare): assign a parameter 
CC                                        to name individual tracers
CC                  01-03 (M. Levy) LNO3 + dia2d
CC----------------------------------------------------------------------
CC----------------------------------------------------------------------
      USE oce_trc
      USE trp_trc
      USE sms
      USE lbclnk
      IMPLICIT NONE
CC local declarations
CC ==================
      INTEGER kt
      INTEGER ji,jj,jk,jn
      REAL ztot(jpi), ze3t(jpk)
      REAL zdet,zzoo,zphy,zno3,znh4,zdom,zlno3,zlnh4,zle,zlt
      REAL zno3phy, znh4phy, zphynh4, zphydom, zphydet, zphyzoo, zdetzoo
     $    ,zzoonh4, zzoodom, zzoodet, zdetnh4, zdetdom, znh4no3, zdomnh4
     $    ,zppz,zpdz,zpppz,zppdz,zfood,zfilpz,zfildz,zphya,zzooa,zno3a
     $    ,znh4a,zdeta,zdoma, ztra, zzoobod

CC----------------------------------------------------------------------
CC statement functions
CC ===================
CDIR$ NOLIST
#include "domzgr_substitute.h90"
CDIR$ LIST
CCC---------------------------------------------------------------------
CCC  OPA8, LODYC (07/99)
CCC---------------------------------------------------------------------
C   | --------------|
C   | LOBSTER1 MODEL| 
C   | --------------|

#if defined key_trc_diaadd
C convert fluxes in per day
      DO jk=1,jpkbm1
        ze3t(jk)=e3t(jk)*86400.
      END DO 
      DO jk=jpkb,jpk
        ze3t(jk)=0.
      END DO 
#endif
C
C vertical slab
C =============
C
      DO 1000 jj = 2,jpjm1
C
C 1. biological level
C ===================
C
        DO ji = 2,jpim1
          fbod(ji,jj)=0.
#if defined key_trc_diaadd
          DO jn=1,jpdia2d
            trc2d(ji,jj,jn)=0.          
          END DO 
#endif
        END DO 

        DO jk=1,jpkbm1
          DO ji = 2,jpim1
C
C
C 1.1 trophic variables( det, zoo, phy, no3, nh4, dom)
C ---------------------------------------------------
C
C negative trophic variables DO not contribute to the fluxes
C
            zdet = max(0.,trn(ji,jj,jk,jpdet))
            zzoo = max(0.,trn(ji,jj,jk,jpzoo))
            zphy = max(0.,trn(ji,jj,jk,jpphy))
            zno3 = max(0.,trn(ji,jj,jk,jpno3))
            znh4 = max(0.,trn(ji,jj,jk,jpnh4))
            zdom = max(0.,trn(ji,jj,jk,jpdom))
C
C
C 1.2  Limitations
C ----------------
C
            zlt = 1.
            zle = 1. - exp( -xpar(ji,jj,jk)/aki/zlt) 
c            zlno3 = (aknut * zno3)/ (( aknut+znh4)*(aknut+zno3+znh4))
            zlno3 = zno3* exp(-1.5*znh4) / (aknut+zno3)
            zlnh4 = znh4 / (znh4+aknut) 

C
C
C 1.3 sinks and sources
C ---------------------
C
C
C 1. phytoplankton production and exsudation
C
            zno3phy = tmumax * zle * zlt * zlno3 * zphy
            znh4phy = tmumax * zle * zlt * zlnh4 * zphy

            zphydom = rgamma * (1-fdoml) * (zno3phy + znh4phy)
            zphynh4 = rgamma * fdoml     * (zno3phy + znh4phy)

C
C 2. zooplankton production
C
C preferences
C
            zppz = rppz
            zpdz = 1. - rppz
            zpppz = ( zppz * zphy ) /
     $          ( ( zppz * zphy + zpdz * zdet ) + 1.e-13 )
            zppdz = ( zpdz * zdet ) /
     $          ( ( zppz * zphy + zpdz * zdet ) + 1.e-13 )
            zfood = zpppz * zphy + zppdz * zdet
C
C filtration
C
            zfilpz = taus * zpppz / (aks + zfood)
            zfildz = taus * zppdz / (aks + zfood)
C
C grazing
C
            zphyzoo = zfilpz * zphy * zzoo
            zdetzoo = zfildz * zdet * zzoo
C
C 3. fecal pellets production
C
            zzoodet = rpnaz * zphyzoo + rdnaz * zdetzoo
C
C 4. zooplankton liquide excretion
C
            zzoonh4 = tauzn * zzoo * fdoml 
            zzoodom = tauzn * zzoo * (1-fdoml) 

C 5. mortality
C
C phytoplankton mortality 
C
            zphydet = tmminp * zphy
C
C
C zooplankton mortality
c closure : flux fbod is redistributed below level jpkbio
C
            zzoobod = tmminz * zzoo * zzoo
            fbod(ji,jj) = fbod(ji,jj) + zzoobod * fse3t(ji,jj,jk)
C
C
C 6. detritus and dom breakdown
C
C
            zdetnh4 = taudn * fdoml * zdet
            zdetdom = taudn * (1 - fdoml) * zdet
            zdomnh4 = taudomn * zdom
C
C
C 7. Nitrification
C
            znh4no3 = taunn * znh4
C
C
C
C 1.4 determination of trends
C ---------------------------
C
C total trend for each biological tracer
C
            zphya =   zno3phy + znh4phy - zphynh4 - zphydom - zphyzoo
     $          - zphydet
            zzooa =   zphyzoo + zdetzoo - zzoodet - zzoodom - zzoonh4
     $          - zzoobod
            zno3a = - zno3phy + znh4no3
            znh4a = - znh4phy - znh4no3 + zphynh4 + zzoonh4 + zdomnh4
     $          + zdetnh4
            zdeta = zphydet + zzoodet  - zdetzoo - zdetnh4 - zdetdom
            zdoma = zphydom + zzoodom + zdetdom - zdomnh4
C
#if defined key_trc_diabio
            trbio(ji,jj,jk,1) = zno3phy
            trbio(ji,jj,jk,2) = znh4phy
            trbio(ji,jj,jk,3) = zphynh4
            trbio(ji,jj,jk,4) = zphydom
            trbio(ji,jj,jk,5) = zphyzoo
            trbio(ji,jj,jk,6) = zphydet
            trbio(ji,jj,jk,7) = zdetzoo
            trbio(ji,jj,jk,9) = zzoodet
            trbio(ji,jj,jk,10) = zzoobod
            trbio(ji,jj,jk,11) = zzoonh4
            trbio(ji,jj,jk,12) = zzoodom
            trbio(ji,jj,jk,13) = znh4no3
            trbio(ji,jj,jk,14) = zdomnh4
            trbio(ji,jj,jk,15) = zdetnh4
#endif
#if defined key_trc_diaadd
            trc2d(ji,jj,1)=trc2d(ji,jj,1)+zno3phy*ze3t(jk)          
            trc2d(ji,jj,2)=trc2d(ji,jj,2)+znh4phy*ze3t(jk)
            trc2d(ji,jj,3)=trc2d(ji,jj,3)+zphydom*ze3t(jk)
            trc2d(ji,jj,4)=trc2d(ji,jj,4)+zphynh4*ze3t(jk)
            trc2d(ji,jj,5)=trc2d(ji,jj,5)+zphyzoo*ze3t(jk)
            trc2d(ji,jj,6)=trc2d(ji,jj,6)+zphydet*ze3t(jk)
            trc2d(ji,jj,7)=trc2d(ji,jj,7)+zdetzoo*ze3t(jk)
c trend number 8 is in trcsed.F            
            trc2d(ji,jj,9)=trc2d(ji,jj,9)+zzoodet*ze3t(jk)
            trc2d(ji,jj,10)=trc2d(ji,jj,10)+zzoobod*ze3t(jk)
            trc2d(ji,jj,11)=trc2d(ji,jj,11)+zzoonh4*ze3t(jk)
            trc2d(ji,jj,12)=trc2d(ji,jj,12)+zzoodom*ze3t(jk)
            trc2d(ji,jj,13)=trc2d(ji,jj,13)+znh4no3*ze3t(jk)
            trc2d(ji,jj,14)=trc2d(ji,jj,14)+zdomnh4*ze3t(jk)
            trc2d(ji,jj,15)=trc2d(ji,jj,15)+zdetnh4*ze3t(jk)
             
            trc2d(ji,jj,16)=trc2d(ji,jj,16)+(zno3phy+znh4phy-zphynh4
     $          -zphydom-zphyzoo-zphydet)*ze3t(jk)
            trc2d(ji,jj,17)=trc2d(ji,jj,17)+(zphyzoo+zdetzoo-zzoodet
     $          -zzoobod-zzoonh4-zzoodom) *ze3t(jk)
            trc2d(ji,jj,18)=trc2d(ji,jj,18)+zdetdom*ze3t(jk)

            trc3d(ji,jj,jk,1)= zno3phy *86400     
            trc3d(ji,jj,jk,2)= znh4phy *86400     
            trc3d(ji,jj,jk,3)= znh4no3 *86400     
#endif
C
C tracer flux at totox-point added to the general trend
C
            tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) + zdeta
            tra(ji,jj,jk,jpzoo) = tra(ji,jj,jk,jpzoo) + zzooa
            tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) + zphya
            tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) + zno3a
            tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + znh4a
            tra(ji,jj,jk,jpdom) = tra(ji,jj,jk,jpdom) + zdoma
C
          END DO
        END DO
C
C 2. under biological level
C =========================
C
        DO jk = jpkb,jpk
C
C 2.1 compute the remineralisation of all quantities towards nitrate 
C ------------------------------------------------------------------
C
          DO ji=2,jpim1
            ztot(ji) = 0.
          END DO 
          DO jn=1,jptra
            IF (ctrcnm(jn).NE.'NO3') THEN 
                DO ji=2,jpim1
                  ztra = remdmp(jk,jn) * trn(ji,jj,jk,jn) 
                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) - ztra
                  ztot(ji) = ztot(ji) + ztra
                END DO 
            ENDIF
          END DO 
          DO jn=1,jptra
            IF (ctrcnm(jn).EQ.'NO3') THEN 
                DO ji=2,jpim1
                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztot(ji)
                END DO
#if defined key_trc_diabio
                trbio(ji,jj,jk,1)=ztot(ji)
#endif 
            ENDIF
          END DO
        END DO 


#if defined key_trc_diabio
        DO jk=jpkb,jpkm1
          DO ji=2,jpim1
            trbio(ji,jj,jk,2)=tra(ji,jj,jk,jptra)
          END DO
        END DO 
#endif
C
C
C END of slab
C ===========
C
 1000 CONTINUE

#if defined key_trc_diaadd

C Lateral boundary conditions on trc2d
      DO jn=1,jpdia2d
          CALL lbc_lnk(trc2d(:,:,jn),'T',1. )
      END DO 

C Lateral boundary conditions on trc3d
      DO jn=1,jpdia3d
          CALL lbc_lnk(trc3d(:,:,1,jn),'T',1. )
      END DO 

#endif

#if defined key_trc_diabio
C Lateral boundary conditions on trcbio
      DO jn=1,jpdiabio
          CALL lbc_lnk(trbio(:,:,1,jn),'T',1. )
      END DO 
#endif