Index: trunk/NEMO/TOP_SRC/TRP/trcrad.F90
===================================================================
--- trunk/NEMO/TOP_SRC/TRP/trcrad.F90	(revision 941)
+++ trunk/NEMO/TOP_SRC/TRP/trcrad.F90	(revision 1003)
@@ -49,8 +49,4 @@
       !!----------------------------------------------------------------------
       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
       !!----------------------------------------------------------------------
@@ -62,23 +58,63 @@
       ENDIF
 
-      IF( lk_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_cfc     )   CALL trc_rad_sms( trb, trn, jp_cfc0, jp_cfc1               ) ! CFC model
+      IF( lk_lobster )   CALL trc_rad_sms( trb, trn, jp_lob0, jp_lob1, cpreserv='Y' )  ! LOBSTER model
+      IF( lk_pisces  )   CALL trc_rad_sms( trb, trn, jp_pcs0, jp_pcs1, cpreserv='Y' )  ! PISCES model
+      IF( lk_my_trc  )   CALL trc_rad_sms( trb, trn, jp_myt0, jp_myt1               ) ! MY_TRC model
 
-      IF( lk_pisces .OR. lk_lobster ) THEN      ! PISCES or LOBSTER bio-model
-         DO jn = 1, jptra
+
+      !
+      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
+
+   SUBROUTINE trc_rad_sms( ptrb, ptrn, jp_sms0, jp_sms1, cpreserv )
+      !!-----------------------------------------------------------------------------
+      !!                  ***  ROUTINE trc_rad_sms  ***
+      !!
+      !! ** Purpose :   "crappy" routine to correct artificial negative
+      !!              concentrations due to isopycnal scheme
+      !!
+      !! ** Method  : 2 cases :
+      !!                - 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.
+      !!                - simply set to zero the negative CFC concentration
+      !!                  (the total content of concentration is not strictly preserved)
+      !!--------------------------------------------------------------------------------
+      !! Arguments
+      INTEGER  , INTENT( in ) ::  &
+         jp_sms0, &       !: First index of the passive tracer model
+         jp_sms1          !: Last  index of  the passive tracer model
+
+      REAL(wp), DIMENSION (jpi,jpj,jpk,jptra), INTENT( inout )  :: &
+         ptrb, ptrn       !: before and now traceur concentration
+
+      CHARACTER( len = 1) , INTENT(in), OPTIONAL  :: &
+         cpreserv          !: flag to preserve content or not
+      
+      ! Local declarations
+      INTEGER  ::  ji, jj, jk, jn     ! dummy loop indices
+      REAL(wp) :: zvolk, ztrcorb, ztrmasb   ! temporary scalars
+      REAL(wp) :: zcoef, ztrcorn, ztrmasn   !    "         "
+
+      !!----------------------------------------------------------------------
+
+      
+      IF( PRESENT( cpreserv )  ) THEN   !  total tracer concentration is preserved 
+      
+         DO jn = jp_sms0, jp_sms1
+            
             ztrcorb = 0.e0
             ztrmasb = 0.e0
             ztrcorn = 0.e0
             ztrmasn = 0.e0
+
             DO jk = 1, jpkm1
                DO jj = 1, jpj
@@ -90,15 +126,16 @@
                         &   * 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
+                     ztrcorb = ztrcorb + MIN( 0., ptrb(ji,jj,jk,jn) )  * zvolk
+                     ztrcorn = ztrcorn + MIN( 0., ptrn(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) )
+                     ptrb(ji,jj,jk,jn) = MAX( 0. , ptrb(ji,jj,jk,jn) )
+                     ptrn(ji,jj,jk,jn) = MAX( 0. , ptrn(ji,jj,jk,jn) )
 
-                     ztrmasb = ztrmasb + trb(ji,jj,jk,jn) * zvolk
-                     ztrmasn = ztrmasn + trn(ji,jj,jk,jn) * zvolk
+                     ztrmasb = ztrmasb + ptrb(ji,jj,jk,jn) * zvolk
+                     ztrmasn = ztrmasn + ptrn(ji,jj,jk,jn) * zvolk
                   END DO
                END DO
             END DO
+
             IF( lk_mpp ) THEN
                CALL mpp_sum( ztrcorb )      ! sum over the global domain
@@ -111,5 +148,5 @@
                zcoef = 1. + ztrcorb / ztrmasb
                DO jk = 1, jpkm1
-                  trb(:,:,jk,jn) = trb(:,:,jk,jn) * zcoef * tmask(:,:,jk)
+                  ptrb(:,:,jk,jn) = ptrb(:,:,jk,jn) * zcoef * tmask(:,:,jk)
                END DO
             ENDIF
@@ -118,5 +155,5 @@
                zcoef = 1. + ztrcorn / ztrmasn
                DO jk = 1, jpkm1
-                  trn(:,:,jk,jn) = trn(:,:,jk,jn) * zcoef * tmask(:,:,jk)
+                  ptrn(:,:,jk,jn) = ptrn(:,:,jk,jn) * zcoef * tmask(:,:,jk)
                END DO
             ENDIF
@@ -124,14 +161,21 @@
          END DO
          !
+         !
+      ELSE  ! total CFC content is not strictly preserved
+
+         DO jn = jp_sms0, jp_sms1  
+            DO jk = 1, jpkm1
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ptrn(ji,jj,jk,jn) = MAX( 0. , ptrn(ji,jj,jk,jn) )
+                     ptrb(ji,jj,jk,jn) = MAX( 0. , ptrb(ji,jj,jk,jn) )
+                  END DO
+               END DO
+            END DO
+         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
 
+   END SUBROUTINE trc_rad_sms
 #else
    !!----------------------------------------------------------------------
