source: trunk/NEMO/TOP_SRC/SMS/p4zsed.F @ 247

CCC $Header$ 
CCC  TOP 1.0 , LOCEAN-IPSL (2005) 
C This software is governed by CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
C ---------------------------------------------------------------------------
CDIR$ LIST
      SUBROUTINE p4zsed
#if defined key_passivetrc && defined key_trc_pisces
CCC---------------------------------------------------------------------
CCC
CCC          ROUTINE p4zsed : PISCES MODEL
CCC          *****************************
CCC
CCC  PURPOSE :
CCC  ---------
CCC         Compute loss of organic matter in the sediments. This
CCC         is by no way a sediment model. The loss is simply 
CCC         computed to balance the inout from rivers and dust
CCC
CC   INPUT :
CC   -----
CC      common
CC              all the common defined in opa
CC
CC
CC   OUTPUT :                   : no
CC   ------
CC
CC   EXTERNAL :
CC   --------
CC             None
CC
CC   MODIFICATIONS:
CC   --------------
CC      original  : 2004 - O. Aumont 
CC----------------------------------------------------------------------
CC parameters and commons
CC ======================
      USE oce_trc
      USE trp_trc
      USE sms
      IMPLICIT NONE
CC----------------------------------------------------------------------
CC local declarations
CC ==================
      INTEGER ji, jj, jk, ikt
      REAL sumsedsi,sumsedpo4,sumsedcal
      REAL xconctmp,denitot,nitrpottot,nitrpot(jpi,jpj)

CC
CC----------------------------------------------------------------------
CC statement functions
CC ===================
CDIR$ NOLIST
#include "domzgr_substitute.h90"
CDIR$ LIST
C
C
C     Initialisation of variables used to compute Sinking Speed
C     ---------------------------------------------------------
C
        sumsedsi = 0.
        sumsedpo4 = 0.
        sumsedcal = 0.
C
C    Loss of biogenic silicon, Caco3 organic carbon in the sediments. 
C    First, the total loss is computed.
C    The factor for calcite comes from the alkalinity effect
C    -------------------------------------------------------------
C
        DO jj=2,jpjm1
          DO ji=2,jpim1
        ikt=max(mbathy(ji,jj)-1,1)
        sumsedsi=sumsedsi+trn(ji,jj,ikt,jpdsi)*e1t(ji,jj)
     &    *e2t(ji,jj)*wsbio4(ji,jj,ikt)/rjjss
        sumsedcal=sumsedcal+trn(ji,jj,ikt,jpcal)*e1t(ji,jj)
     &    *e2t(ji,jj)*wsbio4(ji,jj,ikt)*2./rjjss
        sumsedpo4=sumsedpo4+(trn(ji,jj,ikt,jpgoc)*wsbio4(ji,jj,ikt)
     &    +trn(ji,jj,ikt,jppoc)*wsbio3(ji,jj,ikt))/rjjss
     &    *e1t(ji,jj)*e2t(ji,jj)
          END DO
        END DO

