Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trcrad.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trcrad.F90	(revision 11021)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trcrad.F90	(revision 11434)
@@ -37,5 +37,5 @@
 CONTAINS
 
-   SUBROUTINE trc_rad( kt, Kbb, Kmm, Krhs, ptr )
+   SUBROUTINE trc_rad( kt, Kbb, Kmm, ptr )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trc_rad  ***
@@ -52,7 +52,7 @@
       !!                (the total CFC content is not strictly preserved)
       !!----------------------------------------------------------------------
-      INTEGER,                                    INTENT(in   ) :: kt              ! ocean time-step index
-      INTEGER,                                    INTENT(in   ) :: Kbb, Kmm, Krhs  ! time level indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr             ! passive tracers and RHS of tracer equation
+      INTEGER,                                    INTENT(in   ) :: kt         ! ocean time-step index
+      INTEGER,                                    INTENT(in   ) :: Kbb, Kmm   ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr        ! passive tracers and RHS of tracer equation
       !
       CHARACTER (len=22) :: charout
@@ -61,14 +61,14 @@
       IF( ln_timing )   CALL timing_start('trc_rad')
       !
-      IF( ln_age     )   CALL trc_rad_sms( kt, Kmm, Krhs, ptr(:,:,:,:,Kbb), ptr(:,:,:,:,Kmm), jp_age , jp_age                )  !  AGE
-      IF( ll_cfc     )   CALL trc_rad_sms( kt, Kmm, Krhs, ptr(:,:,:,:,Kbb), ptr(:,:,:,:,Kmm), jp_cfc0, jp_cfc1               )  !  CFC model
-      IF( ln_c14     )   CALL trc_rad_sms( kt, Kmm, Krhs, ptr(:,:,:,:,Kbb), ptr(:,:,:,:,Kmm), jp_c14 , jp_c14                )  !  C14
-      IF( ln_pisces  )   CALL trc_rad_sms( kt, Kmm, Krhs, ptr(:,:,:,:,Kbb), ptr(:,:,:,:,Kmm), jp_pcs0, jp_pcs1, cpreserv='Y' )  !  PISCES model
-      IF( ln_my_trc  )   CALL trc_rad_sms( kt, Kmm, Krhs, ptr(:,:,:,:,Kbb), ptr(:,:,:,:,Kmm), jp_myt0, jp_myt1               )  !  MY_TRC model
+      IF( ln_age     )   CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_age , jp_age                )  !  AGE
+      IF( ll_cfc     )   CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_cfc0, jp_cfc1               )  !  CFC model
+      IF( ln_c14     )   CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_c14 , jp_c14                )  !  C14
+      IF( ln_pisces  )   CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_pcs0, jp_pcs1, cpreserv='Y' )  !  PISCES model
+      IF( ln_my_trc  )   CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_myt0, jp_myt1               )  !  MY_TRC model
       !
       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=ptr(:,:,:,:,Kmm), mask=tmask, clinfo=ctrcnm )
+         CALL prt_ctl_trc( tab4d=ptr(:,:,:,:,Kbb), mask=tmask, clinfo=ctrcnm )
       ENDIF
       !
