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 sms_pisces      ! 
   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

   REAL(wp) ::   &
     xstep            !: Time step duration for biology

   !!* Substitution
#  include "top_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
   !! $Id: p4zrem.F90 1808 2010-03-11 09:17:56Z cetlod $ 
   !! 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
      REAL(wp) ::   zosil, zdenom1, zscave, zaggdfe
#if ! defined key_kriest
      REAL(wp) ::   zofer2, zdenom, zdenom2
#endif
      REAL(wp) ::   zlamfac, zonitr
      REAL(wp), DIMENSION(jpi,jpj)     ::   ztempbac
      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdepbac, zfesatur, zolimi
      CHARACTER (len=25) :: charout

      !!---------------------------------------------------------------------


      IF( ( kt * jnt ) == nittrc000  )  THEN
         CALL p4z_rem_init                ! Initialization (first time-step only)
         xstep = rfact2 / rday            ! Time step duration for the biology
         nitrfac(:,:,:) = 0.0
         denitr (:,:,:) = 0.0  
      ENDIF


       ! Initialisation of temprary arrys
       zdepbac (:,:,:) = 0.0
       zfesatur(:,:,:) = 0.0
       zolimi  (:,:,:) = 0.0
       ztempbac(:,:)   = 0.0

!      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 * xstep / 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 * xstep )

!     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

      DO jk = 1, jpkm1
         DO jj = 1, jpj
            DO ji = 1, jpi
               zolimi (ji,jj,jk) = MAX( 0.e0, zolimi (ji,jj,jk) )
               denitr (ji,jj,jk) = MAX( 0.e0, denitr (ji,jj,jk) )
            END DO
         END DO
      END DO

      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 * xstep * 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 * xstep * 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 * xstep * znusil * facvol(ji,jj,jk)
# else
               zsiremin = xsirem * xstep * 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 * xstep  * facvol(ji,jj,jk)
# else
               zscave = zfeequi * zlam1b * xstep
# 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 * xstep * 0.5 * ( trn(ji,jj,jk,jpfer) - zfeequi ) * facvol(ji,jj,jk)
# else
               zaggdfe = zlam1b * xstep * 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 nampisrem namelist and check the parameters
      !!      called at the first timestep (nittrc000)
      !!
      !! ** input   :   Namelist nampisrem
      !!
      !!----------------------------------------------------------------------

      NAMELIST/nampisrem/ xremik, xremip, nitrif, xsirem, xlam1, oxymin

      REWIND( numnat )                     ! read numnat
      READ  ( numnat, nampisrem )

      IF(lwp) THEN                         ! control print
         WRITE(numout,*) ' '
         WRITE(numout,*) ' Namelist parameters for remineralization, nampisrem'
         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