C
C    Then this loss is scaled at each bottom grid cell for
C    equilibrating the total budget of silica in the ocean.
C    Thus, the amount of silica lost in the sediments equal
C    the supply at the surface (dust+rivers)
C    ------------------------------------------------------
C
        DO jj=1,jpj
          DO ji=1,jpi
        ikt=max(mbathy(ji,jj)-1,1)
        xconctmp=trn(ji,jj,ikt,jpdsi)
        trn(ji,jj,ikt,jpdsi)=trn(ji,jj,ikt,jpdsi)
     &    -xconctmp*wsbio4(ji,jj,ikt)
     &    *rfact2/rjjss/fse3t(ji,jj,ikt)
        trn(ji,jj,ikt,jpsil)=trn(ji,jj,ikt,jpsil)
     &    +xconctmp*wsbio4(ji,jj,ikt)
     &    *rfact2/rjjss/fse3t(ji,jj,ikt)*(1.-(sumdepsi+rivalkinput
     &    /raass/6.)/sumsedsi)
          END DO
        END DO

        DO jj=1,jpj
          DO ji=1,jpi
        ikt=max(mbathy(ji,jj),1)
        xconctmp=trn(ji,jj,ikt,jpcal)
        trn(ji,jj,ikt,jpcal)=trn(ji,jj,ikt,jpcal)
     &    -xconctmp*wsbio4(ji,jj,ikt)
     &    *rfact2/rjjss/fse3t(ji,jj,ikt)
        trn(ji,jj,ikt,jptal)=trn(ji,jj,ikt,jptal)
     &    +xconctmp*wsbio4(ji,jj,ikt)
     &    *rfact2/rjjss/fse3t(ji,jj,ikt)*(1.-(rivalkinput
     &    /raass)/sumsedcal)*2.
        trn(ji,jj,ikt,jpdic)=trn(ji,jj,ikt,jpdic)
     &    +xconctmp*wsbio4(ji,jj,ikt)
     &    *rfact2/rjjss/fse3t(ji,jj,ikt)*(1.-(rivalkinput
     &    /raass)/sumsedcal)
         END DO
       END DO

        DO jj=1,jpj
          DO ji=1,jpi
        ikt=max(mbathy(ji,jj),1)
        trn(ji,jj,ikt,jpgoc)=trn(ji,jj,ikt,jpgoc)
     &    -trn(ji,jj,ikt,jpgoc)*wsbio4(ji,jj,ikt)*rfact2
     &    /fse3t(ji,jj,ikt)/rjjss*rivpo4input/(raass*sumsedpo4)
        trn(ji,jj,ikt,jppoc)=trn(ji,jj,ikt,jppoc)
     &    -trn(ji,jj,ikt,jppoc)*wsbio3(ji,jj,ikt)*rfact2
     &    /fse3t(ji,jj,ikt)/rjjss*rivpo4input/(raass*sumsedpo4)
        trn(ji,jj,ikt,jpbfe)=trn(ji,jj,ikt,jpbfe)
     &    -trn(ji,jj,ikt,jpbfe)*wsbio4(ji,jj,ikt)*rfact2
     &    /fse3t(ji,jj,ikt)/rjjss*rivpo4input/(raass*sumsedpo4)
        trn(ji,jj,ikt,jpsfe)=trn(ji,jj,ikt,jpsfe)
     &    -trn(ji,jj,ikt,jpsfe)*wsbio3(ji,jj,ikt)*rfact2
     &    /fse3t(ji,jj,ikt)/rjjss*rivpo4input/(raass*sumsedpo4)
          END DO
        END DO

C
C  Nitrogen fixation (simple parameterization). The total gain
C  from nitrogen fixation is scaled to balance the loss by 
C  denitrification
C  -------------------------------------------------------------
C
        denitot=0.
        DO jk=1,jpk-1
          DO jj=2,jpj-1
            DO ji=2,jpi-1
        denitot=denitot+denitr(ji,jj,jk)*rdenit*e1t(ji,jj)*e2t(ji,jj)
     &    *fse3t(ji,jj,jk)*tmask(ji,jj,jk)
            END DO
          END DO
        END DO
C
C  Potential nitrogen fication dependant on temperature
C  and iron
C  ----------------------------------------------------
C
        nitrpot(:,:)= 0.
        nitrpottot=0.
        DO jj=2,jpj-1
          DO ji=2,jpi-1
        nitrpot(ji,jj)=prmax(ji,jj,1)*max(0.,(0.1*tn(ji,jj,1)
     &    -2.))*conc0/(trn(ji,jj,1,jpno3)+conc0)*rfact2
     &    *trn(ji,jj,1,jpfer)/(conc3+trn(ji,jj,1,jpfer))
     &    *trn(ji,jj,1,jppo4)/(conc0+trn(ji,jj,1,jppo4))
        nitrpottot=nitrpottot+nitrpot(ji,jj)*e1t(ji,jj)
     &    *e2t(ji,jj)*tmask(ji,jj,1)*fse3t(ji,jj,1)
          END DO
        END DO 
C
C  Nitrogen change due to nitrogen fixation
C  ----------------------------------------
C
 

        DO jj=1,jpj
          DO ji=1,jpi
        trn(ji,jj,1,jpnh4)=trn(ji,jj,1,jpnh4)+nitrpot(ji,jj)
     &    *(denitot-rivnitinput/raass*rfact2)/(nitrpottot+rtrn)
        trn(ji,jj,1,jpoxy)=trn(ji,jj,1,jpoxy)+nitrpot(ji,jj)
     &    *(denitot-rivnitinput/raass*rfact2)/(nitrpottot+rtrn)
     &    *o2nit
          END DO
        END DO
 

C
#    if defined key_trc_diaadd
        DO jj = 1,jpj
          DO ji = 1,jpi
        trc2d(ji,jj,13) = nitrpot(ji,jj)
     &    *(denitot-rivnitinput/raass*rfact2)/(nitrpottot+rtrn)
     &    /rfact2*fse3t(ji,jj,1)
          END DO
        END DO
#    endif
C
#endif
      RETURN
      END