MODULE trcrad
   !!======================================================================
   !!                       ***  MODULE  trcrad  ***
   !! Ocean passive tracers:  correction of negative concentrations
   !!======================================================================
   !! History :   -   !  01-01  (O. Aumont & E. Kestenare)  Original code
   !!            1.0  !  04-03  (C. Ethe)  free form F90
   !!----------------------------------------------------------------------
#if defined key_passivetrc
   !!----------------------------------------------------------------------
   !!   'key_passivetrc'                                    Passive tracers
   !!----------------------------------------------------------------------
   !!   trc_rad    : correction of negative concentrations
   !!----------------------------------------------------------------------
   USE oce_trc             ! ocean dynamics and tracers variables
   USE trc                 ! ocean passive tracers variables
   USE lib_mpp
   USE prtctl_trc          ! Print control for debbuging

   IMPLICIT NONE
   PRIVATE

   PUBLIC trc_rad        ! routine called by trcstp.F90

   !! * Substitutions
#  include "passivetrc_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/TOP 1.0 , LOCEAN-IPSL (2005) 
   !! $Id:$ 
   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
   
CONTAINS

   SUBROUTINE trc_rad( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE trc_rad  ***
      !!
      !! ** Purpose :   "crappy" routine to correct artificial negative
      !!              concentrations due to isopycnal scheme
      !!
      !! ** Method  : - PISCES or LOBSTER: Set negative concentrations to zero
      !!                while computing the corresponding tracer content that
      !!                is added to the tracers. Then, adjust the tracer 
      !!                concentration using a multiplicative factor so that 
      !!                the total tracer concentration is preserved.
      !!              - CFC: simply set to zero the negative CFC concentration
      !!                (the total CFC content is not strictly preserved)
      !!----------------------------------------------------------------------
      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index      
      !!
      INTEGER  ::  ji, jj, jk, jn     ! dummy loop indices
      REAL(wp) :: zvolk, ztrcorb, ztrmasb   ! temporary scalars
      REAL(wp) :: zcoef, ztrcorn, ztrmasn   !    "         "
      CHARACTER (len=22) :: charout
      !!----------------------------------------------------------------------

      IF( kt == nittrc000 ) THEN
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) 'trc_rad : Correct artificial negative concentrations '
         IF(lwp) WRITE(numout,*) '~~~~~~~ '
      ENDIF

      IF( lk_trc_cfc ) THEN                             ! CFC model
         DO jn = 1, jptra
            DO jk = 1, jpkm1
               DO jj = 1, jpj
                  DO ji = 1, jpi
                     trn(ji,jj,jk,jn) = MAX( 0. , trn(ji,jj,jk,jn) )
                     trb(ji,jj,jk,jn) = MAX( 0. , trb(ji,jj,jk,jn) )
                  END DO
               END DO
            END DO
         END DO
      ENDIF 

      IF( lk_trc_pisces .OR. lk_trc_lobster ) THEN      ! PISCES or LOBSTER bio-model
         DO jn = 1, jptra
            ztrcorb = 0.e0
            ztrmasb = 0.e0
            ztrcorn = 0.e0
            ztrmasn = 0.e0
            DO jk = 1, jpkm1
               DO jj = 1, jpj
                  DO ji = 1, jpi
                     zvolk  = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)   &
# if defined key_off_degrad
                        &   * facvol(ji,jj,jk)   &
# endif
                        &   * tmask(ji,jj,jk) * tmask_i(ji,jj)

                     ztrcorb = ztrcorb + MIN( 0., trb(ji,jj,jk,jn) )  * zvolk
                     ztrcorn = ztrcorn + MIN( 0., trn(ji,jj,jk,jn) )  * zvolk

                     trb(ji,jj,jk,jn) = MAX( 0. , trb(ji,jj,jk,jn) )
                     trn(ji,jj,jk,jn) = MAX( 0. , trn(ji,jj,jk,jn) )

                     ztrmasb = ztrmasb + trb(ji,jj,jk,jn) * zvolk
                     ztrmasn = ztrmasn + trn(ji,jj,jk,jn) * zvolk
                  END DO
               END DO
            END DO
            IF( lk_mpp ) THEN
               CALL mpp_sum( ztrcorb )      ! sum over the global domain
               CALL mpp_sum( ztrcorn )      ! sum over the global domain
               CALL mpp_sum( ztrmasb )      ! sum over the global domain
               CALL mpp_sum( ztrmasn )      ! sum over the global domain
            ENDIF

            IF( ztrcorb /= 0 ) THEN
               zcoef = 1. + ztrcorb / ztrmasb
               DO jk = 1, jpkm1
                  DO jj = 1, jpj
                     DO ji = 1, jpi
                         trb(ji,jj,jk,jn) = MAX( 0. , trb(ji,jj,jk,jn) )               !!gm bug already done just above
                        trb(ji,jj,jk,jn) = trb(ji,jj,jk,jn) * zcoef * tmask(ji,jj,jk)
                     END DO
                  END DO
               END DO
            ENDIF

            IF( ztrcorn /= 0 ) THEN
               zcoef = 1. + ztrcorn / ztrmasn
               DO jk = 1, jpkm1
                  DO jj = 1, jpj
                     DO ji = 1, jpi
                        trn(ji,jj,jk,jn) = MAX( 0. , trn(ji,jj,jk,jn) )                !!gm bug already done just above
                        trn(ji,jj,jk,jn) = trn(ji,jj,jk,jn) * zcoef * tmask(ji,jj,jk)
                     END DO
                  END DO
               END DO
            ENDIF
            !
         END DO
         !
      ENDIF
      !
      IF(ln_ctl) THEN      ! print mean trends (used for debugging)
         WRITE(charout, FMT="('rad')")
         CALL prt_ctl_trc_info( charout )
         CALL prt_ctl_trc( tab4d=trn, mask=tmask, clinfo=ctrcnm )
      ENDIF
      !
   END SUBROUTINE trc_rad

#else
   !!----------------------------------------------------------------------
   !!   Dummy module :                                         NO TOP model
   !!----------------------------------------------------------------------
CONTAINS
   SUBROUTINE trc_rad( kt )              ! Empty routine
      INTEGER, INTENT(in) ::   kt
      WRITE(*,*) 'trc_rad: You should not have seen this print! error?', kt
   END SUBROUTINE trc_rad
#endif
   
   !!======================================================================
END MODULE trcrad