Index: /NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/lib_fortran.F90
===================================================================
--- /NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/lib_fortran.F90	(revision 10407)
+++ /NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/lib_fortran.F90	(revision 10408)
@@ -21,4 +21,5 @@
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! distributed memory computing
+   USE lbclnk          ! ocean lateral boundary conditions
 
    IMPLICIT NONE
@@ -28,4 +29,5 @@
    PUBLIC   glob_sum_full ! used in many places (masked with tmask_h, ie only over the halos)
    PUBLIC   local_sum     ! used in trcrad, local operation before glob_sum_delay
+   PUBLIC   sum3x3        ! used in trcrad, do a sum over 3x3 boxes
    PUBLIC   DDPDD         ! also used in closea module
    PUBLIC   glob_min, glob_max
@@ -42,4 +44,7 @@
    INTERFACE local_sum
       MODULE PROCEDURE local_sum_2d, local_sum_3d
+   END INTERFACE
+   INTERFACE sum3x3
+      MODULE PROCEDURE sum3x3_2d, sum3x3_3d
    END INTERFACE
    INTERFACE glob_min
@@ -194,4 +199,102 @@
    END FUNCTION local_sum_3d
 
+!                          ! FUNCTION sum3x3 !
+
+   SUBROUTINE sum3x3_2d( p2d )
+      !!-----------------------------------------------------------------------
+      !!                  ***  routine sum3x3_2d  ***
+      !!
+      !! ** Purpose : sum over 3x3 boxes
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION (:,:), INTENT(inout) ::   p2d
+      !
+      INTEGER ::   ji, ji2, jj, jj2     ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      IF( SIZE(p2d,1) /= jpi ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_2d, the first dimension is not equal to jpi' ) 
+      IF( SIZE(p2d,2) /= jpj ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_2d, the second dimension is not equal to jpj' ) 
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi 
+            IF( MOD(mig(ji), 3) == 1 .AND. MOD(mjg(jj), 3) == 1 ) THEN   ! bottom left corber of a 3x3 box
+               ji2 = MIN(mig(ji)+2, jpiglo) - nimpp + 1                  ! right position of the box
+               jj2 = MIN(mjg(jj)+2, jpjglo) - njmpp + 1                  ! upper position of the box
+               IF( ji2 <= jpi .AND. jj2 <= jpj ) THEN                    ! the box is fully included in the local mpi domain
+                  p2d(ji:ji2,jj:jj2) = SUM(p2d(ji:ji2,jj:jj2))
+               ENDIF
+            ENDIF
+         END DO
+      END DO
+      CALL lbc_lnk( 'lib_fortran', p2d, 'T', 1. )
+      IF( nbondi /= -1 ) THEN
+         IF( MOD(mig(    1), 3) == 1 )   p2d(    1,:) = p2d(    2,:)
+         IF( MOD(mig(    1), 3) == 2 )   p2d(    2,:) = p2d(    1,:)
+      ENDIF
+      IF( nbondi /=  1 ) THEN
+         IF( MOD(mig(jpi-2), 3) == 1 )   p2d(  jpi,:) = p2d(jpi-1,:)
+         IF( MOD(mig(jpi-2), 3) == 0 )   p2d(jpi-1,:) = p2d(  jpi,:)
+      ENDIF
+      IF( nbondj /= -1 ) THEN
+         IF( MOD(mjg(    1), 3) == 1 )   p2d(:,    1) = p2d(:,    2)
+         IF( MOD(mjg(    1), 3) == 2 )   p2d(:,    2) = p2d(:,    1)
+      ENDIF
+      IF( nbondj /=  1 ) THEN
+         IF( MOD(mjg(jpj-2), 3) == 1 )   p2d(:,  jpj) = p2d(:,jpj-1)
+         IF( MOD(mjg(jpj-2), 3) == 0 )   p2d(:,jpj-1) = p2d(:,  jpj)
+      ENDIF
+      CALL lbc_lnk( 'lib_fortran', p2d, 'T', 1. )
+
+   END SUBROUTINE sum3x3_2d
+
+   SUBROUTINE sum3x3_3d( p3d )
+      !!-----------------------------------------------------------------------
+      !!                  ***  routine sum3x3_3d  ***
+      !!
+      !! ** Purpose : sum over 3x3 boxes
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   p3d
+      !
+      INTEGER ::   ji, ji2, jj, jj2, jn     ! dummy loop indices
+      INTEGER ::   ipn                      ! Third dimension size
+      !!----------------------------------------------------------------------
+      !
+      IF( SIZE(p3d,1) /= jpi ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_3d, the first dimension is not equal to jpi' ) 
+      IF( SIZE(p3d,2) /= jpj ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_3d, the second dimension is not equal to jpj' ) 
+      ipn = SIZE(p3d,3)
+      !
+      DO jn = 1, ipn
+         DO jj = 1, jpj
+            DO ji = 1, jpi 
+               IF( MOD(mig(ji), 3) == 1 .AND. MOD(mjg(jj), 3) == 1 ) THEN   ! bottom left corber of a 3x3 box
+                  ji2 = MIN(mig(ji)+2, jpiglo) - nimpp + 1                  ! right position of the box
+                  jj2 = MIN(mjg(jj)+2, jpjglo) - njmpp + 1                  ! upper position of the box
+                  IF( ji2 <= jpi .AND. jj2 <= jpj ) THEN                    ! the box is fully included in the local mpi domain
+                     p3d(ji:ji2,jj:jj2,jn) = SUM(p3d(ji:ji2,jj:jj2,jn))
+                  ENDIF
+               ENDIF
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( 'lib_fortran', p3d, 'T', 1. )
+      IF( nbondi /= -1 ) THEN
+         IF( MOD(mig(    1), 3) == 1 )   p3d(    1,:,:) = p3d(    2,:,:)
+         IF( MOD(mig(    1), 3) == 2 )   p3d(    2,:,:) = p3d(    1,:,:)
+      ENDIF
+      IF( nbondi /=  1 ) THEN
+         IF( MOD(mig(jpi-2), 3) == 1 )   p3d(  jpi,:,:) = p3d(jpi-1,:,:)
+         IF( MOD(mig(jpi-2), 3) == 0 )   p3d(jpi-1,:,:) = p3d(  jpi,:,:)
+      ENDIF
+      IF( nbondj /= -1 ) THEN
+         IF( MOD(mjg(    1), 3) == 1 )   p3d(:,    1,:) = p3d(:,    2,:)
+         IF( MOD(mjg(    1), 3) == 2 )   p3d(:,    2,:) = p3d(:,    1,:)
+      ENDIF
+      IF( nbondj /=  1 ) THEN
+         IF( MOD(mjg(jpj-2), 3) == 1 )   p3d(:,  jpj,:) = p3d(:,jpj-1,:)
+         IF( MOD(mjg(jpj-2), 3) == 0 )   p3d(:,jpj-1,:) = p3d(:,  jpj,:)
+      ENDIF
+      CALL lbc_lnk( 'lib_fortran', p3d, 'T', 1. )
+
+   END SUBROUTINE sum3x3_3d
+
 
    SUBROUTINE DDPDD( ydda, yddb )
Index: /NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/TRP/trcrad.F90
===================================================================
--- /NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/TRP/trcrad.F90	(revision 10407)
+++ /NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/TRP/trcrad.F90	(revision 10408)
@@ -19,4 +19,5 @@
    USE trdtra
    USE prtctl_trc          ! Print control for debbuging
+   USE lib_fortran
 
    IMPLICIT NONE
@@ -27,4 +28,5 @@
 
    LOGICAL , PUBLIC ::   ln_trcrad           !: flag to artificially correct negative concentrations
+   REAL(wp), DIMENSION(:,:), ALLOCATABLE::   gainmass
 
    !!----------------------------------------------------------------------
@@ -104,4 +106,7 @@
          ENDIF
       ENDIF
+      !
+      ALLOCATE( gainmass(jptra,2) )
+      gainmass(:,:) = 0.
       !
    END SUBROUTINE trc_rad_ini
@@ -129,72 +134,47 @@
       CHARACTER( len = 1), OPTIONAL          , INTENT(in   ) ::   cpreserv           ! flag to preserve content or not
       !
-      INTEGER ::   ji, jj, jk, jn     ! dummy loop indices
-      LOGICAL ::   lldebug = .FALSE.           ! local logical
-      REAL(wp)::   ztrcorb, ztrmasb, zs2rdt    ! temporary scalars
-      REAL(wp)::   zcoef  , ztrcorn, ztrmasn   !    -         -
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrtrdb, ztrtrdn   ! workspace arrays
-      !!----------------------------------------------------------------------
-      !
-      IF( l_trdtrc )   ALLOCATE( ztrtrdb(jpi,jpj,jpk), ztrtrdn(jpi,jpj,jpk) )
+      INTEGER ::   ji, ji2, jj, jj2, jk, jn     ! dummy loop indices
+      INTEGER ::   icnt
+      LOGICAL ::   lldebug = .FALSE.            ! local logical
+      REAL(wp)::   zcoef, zs2rdt, ztotmass
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrneg, ztrpos
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrtrd   ! workspace arrays
+      !!----------------------------------------------------------------------
+      !
+      IF( l_trdtrc )   ALLOCATE( ztrtrd(jpi,jpj,jpk) )
+      zs2rdt = 1. / ( 2. * rdt * REAL( nn_dttrc, wp ) )
       !
       IF( PRESENT( cpreserv )  ) THEN     !==  total tracer concentration is preserved  ==!
          !
-         DO jn = jp_sms0, jp_sms1
-            !
-            ztrcorb = 0._wp   ;   ztrmasb = 0._wp
-            ztrcorn = 0._wp   ;   ztrmasn = 0._wp
-            !
-            IF( l_trdtrc ) THEN
-               ztrtrdb(:,:,:) = ptrb(:,:,:,jn)                        ! save input trb for trend computation
-               ztrtrdn(:,:,:) = ptrn(:,:,:,jn)                        ! save input trn for trend computation
-            ENDIF
-            !                                                         ! sum over the global domain 
-            ztrcorb = glob_sum( 'trcrad', MIN( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:) )
-            ztrcorn = glob_sum( 'trcrad', MIN( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:) )
-            !
-            IF( ztrcorb /= 0 ) THEN
-               ztrmasb = glob_sum( 'trcrad', MAX( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:) )
-               zcoef = 1. + ztrcorb / ztrmasb
-               DO jk = 1, jpkm1
-                  ptrb(:,:,jk,jn) = MAX( 0., ptrb(:,:,jk,jn) )
-                  ptrb(:,:,jk,jn) = ptrb(:,:,jk,jn) * zcoef * tmask(:,:,jk)
-               END DO
-            ENDIF
-            !
-            IF( ztrcorn /= 0 ) THEN
-               ztrmasn = glob_sum( 'trcrad', MAX( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:) )
-               zcoef = 1. + ztrcorn / ztrmasn
-               DO jk = 1, jpkm1
-                  ptrn(:,:,jk,jn) = MAX( 0., ptrn(:,:,jk,jn) )
-                  ptrn(:,:,jk,jn) = ptrn(:,:,jk,jn) * zcoef * tmask(:,:,jk)
-               END DO
-            ENDIF
-            !
-            IF( l_trdtrc ) THEN
-               !
-               zs2rdt = 1. / ( 2. * rdt )
-               ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt
-               ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt 
-               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb )       ! Asselin-like trend handling
-               CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn )       ! standard     trend handling
-              !
-            ENDIF
-            !
-         END DO
-         !
-      ELSE                                !==  total CFC content is NOT strictly preserved  ==!
-         !
-         DO jn = jp_sms0, jp_sms1  
-            !
-            IF( l_trdtrc ) THEN
-               ztrtrdb(:,:,:) = ptrb(:,:,:,jn)                        ! save input trb for trend computation
-               ztrtrdn(:,:,:) = ptrn(:,:,:,jn)                        ! save input trn for trend computation
-            ENDIF
+         ALLOCATE( ztrneg(1:jpi,1:jpj,jp_sms0:jp_sms1), ztrpos(1:jpi,1:jpj,jp_sms0:jp_sms1) )
+
+         DO jn = jp_sms0, jp_sms1
+            ztrneg(:,:,jn) = SUM( MIN( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:), dim = 3 )   ! sum of the negative values
+            ztrpos(:,:,jn) = SUM( MAX( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:), dim = 3 )   ! sum of the positive values
+         END DO
+         CALL sum3x3( ztrneg )
+         CALL sum3x3( ztrpos )
+         
+         DO jn = jp_sms0, jp_sms1
+            !
+            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrb(:,:,:,jn)                            ! save input trb for trend computation           
             !
             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) )