@@ -115,168 +115,115 @@
 
 
-   SUBROUTINE trc_rad_sms( kt, Kmm, Krhs, 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)
-      !!--------------------------------------------------------------------------------
-      INTEGER                                , INTENT(in   ) ::   kt                 ! ocean time-step index
-      INTEGER                                , INTENT(in   ) ::   Kmm, Krhs          ! time level indices
-      INTEGER                                , INTENT(in   ) ::   jp_sms0, jp_sms1   ! First & 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), OPTIONAL          , INTENT(in   ) ::   cpreserv           ! flag to preserve content or not
-      !
-      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  ==!
-         !
-         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 tr(:,:,:,:,Kbb) 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
-                        !
-                        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
-            END DO
-            !
-            IF( l_trdtrc ) THEN
-               ztrtrd(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt
-               CALL trd_tra( kt, Kmm, Krhs, '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 tr 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, Kmm, Krhs, '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 tr for trend computation
-            !
-            WHERE( ptrb(:,:,:,jn) < 0. )   ptrb(:,:,:,jn) = 0.
-            !
-            IF( l_trdtrc ) THEN
-               ztrtrd(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt
-               CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_radb, ztrtrd )       ! Asselin-like trend handling
-            ENDIF
-            !
-            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrn(:,:,:,jn)                        ! save input tr for trend computation
-            !
-            WHERE( ptrn(:,:,:,jn) < 0. )   ptrn(:,:,:,jn) = 0.
-            !
-            IF( l_trdtrc ) THEN
-               ztrtrd(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt
-               CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_radn, ztrtrd )       ! standard     trend handling
-            ENDIF
-            !
-         END DO
-         !
-      ENDIF
+   SUBROUTINE trc_rad_sms( kt, Kbb, Kmm, ptr, 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)
+     !!--------------------------------------------------------------------------------
+     INTEGER                                    , INTENT(in   ) ::   kt                 ! ocean time-step index
+     INTEGER                                    , INTENT(in   ) ::   Kbb, Kmm           ! time level indices
+     INTEGER                                    , INTENT(in   ) ::   jp_sms0, jp_sms1   ! First & last index of the passive tracer model
+     REAL(wp), DIMENSION (jpi,jpj,jpk,jptra,jpt), INTENT(inout) ::   ptr                ! before and now traceur concentration
+     CHARACTER( len = 1), OPTIONAL              , INTENT(in   ) ::   cpreserv           ! flag to preserve content or not
+     !
+     INTEGER ::   ji, ji2, jj, jj2, jk, jn, jt ! dummy loop indices
+     INTEGER ::   icnt, itime
+     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 ) )
+     !
+     DO jt = 1,2  ! Loop over time indices since exactly the same fix is applied to "now" and "after" fields
+        IF( jt == 1 ) itime = Kbb
+        IF( jt == 2 ) itime = Kmm
+
+        IF( PRESENT( cpreserv )  ) THEN     !==  total tracer concentration is preserved  ==!
+           !
+           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., ptr(:,:,:,jn,itime) ) * cvol(:,:,:), dim = 3 )   ! sum of the negative values
+              ztrpos(:,:,jn) = SUM( MAX( 0., ptr(:,:,:,jn,itime) ) * 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(:,:,:) = ptr(:,:,:,jn,itime)                       ! save input tr(:,:,:,:,Kbb) 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
+                          !
+                          ptr(ji,jj,jk,jn,itime) = ptr(ji,jj,jk,jn,itime) * tmask(ji,jj,jk)   ! really needed?
+                          IF( ptr(ji,jj,jk,jn,itime) < 0. ) ptr(ji,jj,jk,jn,itime) = 0.       ! suppress negative values
+                          IF( ptr(ji,jj,jk,jn,itime) > 0. ) THEN                    ! use positive values to compensate mass gain
+                             zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn)       ! ztrpos > 0 as ptr > 0
+                             ptr(ji,jj,jk,jn,itime) = ptr(ji,jj,jk,jn,itime) * zcoef
+                             IF( zcoef < 0. ) THEN                                  ! if the compensation exceed the positive value
+                                gainmass(jn,1) = gainmass(jn,1) - ptr(ji,jj,jk,jn,itime) * cvol(ji,jj,jk)   ! we are adding mass...
+                                ptr(ji,jj,jk,jn,itime) = 0.                         ! limit the compensation to keep positive value
+                             ENDIF
+                          ENDIF
+                          !
+                       ENDIF
+                    END DO
+                 END DO
+              END DO
+              !
+              IF( l_trdtrc ) THEN
+                 ztrtrd(:,:,:) = ( ptr(:,:,:,jn,itime) - ztrtrd(:,:,:) ) * zs2rdt
+                 CALL trd_tra( kt, Kbb, Kmm, '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', ptr(:,:,:,jn,itime) * 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
+
+           DEALLOCATE( ztrneg, ztrpos )
+           !
+        ELSE                                !==  total CFC content is NOT strictly preserved  ==!
+           !
+           DO jn = jp_sms0, jp_sms1  
+              !
+              IF( l_trdtrc )   ztrtrd(:,:,:) = ptr(:,:,:,jn,itime)                 ! save input tr for trend computation
+              !
+              WHERE( ptr(:,:,:,jn,itime) < 0. )   ptr(:,:,:,jn,itime) = 0.
+              !
+              IF( l_trdtrc ) THEN
+                 ztrtrd(:,:,:) = ( ptr(:,:,:,jn,itime) - ztrtrd(:,:,:) ) * zs2rdt
+                 CALL trd_tra( kt, Kbb, Kmm, 'TRC', jn, jptra_radb, ztrtrd )       ! Asselin-like trend handling
+              ENDIF
+              !
+           END DO
+           !
+        ENDIF
+        !
+      END DO
       !
       IF( l_trdtrc )  DEALLOCATE( ztrtrd )
@@ -289,7 +236,7 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE trc_rad( kt, Kbb, Kmm, Krhs )              ! Empty routine
+   SUBROUTINE trc_rad( kt, Kbb, Kmm )              ! Empty routine
       INTEGER, INTENT(in) ::   kt
-      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs  ! time level indices
+      INTEGER, INTENT(in) ::   Kbb, Kmm  ! time level indices
       WRITE(*,*) 'trc_rad: You should not have seen this print! error?', kt
    END SUBROUTINE trc_rad