+                     IF( ztrneg(ji,jj,jn) /= 0. ) THEN                                 ! if negative values over the 3x3 box
+                        !
+                        ptrb(ji,jj,jk,jn) = ptrb(ji,jj,jk,jn) * tmask(ji,jj,jk)   ! really needed?
+                        IF( ptrb(ji,jj,jk,jn) < 0. ) ptrb(ji,jj,jk,jn) = 0.       ! supress negative values
+                        IF( ptrb(ji,jj,jk,jn) > 0. ) THEN                         ! use positive values to compensate mass gain
+                           zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn)       ! ztrpos > 0 as ptrb > 0
+                           ptrb(ji,jj,jk,jn) = ptrb(ji,jj,jk,jn) * zcoef
+                           IF( zcoef < 0. ) THEN                                  ! if the compensation exceed the positive value
+                              gainmass(jn,1) = gainmass(jn,1) - ptrb(ji,jj,jk,jn) * cvol(ji,jj,jk)   ! we are adding mass...
+                              ptrb(ji,jj,jk,jn) = 0.                              ! limit the compensation to keep positive value
+                           ENDIF
+                        ENDIF
+                        !
+                     ENDIF
                   END DO
                END DO
@@ -202,11 +182,93 @@
             !
             IF( l_trdtrc ) THEN
-               !
-               zs2rdt = 1. / ( 2. * rdt * REAL( nn_dttrc, wp ) )
-               ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt
-               ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt 
-               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb )       ! Asselin-like trend handling
-               CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn )       ! standard     trend handling
-              !
+               ztrtrd(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt
+               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrd )       ! Asselin-like trend handling
+            ENDIF
+            !
+         END DO
+ 
+         IF( kt == nitend ) THEN
+            CALL mpp_sum( 'trcrad', gainmass(:,1) )
+            DO jn = jp_sms0, jp_sms1
+               IF( gainmass(jn,1) > 0. ) THEN
+                  ztotmass = glob_sum( 'trcrad', ptrb(:,:,:,jn) * cvol(:,:,:) )
+                  IF(lwp) WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrb, traceur ', jn  &
+                     &        , ' total mass : ', ztotmass, ', mass gain : ',  gainmass(jn,1)
+               END IF
+            END DO
+         ENDIF
+
+         DO jn = jp_sms0, jp_sms1
+            ztrneg(:,:,jn) = SUM( MIN( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:), dim = 3 )   ! sum of the negative values
+            ztrpos(:,:,jn) = SUM( MAX( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:), dim = 3 )   ! sum of the positive values
+         END DO
+         CALL sum3x3( ztrneg )
+         CALL sum3x3( ztrpos )
+         
+         DO jn = jp_sms0, jp_sms1
+            !
+            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrn(:,:,:,jn)                            ! save input trb for trend computation
+            !
+            DO jk = 1, jpkm1
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     IF( ztrneg(ji,jj,jn) /= 0. ) THEN                                 ! if negative values over the 3x3 box
+                        !
+                        ptrn(ji,jj,jk,jn) = ptrn(ji,jj,jk,jn) * tmask(ji,jj,jk)   ! really needed?
+                        IF( ptrn(ji,jj,jk,jn) < 0. ) ptrn(ji,jj,jk,jn) = 0.       ! supress negative values
+                        IF( ptrn(ji,jj,jk,jn) > 0. ) THEN                         ! use positive values to compensate mass gain
+                           zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn)       ! ztrpos > 0 as ptrb > 0
+                           ptrn(ji,jj,jk,jn) = ptrn(ji,jj,jk,jn) * zcoef
+                           IF( zcoef < 0. ) THEN                                  ! if the compensation exceed the positive value
+                              gainmass(jn,2) = gainmass(jn,2) - ptrn(ji,jj,jk,jn) * cvol(ji,jj,jk)   ! we are adding mass...
+                              ptrn(ji,jj,jk,jn) = 0.                              ! limit the compensation to keep positive value
+                           ENDIF
+                        ENDIF
+                        !
+                     ENDIF
+                  END DO
+               END DO
+            END DO
+            !
+            IF( l_trdtrc ) THEN
+               ztrtrd(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt
+               CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrd )       ! standard     trend handling
+            ENDIF
+            !
+         END DO
+ 
+         IF( kt == nitend ) THEN
+            CALL mpp_sum( 'trcrad', gainmass(:,2) )
+            DO jn = jp_sms0, jp_sms1
+               IF( gainmass(jn,2) > 0. ) THEN
+                  ztotmass = glob_sum( 'trcrad', ptrn(:,:,:,jn) * cvol(:,:,:) )
+                  WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrn, traceur ', jn  &
+                     &        , ' total mass : ', ztotmass, ', mass gain : ',  gainmass(jn,1)
+               END IF
+            END DO
+         ENDIF
+
+         DEALLOCATE( ztrneg, ztrpos )
+         !
+      ELSE                                !==  total CFC content is NOT strictly preserved  ==!
+         !
+         DO jn = jp_sms0, jp_sms1  
+            !
+            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrb(:,:,:,jn)                        ! save input trb for trend computation
+            !
+            WHERE( ptrb(:,:,:,jn) < 0. )   ptrb(:,:,:,jn) = 0.
+            !
+            IF( l_trdtrc ) THEN
+               ztrtrd(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt
+               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrd )       ! Asselin-like trend handling
+            ENDIF
+            !
+            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrn(:,:,:,jn)                        ! save input trn for trend computation
+            !
+            WHERE( ptrn(:,:,:,jn) < 0. )   ptrn(:,:,:,jn) = 0.
+            !
+            IF( l_trdtrc ) THEN
+               ztrtrd(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt
+               CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrd )       ! standard     trend handling
             ENDIF
             !
@@ -215,5 +277,5 @@
       ENDIF
       !
-      IF( l_trdtrc )  DEALLOCATE( ztrtrdb, ztrtrdn )
+      IF( l_trdtrc )  DEALLOCATE( ztrtrd )
       !
    END SUBROUTINE trc_rad_sms
