Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/cfgs/SHARED/field_def_nemo-oce.xml	(revision 11056)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/cfgs/SHARED/field_def_nemo-oce.xml	(revision 11057)
@@ -62,14 +62,4 @@
         <field id="mldr10_1max"  long_name="Max of Mixed Layer Depth (dsigma = 0.01 wrt 10m)"   field_ref="mldr10_1"   operation="maximum"                                                                          />
         <field id="mldr10_1min"  long_name="Min of Mixed Layer Depth (dsigma = 0.01 wrt 10m)"   field_ref="mldr10_1"   operation="minimum"                                                                          />
-         <field id="mldzint_1"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          />
-         <field id="mldzint_2"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          />
-         <field id="mldzint_3"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          />
-         <field id="mldzint_4"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          />
-         <field id="mldzint_5"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          />
-         <field id="mldhtc_1"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       />
-         <field id="mldhtc_2"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       />
-         <field id="mldhtc_3"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       />
-         <field id="mldhtc_4"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       />
-         <field id="mldhtc_5"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       />
         <field id="heatc"        long_name="Heat content vertically integrated"                 standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       />
         <field id="saltc"        long_name="Salt content vertically integrated"                                                                                                                   unit="1e-3*kg/m2" />
@@ -87,4 +77,10 @@
         <!-- variables available with MLE -->
         <field id="Lf_NHpf"      long_name="MLE: Lf = N H / f"   unit="m" />
+
+        <!-- next variables available with ln_zad_Aimp=.true. -->
+        <field id="Courant"      long_name="Courant number"                                                                 unit="#"   grid_ref="grid_T_3D" />
+        <field id="wimp"         long_name="Implicit vertical velocity"                                                     unit="m/s" grid_ref="grid_T_3D" />
+        <field id="wexp"         long_name="Explicit vertical velocity"                                                     unit="m/s" grid_ref="grid_T_3D" />
+        <field id="wi_cff"       long_name="Fraction of implicit vertical velocity"                                         unit="#"   grid_ref="grid_T_3D" />
 
         <!-- next variables available with key_diahth -->
@@ -355,5 +351,4 @@
         <field id="uoce"         long_name="ocean current along i-axis"                             standard_name="sea_water_x_velocity"        unit="m/s"        grid_ref="grid_U_3D" />
         <field id="uoce_e3u"     long_name="ocean current along i-axis  (thickness weighted)"                                                   unit="m/s"        grid_ref="grid_U_3D"  > uoce * e3u </field>
-        <field id="uoce2_e3u"    long_name="ocean current along i-axis squared (thickness weighted)"                                            unit="m3/s2"      grid_ref="grid_U_3D"  > uoce * uoce * e3u </field>
         <field id="ssu"          long_name="ocean surface current along i-axis"                                                                 unit="m/s"                             />
         <field id="sbu"          long_name="ocean bottom current along i-axis"                                                                  unit="m/s"                             />
@@ -410,5 +405,4 @@
         <field id="voce"         long_name="ocean current along j-axis"                             standard_name="sea_water_y_velocity"        unit="m/s"        grid_ref="grid_V_3D" />
         <field id="voce_e3v"     long_name="ocean current along j-axis  (thickness weighted)"                                                   unit="m/s"        grid_ref="grid_V_3D"  > voce * e3v </field>
-        <field id="voce2_e3v"    long_name="ocean current along j-axis squared (thickness weighted)"                                            unit="m3/s2"      grid_ref="grid_V_3D"  > voce * voce * e3v </field>
         <field id="ssv"          long_name="ocean surface current along j-axis"                                                                 unit="m/s"                             />
         <field id="sbv"          long_name="ocean bottom current along j-axis"                                                                  unit="m/s"                             />
@@ -981,5 +975,5 @@
 
     <field_group id="25h_grid_W" grid_ref="grid_W_3D" operation="instant">
-      <field id="vomecrtz25h"         name="k current 25h mean"                 unit="m/s"      />
+      <field id="vovecrtz25h"         name="k current 25h mean"                 unit="m/s"      />
       <field id="avt25h"              name="vertical diffusivity25h mean"       unit="m2/s" />
       <field id="avm25h"              name="vertical viscosity 25h mean"        unit="m2/s" />
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/C1D/step_c1d.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/C1D/step_c1d.F90	(revision 11056)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/C1D/step_c1d.F90	(revision 11057)
@@ -107,8 +107,8 @@
       IF(.NOT.ln_linssh)CALL tra_adv( kstp, Nbb, Nnn, ts, Nrhs  )  ! horizontal & vertical advection
       IF( ln_zdfosm  )  CALL tra_osm( kstp, Nnn     , ts, Nrhs  )  ! OSMOSIS non-local tracer fluxes
-                        CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa   ) ! vertical mixing
+                        CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa   )         ! vertical mixing
                         CALL eos( ts(:,:,:,:,Nnn), rhd, rhop, gdept_0(:,:,:) )  ! now potential density for zdfmxl
-      IF( ln_zdfnpc )   CALL tra_npc( kstp,      Nnn, Nrhs, ts, Naa   ) ! applied non penetrative convective adjustment on (t,s)
-                        CALL tra_nxt( kstp, Nbb, Nnn, Nrhs,     Naa   ) ! tracer fields at next time step
+      IF( ln_zdfnpc )   CALL tra_npc( kstp,      Nnn, Nrhs, ts, Naa   )         ! applied non penetrative convective adjustment on (t,s)
+                        CALL tra_atf( kstp, Nbb, Nnn, Nrhs,     Naa, ts   )     ! time filtering of "now" tracer fields
 
 
@@ -124,8 +124,14 @@
       IF( ln_zdfosm  )  CALL dyn_osm    ( kstp,      Nnn      , uu, vv, Nrhs )  ! OSMOSIS non-local velocity fluxes
                         CALL dyn_zdf    ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion
-                        CALL dyn_nxt    ( kstp, Nbb, Nnn, Naa )                 ! lateral velocity at next time step
-      IF(.NOT.ln_linssh)CALL ssh_swp    ( kstp, Nbb, Nnn, Naa )                 ! swap of sea surface height
-
-      IF(.NOT.ln_linssh)CALL dom_vvl_sf_swp( kstp, Nbb, Nnn, Naa )              ! swap of vertical scale factors
+                        CALL dyn_atf    ( kstp, Nbb, Nnn, Naa , uu, vv, e3t, e3u, e3v )  ! time filtering of "now" fields
+      IF(.NOT.ln_linssh)CALL ssh_atf    ( kstp, Nbb, Nnn, Naa , ssh )                    ! time filtering of "now" sea surface height
+      !
+      ! Swap time levels
+      Nrhs = Nbb
+      Nbb = Nnn
+      Nnn = Naa
+      Naa = Nrhs
+      !
+      IF(.NOT.ln_linssh)CALL dom_vvl_sf_swp( kstp, Nbb, Nnn, Naa )                       ! swap of vertical scale factors
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DYN/sshwzv.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DYN/sshwzv.F90	(revision 11056)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DYN/sshwzv.F90	(revision 11057)
@@ -13,5 +13,5 @@
    !!----------------------------------------------------------------------
    !!   ssh_nxt       : after ssh
-   !!   ssh_swp       : filter ans swap the ssh arrays
+   !!   ssh_atf       : filter ans swap the ssh arrays
    !!   wzv           : compute now vertical velocity
    !!----------------------------------------------------------------------
@@ -43,5 +43,5 @@
    PUBLIC   wzv        ! called by step.F90
    PUBLIC   wAimp      ! called by step.F90
-   PUBLIC   ssh_swp    ! called by step.F90
+   PUBLIC   ssh_atf    ! called by step.F90
 
    !! * Substitutions
@@ -218,7 +218,9 @@
 
 
-   SUBROUTINE ssh_swp( kt, Kbb, Kmm, Kaa )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE ssh_nxt  ***
+   SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE ssh_atf  ***
+      !!
+      !!      !!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!
       !!
       !! ** Purpose :   achieve the sea surface  height time stepping by 
@@ -237,38 +239,35 @@
       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt              ! ocean time-step index
-      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa   ! ocean time-step index
+      INTEGER                         , INTENT(in   ) ::   kt             ! ocean time-step index
+      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! SSH field
       !
       REAL(wp) ::   zcoef   ! local scalar
       !!----------------------------------------------------------------------
       !
-      IF( ln_timing )   CALL timing_start('ssh_swp')
+      IF( ln_timing )   CALL timing_start('ssh_atf')
       !
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
+         IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height'
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
       ENDIF
       !              !==  Euler time-stepping: no filter, just swap  ==!
-      IF ( neuler == 0 .AND. kt == nit000 ) THEN
-         ssh(:,:,Kmm) = ssh(:,:,Kaa)                              ! now    <-- after  (before already = now)
-         !
-      ELSE           !==  Leap-Frog time-stepping: Asselin filter + swap  ==!
-         !                                                  ! before <-- now filtered
-         ssh(:,:,Kbb) = ssh(:,:,Kmm) + atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa) )
-         IF( .NOT.ln_linssh ) THEN                          ! before <-- with forcing removed
+      IF ( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN   ! Only do time filtering for leapfrog timesteps
+         !                                                  ! filtered "now" field
+         ssh(:,:,Kmm) = ssh(:,:,Kmm) + atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa) )
+         IF( .NOT.ln_linssh ) THEN                          ! "now" <-- with forcing removed
             zcoef = atfp * rdt * r1_rau0
-            ssh(:,:,Kbb) = ssh(:,:,Kbb) - zcoef * (     emp_b(:,:) - emp   (:,:)   &
+            ssh(:,:,Kmm) = ssh(:,:,Kmm) - zcoef * (     emp_b(:,:) - emp   (:,:)   &
                &                             -    rnf_b(:,:) + rnf   (:,:)   &
                &                             + fwfisf_b(:,:) - fwfisf(:,:)   ) * ssmask(:,:)
          ENDIF
-         ssh(:,:,Kmm) = ssh(:,:,Kaa)                              ! now <-- after
-      ENDIF
-      !
-      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssh(:,:,Kbb), clinfo1=' ssh(:,:,Kbb)  - : ', mask1=tmask )
-      !
-      IF( ln_timing )   CALL timing_stop('ssh_swp')
-      !
-   END SUBROUTINE ssh_swp
+      ENDIF
+      !
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssh(:,:,Kmm), clinfo1=' ssh(:,:,Kmm)  - : ', mask1=tmask )
+      !
+      IF( ln_timing )   CALL timing_stop('ssh_atf')
+      !
+   END SUBROUTINE ssh_atf
 
    SUBROUTINE wAimp( kt, Kmm )
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/TRA/traatf.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/TRA/traatf.F90	(revision 11057)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/TRA/traatf.F90	(revision 11057)
@@ -0,0 +1,382 @@
+MODULE traatf
+   !!======================================================================
+   !!                       ***  MODULE  traatf  ***
+   !! Ocean active tracers:  time stepping on temperature and salinity
+   !!======================================================================
+   !! History :  OPA  !  1991-11  (G. Madec)  Original code
+   !!            7.0  !  1993-03  (M. Guyon)  symetrical conditions
+   !!            8.0  !  1996-02  (G. Madec & M. Imbard)  opa release 8.0
+   !!             -   !  1996-04  (A. Weaver)  Euler forward step
+   !!            8.2  !  1999-02  (G. Madec, N. Grima)  semi-implicit pressure grad.
+   !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module
+   !!             -   !  2002-11  (C. Talandier, A-M Treguier) Open boundaries
+   !!             -   !  2005-04  (C. Deltel) Add Asselin trend in the ML budget
+   !!            2.0  !  2006-02  (L. Debreu, C. Mazauric) Agrif implementation
+   !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf
+   !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  semi-implicit hpg with asselin filter + modified LF-RA
+   !!             -   !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   tra_atf       : time stepping on tracers
+   !!   tra_atf_fix   : time stepping on tracers : fixed    volume case
+   !!   tra_atf_vvl   : time stepping on tracers : variable volume case
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain variables 
+   USE sbc_oce         ! surface boundary condition: ocean
+   USE sbcrnf          ! river runoffs
+   USE sbcisf          ! ice shelf melting
+   USE zdf_oce         ! ocean vertical mixing
+   USE domvvl          ! variable volume
+   USE trd_oce         ! trends: ocean variables
+   USE trdtra          ! trends manager: tracers 
+   USE traqsr          ! penetrative solar radiation (needed for nksr)
+   USE phycst          ! physical constant
+   USE ldftra          ! lateral physics : tracers
+   USE ldfslp          ! lateral physics : slopes
+   USE bdy_oce  , ONLY : ln_bdy
+   USE bdytra          ! open boundary condition (bdy_tra routine)
+   !
+   USE in_out_manager  ! I/O manager
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl          ! Print control
+   USE timing          ! Timing
+#if defined key_agrif
+   USE agrif_oce_interp
+#endif
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   tra_atf       ! routine called by step.F90
+   PUBLIC   tra_atf_fix   ! to be used in trcnxt
+   PUBLIC   tra_atf_vvl   ! to be used in trcnxt
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE tra_atf( kt, Kbb, Kmm, Kaa, pts )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE traatf  ***
+      !!
+      !!       !!!!!!!!!!!!!!!!!  REWRITE HEADER COMMENTS !!!!!!!!!!!!!!!
+      !!
+      !! ** Purpose :   Apply the boundary condition on the after temperature  
+      !!             and salinity fields, achieved the time stepping by adding
+      !!             the Asselin filter on now fields and swapping the fields.
+      !! 
+      !! ** Method  :   At this stage of the computation, ta and sa are the 
+      !!             after temperature and salinity as the time stepping has
+      !!             been performed in trazdf_imp or trazdf_exp module.
+      !!
+      !!              - Apply lateral boundary conditions on (ta,sa) 
+      !!             at the local domain   boundaries through lbc_lnk call, 
+      !!             at the one-way open boundaries (ln_bdy=T), 
+      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
+      !!
+      !!              - Update lateral boundary conditions on AGRIF children
+      !!             domains (lk_agrif=T)
+      !!
+      !! ** Action  : - ts(Kbb) & ts(Kmm) ready for the next time step
+      !!----------------------------------------------------------------------
+      INTEGER                                  , INTENT(in   ) :: kt             ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) :: Kbb, Kmm, Kaa  ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts            ! active tracers 
+      !!
+      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
+      REAL(wp) ::   zfact            ! local scalars
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start( 'tra_atf')
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_atf : achieve the time stepping by Asselin filter and array swap'
+         IF(lwp) WRITE(numout,*) '~~~~~~~'
+      ENDIF
+
+      ! Update after tracer on domain lateral boundaries
+      ! 
+#if defined key_agrif
+      CALL Agrif_tra                     ! AGRIF zoom boundaries
+#endif
+      !                                              ! local domain boundaries  (T-point, unchanged sign)
+      CALL lbc_lnk_multi( 'traatf', pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. )
+      !
+      IF( ln_bdy )   CALL bdy_tra( kt, Kbb, pts, Kaa )  ! BDY open boundaries
+ 
+      ! set time step size (Euler/Leapfrog)
+      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =        rdt   ! at nit000             (Euler)
+      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt = 2._wp* rdt   ! at nit000 or nit000+1 (Leapfrog)
+      ENDIF
+
+      ! trends computation initialisation
+      IF( l_trdtra )   THEN                    
+         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
+         ztrdt(:,:,jpk) = 0._wp
+         ztrds(:,:,jpk) = 0._wp
+         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend 
+            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_zdfp, ztrdt )
+            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_zdfp, ztrds )
+         ENDIF
+         ! total trend for the non-time-filtered variables. 
+         zfact = 1.0 / rdt
+         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from pts(Kmm) terms
+         DO jk = 1, jpkm1
+            ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - pts(:,:,jk,jp_tem,Kmm)) * zfact
+            ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - pts(:,:,jk,jp_sal,Kmm)) * zfact
+         END DO
+         CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_tot, ztrdt )
+         CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_tot, ztrds )
+         IF( ln_linssh ) THEN       ! linear sea surface height only
+            ! Store now fields before applying the Asselin filter 
+            ! in order to calculate Asselin filter trend later.
+            ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kmm) 
+            ztrds(:,:,:) = pts(:,:,:,jp_sal,Kmm)
+         ENDIF
+      ENDIF
+
+      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
+         DO jn = 1, jpts
+            DO jk = 1, jpkm1
+               pts(:,:,jk,jn,Kmm) = pts(:,:,jk,jn,Kaa)    
+            END DO
+         END DO
+         IF (l_trdtra .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl
+            !                                        ! Asselin filter is output by tra_atf_vvl that is not called on this time step
+            ztrdt(:,:,:) = 0._wp
+            ztrds(:,:,:) = 0._wp
+            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_atf, ztrdt )
+            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_atf, ztrds )
+         END IF
+         !
+      ELSE                                            ! Leap-Frog + Asselin filter time stepping
+         !
+         IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000,      'TRA', pts, jpts )  ! linear free surface 
+         ELSE                   ;   CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nit000, rdt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface
+         ENDIF
+         !
+         CALL lbc_lnk_multi( 'traatf', pts(:,:,:,jp_tem,Kbb) , 'T', 1., pts(:,:,:,jp_sal,Kbb) , 'T', 1., &
+                  &                    pts(:,:,:,jp_tem,Kmm) , 'T', 1., pts(:,:,:,jp_sal,Kmm) , 'T', 1., &
+                  &                    pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1.  )
+         !
+      ENDIF     
+      !
+      IF( l_trdtra .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
+         zfact = 1._wp / r2dt             
+         DO jk = 1, jpkm1
+            ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kbb) - ztrdt(:,:,jk) ) * zfact
+            ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kbb) - ztrds(:,:,jk) ) * zfact
+         END DO
+         CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_atf, ztrdt )
+         CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_atf, ztrds )
+      END IF
+      IF( l_trdtra )   DEALLOCATE( ztrdt , ztrds )
+      !
+      !                        ! control print
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kmm), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
+         &                       tab3d_2=pts(:,:,:,jp_sal,Kmm), clinfo2=       ' Sn: ', mask2=tmask )
+      !
+      IF( ln_timing )   CALL timing_stop('tra_atf')
+      !
+   END SUBROUTINE tra_atf
+
+
+   SUBROUTINE tra_atf_fix( kt, kit000, Kbb, Kmm, Kaa, cdtype, pt, kjpt )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE tra_atf_fix  ***
+      !!
+      !! ** Purpose :   fixed volume: apply the Asselin time filter and 
+      !!                swap the tracer fields.
+      !! 
+      !! ** Method  : - Apply a Asselin time filter on now fields.
+      !!              - swap tracer fields to prepare the next time_step.
+      !!
+      !! ** Action  : - ptb & ptn ready for the next time step
+      !!----------------------------------------------------------------------
+      INTEGER                                  , INTENT(in   ) ::  kt            ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) ::  Kbb, Kmm, Kaa ! time level indices
+      INTEGER                                  , INTENT(in   ) ::  kit000        ! first time step index
+      CHARACTER(len=3)                         , INTENT(in   ) ::  cdtype        ! =TRA or TRC (tracer indicator)
+      INTEGER                                  , INTENT(in   ) ::  kjpt          ! number of tracers
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::  pt            ! tracer fields
+      !
+      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
+      REAL(wp) ::   ztn, ztd         ! local scalars
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == kit000 )  THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_atf_fix : time stepping', cdtype
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      ENDIF
+      !
+      DO jn = 1, kjpt
+         !
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  ztn = pt(ji,jj,jk,jn,Kmm)                                    
+                  ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb)  ! time laplacian on tracers
+                  !
+                  pt(ji,jj,jk,jn,Kmm) = ztn + atfp * ztd                      ! ptb <-- filtered ptn 
+               END DO
+           END DO
+         END DO
+         !
+      END DO
+      !
+   END SUBROUTINE tra_atf_fix
+
+
+   SUBROUTINE tra_atf_vvl( kt, Kbb, Kmm, Kaa, kit000, p2dt, cdtype, pt, psbc_tc, psbc_tc_b, kjpt )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE tra_atf_vvl  ***
+      !!
+      !! ** Purpose :   Time varying volume: apply the Asselin time filter  
+      !!                and swap the tracer fields.
+      !! 
+      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
+      !!              - swap tracer fields to prepare the next time_step.
+      !!             tb  = ( e3t(Kmm)*tn + atfp*[ e3t(Kbb)*tb - 2 e3t(Kmm)*tn + e3t_a*ta ] )
+      !!                  /( e3t(Kmm)    + atfp*[ e3t(Kbb)    - 2 e3t(Kmm)    + e3t(Kaa)    ] )
+      !!             tn  = ta 
+      !!
+      !! ** Action  : - ptb & ptn ready for the next time step
+      !!----------------------------------------------------------------------
+      INTEGER                                  , INTENT(in   ) ::  kt        ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) ::  Kbb, Kmm, Kaa ! time level indices
+      INTEGER                                  , INTENT(in   ) ::  kit000    ! first time step index
+      REAL(wp)                                 , INTENT(in   ) ::  p2dt      ! time-step
+      CHARACTER(len=3)                         , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator)
+      INTEGER                                  , INTENT(in   ) ::  kjpt      ! number of tracers
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::  pt        ! tracer fields
+      REAL(wp), DIMENSION(jpi,jpj    ,kjpt)    , INTENT(in   ) ::  psbc_tc   ! surface tracer content
+      REAL(wp), DIMENSION(jpi,jpj    ,kjpt)    , INTENT(in   ) ::  psbc_tc_b ! before surface tracer content
+      !
+      LOGICAL  ::   ll_traqsr, ll_rnf, ll_isf   ! local logical
+      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
+      REAL(wp) ::   zfact, zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
+      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrd_atf
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == kit000 )  THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_atf_vvl : time stepping', cdtype
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      ENDIF
+      !
+      IF( cdtype == 'TRA' )  THEN   
+         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
+         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
+         ll_isf     = ln_isf           ! active  tracers case  and  ice shelf melting
+      ELSE                          ! passive tracers case
+         ll_traqsr  = .FALSE.          ! NO solar penetration
+         ll_rnf     = .FALSE.          ! NO river runoffs ????          !!gm BUG ?  
+         ll_isf     = .FALSE.          ! NO ice shelf melting/freezing  !!gm BUG ?? 
+      ENDIF
+      !
+      IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) )   THEN
+         ALLOCATE( ztrd_atf(jpi,jpj,jpk,kjpt) )
+         ztrd_atf(:,:,:,:) = 0.0_wp
+      ENDIF
+      zfact = 1._wp / p2dt
+      zfact1 = atfp * p2dt
+      zfact2 = zfact1 * r1_rau0
+      DO jn = 1, kjpt      
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  ze3t_b = e3t(ji,jj,jk,Kbb)
+                  ze3t_n = e3t(ji,jj,jk,Kmm)
+                  ze3t_a = e3t(ji,jj,jk,Kaa)
+                  !                                         ! tracer content at Before, now and after
+                  ztc_b  = pt(ji,jj,jk,jn,Kbb) * ze3t_b
+                  ztc_n  = pt(ji,jj,jk,jn,Kmm) * ze3t_n
+                  ztc_a  = pt(ji,jj,jk,jn,Kaa) * ze3t_a
+                  !
+                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
+                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
+                  !
+                  ze3t_f = ze3t_n + atfp * ze3t_d
+                  ztc_f  = ztc_n  + atfp * ztc_d
+                  !
+                  IF( jk == mikt(ji,jj) ) THEN           ! first level 
+                     ze3t_f = ze3t_f - zfact2 * ( (emp_b(ji,jj)    - emp(ji,jj)   )  &
+                            &                   + (fwfisf_b(ji,jj) - fwfisf(ji,jj))  )
+                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
+                  ENDIF
+                  IF( ln_rnf_depth ) THEN
+                     ! Rivers are not just at the surface must go down to nk_rnf(ji,jj)
+                     IF( mikt(ji,jj) <=jk .and. jk <= nk_rnf(ji,jj)  ) THEN
+                        ze3t_f = ze3t_f - zfact2 * ( - (rnf_b(ji,jj) - rnf(ji,jj)   )  ) &
+                    &                            * ( e3t(ji,jj,jk,Kmm) / h_rnf(ji,jj) ) 
+                     ENDIF
+                  ELSE
+                     IF( jk == mikt(ji,jj) ) THEN           ! first level 
+                        ze3t_f = ze3t_f - zfact2 * ( - (rnf_b(ji,jj)    - rnf(ji,jj)   ) ) 
+                     ENDIF
+                  ENDIF
+
+                  !
+                  ! solar penetration (temperature only)
+                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            & 
+                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
+                     !
+                  ! river runoff
+                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          &
+                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
+                     &                              * e3t(ji,jj,jk,Kmm) / h_rnf(ji,jj)
+                     !
+                  ! ice shelf
+                  IF( ll_isf ) THEN
+                     ! level fully include in the Losch_2008 ice shelf boundary layer
+                     IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) )                          &
+                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
+                               &                 * e3t(ji,jj,jk,Kmm) * r1_hisf_tbl (ji,jj)
+                     ! level partially include in Losch_2008 ice shelf boundary layer 
+                     IF ( jk == misfkb(ji,jj) )                                                   &
+                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
+                               &                 * e3t(ji,jj,jk,Kmm) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj)
+                  END IF
+                  !
+                  ze3t_f = 1.e0 / ze3t_f
+                  pt(ji,jj,jk,jn,Kmm) = ztc_f * ze3t_f    ! time filtered "now" field
+                  !
+                  IF( ( l_trdtra .and. cdtype == 'TRA' ) .OR. ( l_trdtrc .and. cdtype == 'TRC' ) ) THEN
+                     ztrd_atf(ji,jj,jk,jn) = (ztc_f - ztc_n) * zfact/ze3t_n
+                  ENDIF
+                  !
+               END DO
+            END DO
+         END DO
+         ! 
+      END DO
+      !
+      IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) )   THEN
+         IF( l_trdtra .AND. cdtype == 'TRA' ) THEN 
+            CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) )
+            CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) )
+         ENDIF
+         IF( l_trdtrc .AND. cdtype == 'TRC' ) THEN
+            DO jn = 1, kjpt
+               CALL trd_tra( kt, Kmm, Kaa, cdtype, jn, jptra_atf, ztrd_atf(:,:,:,jn) )
+            END DO
+         ENDIF
+         DEALLOCATE( ztrd_atf )
+      ENDIF
+      !
+   END SUBROUTINE tra_atf_vvl
+
+   !!======================================================================
+END MODULE traatf
Index: MO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/TRA/tranxt.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/TRA/tranxt.F90	(revision 11056)
+++ 	(revision )
@@ -1,392 +1,0 @@
-MODULE tranxt
-   !!======================================================================
-   !!                       ***  MODULE  tranxt  ***
-   !! Ocean active tracers:  time stepping on temperature and salinity
-   !!======================================================================
-   !! History :  OPA  !  1991-11  (G. Madec)  Original code
-   !!            7.0  !  1993-03  (M. Guyon)  symetrical conditions
-   !!            8.0  !  1996-02  (G. Madec & M. Imbard)  opa release 8.0
-   !!             -   !  1996-04  (A. Weaver)  Euler forward step
-   !!            8.2  !  1999-02  (G. Madec, N. Grima)  semi-implicit pressure grad.
-   !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module
-   !!             -   !  2002-11  (C. Talandier, A-M Treguier) Open boundaries
-   !!             -   !  2005-04  (C. Deltel) Add Asselin trend in the ML budget
-   !!            2.0  !  2006-02  (L. Debreu, C. Mazauric) Agrif implementation
-   !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf
-   !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option
-   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  semi-implicit hpg with asselin filter + modified LF-RA
-   !!             -   !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   tra_nxt       : time stepping on tracers
-   !!   tra_nxt_fix   : time stepping on tracers : fixed    volume case
-   !!   tra_nxt_vvl   : time stepping on tracers : variable volume case
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain variables 
-   USE sbc_oce         ! surface boundary condition: ocean
-   USE sbcrnf          ! river runoffs
-   USE sbcisf          ! ice shelf melting
-   USE zdf_oce         ! ocean vertical mixing
-   USE domvvl          ! variable volume
-   USE trd_oce         ! trends: ocean variables
-   USE trdtra          ! trends manager: tracers 
-   USE traqsr          ! penetrative solar radiation (needed for nksr)
-   USE phycst          ! physical constant
-   USE ldftra          ! lateral physics : tracers
-   USE ldfslp          ! lateral physics : slopes
-   USE bdy_oce  , ONLY : ln_bdy
-   USE bdytra          ! open boundary condition (bdy_tra routine)
-   !
-   USE in_out_manager  ! I/O manager
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE prtctl          ! Print control
-   USE timing          ! Timing
-#if defined key_agrif
-   USE agrif_oce_interp
-#endif
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   tra_nxt       ! routine called by step.F90
-   PUBLIC   tra_nxt_fix   ! to be used in trcnxt
-   PUBLIC   tra_nxt_vvl   ! to be used in trcnxt
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE tra_nxt( kt, Kbb, Kmm, Krhs, Kaa )
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE tranxt  ***
-      !!
-      !! ** Purpose :   Apply the boundary condition on the after temperature  
-      !!             and salinity fields, achieved the time stepping by adding
-      !!             the Asselin filter on now fields and swapping the fields.
-      !! 
-      !! ** Method  :   At this stage of the computation, ta and sa are the 
-      !!             after temperature and salinity as the time stepping has
-      !!             been performed in trazdf_imp or trazdf_exp module.
-      !!
-      !!              - Apply lateral boundary conditions on (ta,sa) 
-      !!             at the local domain   boundaries through lbc_lnk call, 
-      !!             at the one-way open boundaries (ln_bdy=T), 
-      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
-      !!
-      !!              - Update lateral boundary conditions on AGRIF children
-      !!             domains (lk_agrif=T)
-      !!
-      !! ** Action  : - ts(Kbb) & ts(Kmm) ready for the next time step
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt                   ! ocean time-step index
-      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs, Kaa  ! time level indices
-      !!
-      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp) ::   zfact            ! local scalars
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
-      !!----------------------------------------------------------------------
-      !
-      IF( ln_timing )   CALL timing_start( 'tra_nxt')
-      !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
-         IF(lwp) WRITE(numout,*) '~~~~~~~'
-      ENDIF
-
-      ! Update after tracer on domain lateral boundaries
-      ! 
-#if defined key_agrif
-      CALL Agrif_tra                     ! AGRIF zoom boundaries
-#endif
-      !                                              ! local domain boundaries  (T-point, unchanged sign)
-      CALL lbc_lnk_multi( 'tranxt', ts(:,:,:,jp_tem,Krhs), 'T', 1., ts(:,:,:,jp_sal,Krhs), 'T', 1. )
-      !
-      !! IMMERSE development : Following the general pattern for the new code we want to pass in the 
-      !!                       velocities to bdy_dyn as arguments so here we use "ts" even 
-      !!                       though we haven't converted the tracer names in the rest of tranxt.F90
-      !!                       because it will be completely rewritten. DS.
-      IF( ln_bdy )   CALL bdy_tra( kt, Kbb, ts, Kaa )  ! BDY open boundaries
- 
-      ! set time step size (Euler/Leapfrog)
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =        rdt   ! at nit000             (Euler)
-      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt = 2._wp* rdt   ! at nit000 or nit000+1 (Leapfrog)
-      ENDIF
-
-      ! trends computation initialisation
-      IF( l_trdtra )   THEN                    
-         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
-         ztrdt(:,:,jpk) = 0._wp
-         ztrds(:,:,jpk) = 0._wp
-         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend 
-            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_zdfp, ztrdt )
-            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_zdfp, ztrds )
-         ENDIF
-         ! total trend for the non-time-filtered variables. 
-         zfact = 1.0 / rdt
-         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from ts(Kmm) terms
-         DO jk = 1, jpkm1
-            ztrdt(:,:,jk) = ( ts(:,:,jk,jp_tem,Krhs)*e3t(:,:,jk,Krhs) / e3t(:,:,jk,Kmm) - ts(:,:,jk,jp_tem,Kmm)) * zfact
-            ztrds(:,:,jk) = ( ts(:,:,jk,jp_sal,Krhs)*e3t(:,:,jk,Krhs) / e3t(:,:,jk,Kmm) - ts(:,:,jk,jp_sal,Kmm)) * zfact
-         END DO
-         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_tot, ztrdt )
-         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_tot, ztrds )
-         IF( ln_linssh ) THEN       ! linear sea surface height only
-            ! Store now fields before applying the Asselin filter 
-            ! in order to calculate Asselin filter trend later.
-            ztrdt(:,:,:) = ts(:,:,:,jp_tem,Kmm) 
-            ztrds(:,:,:) = ts(:,:,:,jp_sal,Kmm)
-         ENDIF
-      ENDIF
-
-      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
-         DO jn = 1, jpts
-            DO jk = 1, jpkm1
-               ts(:,:,jk,jn,Kmm) = ts(:,:,jk,jn,Krhs)    
-            END DO
-         END DO
-         IF (l_trdtra .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl
-            !                                        ! Asselin filter is output by tra_nxt_vvl that is not called on this time step
-            ztrdt(:,:,:) = 0._wp
-            ztrds(:,:,:) = 0._wp
-            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_atf, ztrdt )
-            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_atf, ztrds )
-         END IF
-         !
-      ELSE                                            ! Leap-Frog + Asselin filter time stepping
-         !
-         IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt,      Kmm,       nit000,      'TRA',                           &
-           &                                                   ts(:,:,:,:,Kbb), ts(:,:,:,:,Kmm), ts(:,:,:,:,Krhs), jpts )  ! linear free surface 
-         ELSE                   ;   CALL tra_nxt_vvl( kt, Kbb, Kmm, Krhs, nit000, rdt, 'TRA',                           &
-           &                                                   ts(:,:,:,:,Kbb), ts(:,:,:,:,Kmm), ts(:,:,:,:,Krhs),      &
-           &                                                   sbc_tsc        , sbc_tsc_b                        , jpts )  ! non-linear free surface
-         ENDIF
-         !
-         CALL lbc_lnk_multi( 'tranxt', ts(:,:,:,jp_tem,Kbb) , 'T', 1., ts(:,:,:,jp_sal,Kbb) , 'T', 1., &
-                  &                    ts(:,:,:,jp_tem,Kmm) , 'T', 1., ts(:,:,:,jp_sal,Kmm) , 'T', 1., &
-                  &                    ts(:,:,:,jp_tem,Krhs), 'T', 1., ts(:,:,:,jp_sal,Krhs), 'T', 1.  )
-         !
-      ENDIF     
-      !
-      IF( l_trdtra .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
-         zfact = 1._wp / r2dt             
-         DO jk = 1, jpkm1
-            ztrdt(:,:,jk) = ( ts(:,:,jk,jp_tem,Kbb) - ztrdt(:,:,jk) ) * zfact
-            ztrds(:,:,jk) = ( ts(:,:,jk,jp_sal,Kbb) - ztrds(:,:,jk) ) * zfact
-         END DO
-         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_atf, ztrdt )
-         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_atf, ztrds )
-      END IF
-      IF( l_trdtra )   DEALLOCATE( ztrdt , ztrds )
-      !
-      !                        ! control print
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
-         &                       tab3d_2=ts(:,:,:,jp_sal,Kmm), clinfo2=       ' Sn: ', mask2=tmask )
-      !
-      IF( ln_timing )   CALL timing_stop('tra_nxt')
-      !
-   END SUBROUTINE tra_nxt
-
-
-   SUBROUTINE tra_nxt_fix( kt, Kmm, kit000, cdtype, ptb, ptn, pta, kjpt )
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE tra_nxt_fix  ***
-      !!
-      !! ** Purpose :   fixed volume: apply the Asselin time filter and 
-      !!                swap the tracer fields.
-      !! 
-      !! ** Method  : - Apply a Asselin time filter on now fields.
-      !!              - swap tracer fields to prepare the next time_step.
-      !!
-      !! ** Action  : - ptb & ptn ready for the next time step
-      !!----------------------------------------------------------------------
-      INTEGER                              , INTENT(in   ) ::  kt        ! ocean time-step index
-      INTEGER                              , INTENT(in   ) ::  Kmm       ! time level index
-      INTEGER                              , INTENT(in   ) ::  kit000    ! first time step index
-      CHARACTER(len=3)                     , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator)
-      INTEGER                              , INTENT(in   ) ::  kjpt      ! number of tracers
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptb       ! before tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptn       ! now tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  pta       ! tracer trend
-      !
-      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp) ::   ztn, ztd         ! local scalars
-      !!----------------------------------------------------------------------
-      !
-      IF( kt == kit000 )  THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping', cdtype
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
-      !
-      DO jn = 1, kjpt
-         !
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  ztn = ptn(ji,jj,jk,jn)                                    
-                  ztd = pta(ji,jj,jk,jn) - 2._wp * ztn + ptb(ji,jj,jk,jn)  ! time laplacian on tracers
-                  !
-                  ptb(ji,jj,jk,jn) = ztn + atfp * ztd                      ! ptb <-- filtered ptn 
-                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                      ! ptn <-- pta
-               END DO
-           END DO
-         END DO
-         !
-      END DO
-      !
-   END SUBROUTINE tra_nxt_fix
-
-
-   SUBROUTINE tra_nxt_vvl( kt, Kbb, Kmm, Krhs, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt )
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE tra_nxt_vvl  ***
-      !!
-      !! ** Purpose :   Time varying volume: apply the Asselin time filter  
-      !!                and swap the tracer fields.
-      !! 
-      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
-      !!              - swap tracer fields to prepare the next time_step.
-      !!             tb  = ( e3t(Kmm)*tn + atfp*[ e3t(Kbb)*tb - 2 e3t(Kmm)*tn + e3t_a*ta ] )
-      !!                  /( e3t(Kmm)    + atfp*[ e3t(Kbb)    - 2 e3t(Kmm)    + e3t(Krhs)    ] )
-      !!             tn  = ta 
-      !!
-      !! ** Action  : - ptb & ptn ready for the next time step
-      !!----------------------------------------------------------------------
-      INTEGER                              , INTENT(in   ) ::  kt        ! ocean time-step index
-      INTEGER                              , INTENT(in   ) ::  Kbb, Kmm, Krhs ! time level indices
-      INTEGER                              , INTENT(in   ) ::  kit000    ! first time step index
-      REAL(wp)                             , INTENT(in   ) ::  p2dt      ! time-step
-      CHARACTER(len=3)                     , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator)
-      INTEGER                              , INTENT(in   ) ::  kjpt      ! number of tracers
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptb       ! before tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptn       ! now tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  pta       ! tracer trend
-      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::  psbc_tc   ! surface tracer content
-      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::  psbc_tc_b ! before surface tracer content
-      !
-      LOGICAL  ::   ll_traqsr, ll_rnf, ll_isf   ! local logical
-      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
-      REAL(wp) ::   zfact, zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
-      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrd_atf
-      !!----------------------------------------------------------------------
-      !
-      IF( kt == kit000 )  THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
-      !
-      IF( cdtype == 'TRA' )  THEN   
-         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
-         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
-         ll_isf     = ln_isf           ! active  tracers case  and  ice shelf melting
-      ELSE                          ! passive tracers case
-         ll_traqsr  = .FALSE.          ! NO solar penetration
-         ll_rnf     = .FALSE.          ! NO river runoffs ????          !!gm BUG ?  
-         ll_isf     = .FALSE.          ! NO ice shelf melting/freezing  !!gm BUG ?? 
-      ENDIF
-      !
-      IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) )   THEN
-         ALLOCATE( ztrd_atf(jpi,jpj,jpk,kjpt) )
-         ztrd_atf(:,:,:,:) = 0.0_wp
-      ENDIF
-      zfact = 1._wp / p2dt
-      zfact1 = atfp * p2dt
-      zfact2 = zfact1 * r1_rau0
-      DO jn = 1, kjpt      
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  ze3t_b = e3t(ji,jj,jk,Kbb)
-                  ze3t_n = e3t(ji,jj,jk,Kmm)
-                  ze3t_a = e3t(ji,jj,jk,Krhs)
-                  !                                         ! tracer content at Before, now and after
-                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
-                  ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
-                  ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
-                  !
-                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
-                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
-                  !
-                  ze3t_f = ze3t_n + atfp * ze3t_d
-                  ztc_f  = ztc_n  + atfp * ztc_d
-                  !
-                  IF( jk == mikt(ji,jj) ) THEN           ! first level 
-                     ze3t_f = ze3t_f - zfact2 * ( (emp_b(ji,jj)    - emp(ji,jj)   )  &
-                            &                   + (fwfisf_b(ji,jj) - fwfisf(ji,jj))  )
-                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
-                  ENDIF
-                  IF( ln_rnf_depth ) THEN
-                     ! Rivers are not just at the surface must go down to nk_rnf(ji,jj)
-                     IF( mikt(ji,jj) <=jk .and. jk <= nk_rnf(ji,jj)  ) THEN
-                        ze3t_f = ze3t_f - zfact2 * ( - (rnf_b(ji,jj) - rnf(ji,jj)   )  ) &
-                    &                            * ( e3t(ji,jj,jk,Kmm) / h_rnf(ji,jj) ) 
-                     ENDIF
-                  ELSE
-                     IF( jk == mikt(ji,jj) ) THEN           ! first level 
-                        ze3t_f = ze3t_f - zfact2 * ( - (rnf_b(ji,jj)    - rnf(ji,jj)   ) ) 
-                     ENDIF
-                  ENDIF
-
-                  !
-                  ! solar penetration (temperature only)
-                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            & 
-                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
-                     !
-                  ! river runoff
-                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          &
-                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
-                     &                              * e3t(ji,jj,jk,Kmm) / h_rnf(ji,jj)
-                     !
-                  ! ice shelf
-                  IF( ll_isf ) THEN
-                     ! level fully include in the Losch_2008 ice shelf boundary layer
-                     IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) )                          &
-                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
-                               &                 * e3t(ji,jj,jk,Kmm) * r1_hisf_tbl (ji,jj)
-                     ! level partially include in Losch_2008 ice shelf boundary layer 
-                     IF ( jk == misfkb(ji,jj) )                                                   &
-                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
-                               &                 * e3t(ji,jj,jk,Kmm) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj)
-                  END IF
-                  !
-                  ze3t_f = 1.e0 / ze3t_f
-                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
-                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
-                  !
-                  IF( ( l_trdtra .and. cdtype == 'TRA' ) .OR. ( l_trdtrc .and. cdtype == 'TRC' ) ) THEN
-                     ztrd_atf(ji,jj,jk,jn) = (ztc_f - ztc_n) * zfact/ze3t_n
-                  ENDIF
-                  !
-               END DO
-            END DO
-         END DO
-         ! 
-      END DO
-      !
-      IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) )   THEN
-         IF( l_trdtra .AND. cdtype == 'TRA' ) THEN 
-            CALL trd_tra( kt, Kmm, Krhs, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) )
-            CALL trd_tra( kt, Kmm, Krhs, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) )
-         ENDIF
-         IF( l_trdtrc .AND. cdtype == 'TRC' ) THEN
-            DO jn = 1, kjpt
-               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_atf, ztrd_atf(:,:,:,jn) )
-            END DO
-         ENDIF
-         DEALLOCATE( ztrd_atf )
-      ENDIF
-      !
-   END SUBROUTINE tra_nxt_vvl
-
-   !!======================================================================
-END MODULE tranxt
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/step.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/step.F90	(revision 11056)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/step.F90	(revision 11057)
@@ -280,7 +280,7 @@
 !! 
 !!jc2: dynnxt must be the latest call. e3t(:,:,:,Nbb) are indeed updated in that routine
-                         CALL tra_nxt       ( kstp, Nbb, Nnn, Nrhs, Naa )  ! finalize (bcs) tracer fields at next time step and swap
-                         CALL dyn_atf       ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v  )  ! time filtering of "now" arrays
-                         CALL ssh_swp       ( kstp, Nbb, Nnn, Naa )  ! swap of sea surface height
+                         CALL tra_atf       ( kstp, Nbb, Nnn, Naa, ts )                      ! time filtering of "now" tracer arrays
+                         CALL dyn_atf       ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v  )  ! time filtering of "now" velocities and scale factors
+                         CALL ssh_atf       ( kstp, Nbb, Nnn, Naa, ssh )                     ! time filtering of "now" sea surface height
       !
       ! Swap time levels
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/step_oce.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/step_oce.F90	(revision 11056)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/step_oce.F90	(revision 11057)
@@ -30,5 +30,5 @@
    USE traldf          ! lateral mixing                   (tra_ldf routine)
    USE trazdf          ! vertical mixing                  (tra_zdf routine)
-   USE tranxt          ! time-stepping                    (tra_nxt routine)
+   USE traatf          ! time filtering                   (tra_atf routine)
    USE tranpc          ! non-penetrative convection       (tra_npc routine)
 
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trcatf.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trcatf.F90	(revision 11057)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trcatf.F90	(revision 11057)
@@ -0,0 +1,271 @@
+MODULE trcatf
+   !!======================================================================
+   !!                       ***  MODULE  trcatf  ***
+   !! Ocean passive tracers:  time stepping on passives tracers
+   !!======================================================================
+   !! History :  7.0  !  1991-11  (G. Madec)  Original code
+   !!                 !  1993-03  (M. Guyon)  symetrical conditions
+   !!                 !  1995-02  (M. Levy)   passive tracers
+   !!                 !  1996-02  (G. Madec & M. Imbard)  opa release 8.0
+   !!            8.0  !  1996-04  (A. Weaver)  Euler forward step
+   !!            8.2  !  1999-02  (G. Madec, N. Grima)  semi-implicit pressure grad.
+   !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module
+   !!                 !  2002-08  (G. Madec)  F90: Free form and module
+   !!                 !  2002-11  (C. Talandier, A-M Treguier) Open boundaries
+   !!                 !  2004-03  (C. Ethe) passive tracers
+   !!                 !  2007-02  (C. Deltel) Diagnose ML trends for passive tracers
+   !!            2.0  !  2006-02  (L. Debreu, C. Mazauric) Agrif implementation
+   !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf
+   !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option
+   !!            3.3  !  2010-06  (C. Ethe, G. Madec) Merge TRA-TRC
+   !!----------------------------------------------------------------------
+#if defined key_top
+   !!----------------------------------------------------------------------
+   !!   'key_top'                                                TOP models
+   !!----------------------------------------------------------------------
+   !!   trc_atf     : time stepping on passive tracers
+   !!----------------------------------------------------------------------
+   USE oce_trc         ! ocean dynamics and tracers variables
+   USE trc             ! ocean passive tracers variables
+   USE trd_oce
+   USE trdtra
+   USE tranxt
+   USE bdy_oce   , ONLY: ln_bdy
+   USE trcbdy          ! BDY open boundaries
+# if defined key_agrif
+   USE agrif_top_interp
+# endif
+   !
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl_trc      ! Print control for debbuging
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   trc_atf   ! routine called by step.F90
+
+   REAL(wp) ::   rfact1, rfact2
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
+   !! $Id$ 
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE trc_atf( kt, Kbb, Kmm, Kaa, ptr )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE trcatf  ***
+      !!
+      !!        !!!!!!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!!!
+      !!
+      !! ** Purpose :   Compute the passive tracers fields at the 
+      !!      next time-step from their temporal trends and swap the fields.
+      !! 
+      !! ** Method  :   Apply lateral boundary conditions on (uu(Kaa),vv(Kaa)) through 
+      !!      call to lbc_lnk routine
+      !!   default:
+      !!      arrays swap
+      !!         (tr(Kmm)) = (tr(Kaa)) ; (tr(Kaa)) = (0,0)
+      !!         (tr(Kbb)) = (tr(Kmm)) 
+      !!
+      !!   For Arakawa or TVD Scheme : 
+      !!      A Asselin time filter applied on now tracers (tr(:,:,:,:,Kmm)) to avoid
+      !!      the divergence of two consecutive time-steps and tr arrays
+      !!      to prepare the next time_step:
+      !!         (tr(Kbb)) = (tr(Kmm)) + atfp [ (tr(Kbb)) + (tr(Kaa)) - 2 (tr(Kmm)) ]
+      !!         (tr(Kmm)) = (tr(Kaa)) ; (tr(Kaa)) = (0,0)
+      !!
+      !!
+      !! ** Action  : - update tr(Kbb), tr(Kmm)
+      !!----------------------------------------------------------------------
+      INTEGER                                   , INTENT( in )  :: kt             ! ocean time-step index
+      INTEGER                                   , INTENT( in )  :: Kbb, Kmm, Kaa ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr            ! passive tracers
+      !
+      INTEGER  ::   jk, jn   ! dummy loop indices
+      REAL(wp) ::   zfact            ! temporary scalar
+      CHARACTER (len=22) :: charout
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrdt    ! 4D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('trc_atf')
+      !
+      IF( kt == nittrc000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'trc_atf : Asselin time filtering on passive tracers'
+      ENDIF
+      !
+#if defined key_agrif
+      CALL Agrif_trc                   ! AGRIF zoom boundaries
+#endif
+      ! Update after tracer on domain lateral boundaries
+      CALL lbc_lnk( 'trcatf', ptr(:,:,:,:,Kaa), 'T', 1. )   
+
+      IF( ln_bdy )  CALL trc_bdy( kt, Kbb, Kmm, Kaa )
+
+      IF( l_trdtrc )  THEN             ! trends: store now fields before the Asselin filter application
+         ALLOCATE( ztrdt(jpi,jpj,jpk,jptra) )
+         ztrdt(:,:,:,:)  = 0._wp
+         IF( ln_traldf_iso ) THEN                       ! diagnose the "pure" Kz diffusive trend 
+            DO jn = 1, jptra
+               CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_zdfp, ztrdt(:,:,:,jn) )
+            ENDDO
+         ENDIF
+
+         ! total trend for the non-time-filtered variables. 
+         zfact = 1.0 / rdttrc
+         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from ts(Kmm) terms
+         IF( ln_linssh ) THEN       ! linear sea surface height only
+            DO jn = 1, jptra
+               DO jk = 1, jpkm1
+                  ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - ptr(:,:,jk,jn,Kmm)) * zfact
+               END DO
+            END DO
+         ELSE
+            DO jn = 1, jptra
+               DO jk = 1, jpkm1
+                  ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kaa) - ptr(:,:,jk,jn,Kmm) ) * zfact
+               END DO
+            END DO
+         ENDIF
+         !
+         DO jn = 1, jptra
+            CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_tot, ztrdt(:,:,:,jn) )
+         ENDDO
+         !
+         IF( ln_linssh ) THEN       ! linear sea surface height only
+            ! Store now fields before applying the Asselin filter 
+            ! in order to calculate Asselin filter trend later.
+            ztrdt(:,:,:,:) = ptr(:,:,:,:,Kmm) 
+         ENDIF
+
+      ENDIF
+      !                                ! Leap-Frog + Asselin filter time stepping
+      IF( .NOT.( (neuler == 0 .AND. kt == nittrc000) .OR. ln_top_euler ) ) THEN    ! Only time filter if not an Euler timestep
+         IF( .NOT. l_offline ) THEN ! Leap-Frog + Asselin filter time stepping
+            IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt, Kbb, Kmm, Kaa, nittrc000,         'TRC', ptr, jptra )                     !     linear ssh
+            ELSE                   ;   CALL tra_nxt_vvl( kt, Kbb, Kmm, Kaa, nittrc000, rdttrc, 'TRC', ptr, sbc_trc, sbc_trc_b, jptra ) ! non-linear ssh
+            ENDIF
+         ELSE
+                                       CALL trc_atf_off( kt, Kbb, Kmm, Kaa )       ! offline 
+         ENDIF
+         !
+         CALL lbc_lnk_multi( 'trcatf', ptr(:,:,:,:,Kbb), 'T', 1._wp, ptr(:,:,:,:,Kmm), 'T', 1._wp, ptr(:,:,:,:,Kaa), 'T', 1._wp )
+      ENDIF
+      !
+      IF( l_trdtrc .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt )
+         DO jn = 1, jptra
+            DO jk = 1, jpkm1
+               zfact = 1._wp / r2dttrc  
+               ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kbb) - ztrdt(:,:,jk,jn) ) * zfact 
+            END DO
+            CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
+         END DO
+      END IF
+      IF( l_trdtrc ) DEALLOCATE( ztrdt ) 
+      !
+      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+         WRITE(charout, FMT="('nxt')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc(tab4d=ptr(:,:,:,:,Kmm), mask=tmask, clinfo=ctrcnm)
+      ENDIF
+      !
+      IF( ln_timing )   CALL timing_stop('trc_atf')
+      !
+   END SUBROUTINE trc_atf
+
+
+   SUBROUTINE trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE tra_atf_off  ***
+      !!
+      !!          !!!!!!!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!!!
+      !!
+      !! ** Purpose :   Time varying volume: apply the Asselin time filter  
+      !! 
+      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
+      !!              - save in (ta,sa) a thickness weighted average over the three 
+      !!             time levels which will be used to compute rdn and thus the semi-
+      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
+      !!              - swap tracer fields to prepare the next time_step.
+      !!                This can be summurized for tempearture as:
+      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
+      !!                  /( e3t(:,:,:,Kmm)    + rbcp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Kaa)    ] )   
+      !!             ztm = 0                                                       otherwise
+      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
+      !!                  /( e3t(:,:,:,Kmm)    + atfp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Kaa)    ] )
+      !!             tn  = ta 
+      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
+      !!
+      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
+      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
+      !!----------------------------------------------------------------------
+      INTEGER                                   , INTENT(in   ) ::  kt            ! ocean time-step index
+      INTEGER                                   , INTENT(in   ) ::  Kbb, Kmm, Kaa ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) ::  ptr           ! passive tracers
+      !!     
+      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
+      REAL(wp) ::   ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
+      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == nittrc000 )  THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'trc_atf_off : Asselin time filtering'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+         IF( .NOT. ln_linssh ) THEN
+            rfact1 = atfp * rdttrc
+            rfact2 = rfact1 / rau0
+         ENDIF
+        !  
+      ENDIF
+      !
+      DO jn = 1, jptra      
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ze3t_b = e3t(ji,jj,jk,Kbb)
+                  ze3t_n = e3t(ji,jj,jk,Kmm)
+                  ze3t_a = e3t(ji,jj,jk,Kaa)
+                  !                                         ! tracer content at Before, now and after
+                  ztc_b  = ptr(ji,jj,jk,jn,Kbb)  * ze3t_b
+                  ztc_n  = ptr(ji,jj,jk,jn,Kmm)  * ze3t_n
+                  ztc_a  = ptr(ji,jj,jk,jn,Kaa) * ze3t_a
+                  !
+                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
+                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
+                  !
+                  ze3t_f = ze3t_n + atfp * ze3t_d
+                  ztc_f  = ztc_n  + atfp * ztc_d
+                  !
+                  IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN           ! first level 
+                     ze3t_f = ze3t_f - rfact2 * ( emp_b(ji,jj)      - emp(ji,jj)   ) 
+                     ztc_f  = ztc_f  - rfact1 * ( sbc_trc(ji,jj,jn) - sbc_trc_b(ji,jj,jn) )
+                  ENDIF
+
+                  ze3t_f = 1.e0 / ze3t_f
+                  ptr(ji,jj,jk,jn,Kmm) = ztc_f * ze3t_f     ! time filtered "now" field
+                  !
+               END DO
+            END DO
+         END DO
+         ! 
+      END DO
+      !
+   END SUBROUTINE trc_atf_off
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option                                         Empty module
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_atf( kt, Kbb, Kmm, Kaa, ptr )  
+      INTEGER                                   , INTENT(in) :: kt
+      INTEGER                                   , INTENT(in) :: Kbb, Kmm, Kaa  ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr            ! passive tracers
+      WRITE(*,*) 'trc_atf: You should not have seen this print! error?', kt
+   END SUBROUTINE trc_atf
+#endif
+   !!======================================================================
+END MODULE trcatf
Index: MO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trcnxt.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trcnxt.F90	(revision 11056)
+++ 	(revision )
@@ -1,284 +1,0 @@
-MODULE trcnxt
-   !!======================================================================
-   !!                       ***  MODULE  trcnxt  ***
-   !! Ocean passive tracers:  time stepping on passives tracers
-   !!======================================================================
-   !! History :  7.0  !  1991-11  (G. Madec)  Original code
-   !!                 !  1993-03  (M. Guyon)  symetrical conditions
-   !!                 !  1995-02  (M. Levy)   passive tracers
-   !!                 !  1996-02  (G. Madec & M. Imbard)  opa release 8.0
-   !!            8.0  !  1996-04  (A. Weaver)  Euler forward step
-   !!            8.2  !  1999-02  (G. Madec, N. Grima)  semi-implicit pressure grad.
-   !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module
-   !!                 !  2002-08  (G. Madec)  F90: Free form and module
-   !!                 !  2002-11  (C. Talandier, A-M Treguier) Open boundaries
-   !!                 !  2004-03  (C. Ethe) passive tracers
-   !!                 !  2007-02  (C. Deltel) Diagnose ML trends for passive tracers
-   !!            2.0  !  2006-02  (L. Debreu, C. Mazauric) Agrif implementation
-   !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf
-   !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option
-   !!            3.3  !  2010-06  (C. Ethe, G. Madec) Merge TRA-TRC
-   !!----------------------------------------------------------------------
-#if defined key_top
-   !!----------------------------------------------------------------------
-   !!   'key_top'                                                TOP models
-   !!----------------------------------------------------------------------
-   !!   trc_nxt     : time stepping on passive tracers
-   !!----------------------------------------------------------------------
-   USE oce_trc         ! ocean dynamics and tracers variables
-   USE trc             ! ocean passive tracers variables
-   USE trd_oce
-   USE trdtra
-   USE tranxt
-   USE bdy_oce   , ONLY: ln_bdy
-   USE trcbdy          ! BDY open boundaries
-# if defined key_agrif
-   USE agrif_top_interp
-# endif
-   !
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE prtctl_trc      ! Print control for debbuging
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   trc_nxt   ! routine called by step.F90
-
-   REAL(wp) ::   rfact1, rfact2
-
-   !!----------------------------------------------------------------------
-   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
-   !! $Id$ 
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE trc_nxt( kt, Kbb, Kmm, Krhs )
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE trcnxt  ***
-      !!
-      !! ** Purpose :   Compute the passive tracers fields at the 
-      !!      next time-step from their temporal trends and swap the fields.
-      !! 
-      !! ** Method  :   Apply lateral boundary conditions on (uu(Krhs),vv(Krhs)) through 
-      !!      call to lbc_lnk routine
-      !!   default:
-      !!      arrays swap
-      !!         (tr(Kmm)) = (tr(Krhs)) ; (tr(Krhs)) = (0,0)
-      !!         (tr(Kbb)) = (tr(Kmm)) 
-      !!
-      !!   For Arakawa or TVD Scheme : 
-      !!      A Asselin time filter applied on now tracers (tr(:,:,:,:,Kmm)) to avoid
-      !!      the divergence of two consecutive time-steps and tr arrays
-      !!      to prepare the next time_step:
-      !!         (tr(Kbb)) = (tr(Kmm)) + atfp [ (tr(Kbb)) + (tr(Krhs)) - 2 (tr(Kmm)) ]
-      !!         (tr(Kmm)) = (tr(Krhs)) ; (tr(Krhs)) = (0,0)
-      !!
-      !!
-      !! ** Action  : - update tr(Kbb), tr(Kmm)
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
-      INTEGER, INTENT( in ) ::   Kbb, Kmm, Krhs  ! time level indices
-      !
-      INTEGER  ::   jk, jn   ! dummy loop indices
-      REAL(wp) ::   zfact            ! temporary scalar
-      CHARACTER (len=22) :: charout
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrdt    ! 4D workspace
-      !!----------------------------------------------------------------------
-      !
-      IF( ln_timing )   CALL timing_start('trc_nxt')
-      !
-      IF( kt == nittrc000 .AND. lwp ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'trc_nxt : time stepping on passive tracers'
-      ENDIF
-      !
-#if defined key_agrif
-      CALL Agrif_trc                   ! AGRIF zoom boundaries
-#endif
-      ! Update after tracer on domain lateral boundaries
-      CALL lbc_lnk( 'trcnxt', tr(:,:,:,:,Krhs), 'T', 1. )   
-
-      IF( ln_bdy )  CALL trc_bdy( kt, Kbb, Kmm, Krhs )
-
-      IF( l_trdtrc )  THEN             ! trends: store now fields before the Asselin filter application
-         ALLOCATE( ztrdt(jpi,jpj,jpk,jptra) )
-         ztrdt(:,:,:,:)  = 0._wp
-         IF( ln_traldf_iso ) THEN                       ! diagnose the "pure" Kz diffusive trend 
-            DO jn = 1, jptra
-               CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_zdfp, ztrdt(:,:,:,jn) )
-            ENDDO
-         ENDIF
-
-         ! total trend for the non-time-filtered variables. 
-         zfact = 1.0 / rdttrc
-         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from ts(Kmm) terms
-         IF( ln_linssh ) THEN       ! linear sea surface height only
-            DO jn = 1, jptra
-               DO jk = 1, jpkm1
-                  ztrdt(:,:,jk,jn) = ( tr(:,:,jk,jn,Krhs)*e3t(:,:,jk,Krhs) / e3t(:,:,jk,Kmm) - tr(:,:,jk,jn,Kmm)) * zfact
-               END DO
-            END DO
-         ELSE
-            DO jn = 1, jptra
-               DO jk = 1, jpkm1
-                  ztrdt(:,:,jk,jn) = ( tr(:,:,jk,jn,Krhs) - tr(:,:,jk,jn,Kmm) ) * zfact
-               END DO
-            END DO
-         ENDIF
-         !
-         DO jn = 1, jptra
-            CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_tot, ztrdt(:,:,:,jn) )
-         ENDDO
-         !
-         IF( ln_linssh ) THEN       ! linear sea surface height only
-            ! Store now fields before applying the Asselin filter 
-            ! in order to calculate Asselin filter trend later.
-            ztrdt(:,:,:,:) = tr(:,:,:,:,Kmm) 
-         ENDIF
-
-      ENDIF
-      !                                ! Leap-Frog + Asselin filter time stepping
-      IF( (neuler == 0 .AND. kt == nittrc000) .OR. ln_top_euler ) THEN    ! Euler time-stepping (only swap)
-         DO jn = 1, jptra
-            DO jk = 1, jpkm1
-               tr(:,:,jk,jn,Kmm) = tr(:,:,jk,jn,Krhs)
-               tr(:,:,jk,jn,Kbb) = tr(:,:,jk,jn,Kmm)  
-            END DO
-         END DO
-         IF (l_trdtrc .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl
-            !                                        ! Asselin filter is output by tra_nxt_vvl that is not called on this time step
-            ztrdt(:,:,:,:) = 0._wp            
-            DO jn = 1, jptra
-               CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
-            ENDDO
-         END IF
-         !
-      ELSE     
-         IF( .NOT. l_offline ) THEN ! Leap-Frog + Asselin filter time stepping
-            IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt,      Kmm,       nittrc000,         'TRC',                &
-              &                                              tr(:,:,:,:,Kbb), tr(:,:,:,:,Kmm), tr(:,:,:,:,Krhs), jptra )  !     linear ssh
-            ELSE                   ;   CALL tra_nxt_vvl( kt, Kbb, Kmm, Krhs, nittrc000, rdttrc, 'TRC',                &
-              &                                              tr(:,:,:,:,Kbb), tr(:,:,:,:,Kmm), tr(:,:,:,:,Krhs),      &
-              &                                              sbc_trc, sbc_trc_b, jptra )                                  ! non-linear ssh
-            ENDIF
-         ELSE
-                                       CALL trc_nxt_off( kt, Kbb, Kmm, Krhs )       ! offline 
-         ENDIF
-         !
-         CALL lbc_lnk_multi( 'trcnxt', tr(:,:,:,:,Kbb), 'T', 1._wp, tr(:,:,:,:,Kmm), 'T', 1._wp, tr(:,:,:,:,Krhs), 'T', 1._wp )
-      ENDIF
-      !
-      IF( l_trdtrc .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt )
-         DO jn = 1, jptra
-            DO jk = 1, jpkm1
-               zfact = 1._wp / r2dttrc  
-               ztrdt(:,:,jk,jn) = ( tr(:,:,jk,jn,Kbb) - ztrdt(:,:,jk,jn) ) * zfact 
-            END DO
-            CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
-         END DO
-      END IF
-      IF( l_trdtrc ) DEALLOCATE( ztrdt ) 
-      !
-      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
-         WRITE(charout, FMT="('nxt')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Kmm), mask=tmask, clinfo=ctrcnm)
-      ENDIF
-      !
-      IF( ln_timing )   CALL timing_stop('trc_nxt')
-      !
-   END SUBROUTINE trc_nxt
-
-
-   SUBROUTINE trc_nxt_off( kt, Kbb, Kmm, Krhs )
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE tra_nxt_vvl  ***
-      !!
-      !! ** Purpose :   Time varying volume: apply the Asselin time filter  
-      !!                and swap the tracer fields.
-      !! 
-      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
-      !!              - save in (ta,sa) a thickness weighted average over the three 
-      !!             time levels which will be used to compute rdn and thus the semi-
-      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
-      !!              - swap tracer fields to prepare the next time_step.
-      !!                This can be summurized for tempearture as:
-      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
-      !!                  /( e3t(:,:,:,Kmm)    + rbcp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Krhs)    ] )   
-      !!             ztm = 0                                                       otherwise
-      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
-      !!                  /( e3t(:,:,:,Kmm)    + atfp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Krhs)    ] )
-      !!             tn  = ta 
-      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
-      !!
-      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
-      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in   )   ::  kt              ! ocean time-step index
-      INTEGER, INTENT(in   )   ::  Kbb, Kmm, Krhs  ! time level indices
-      !!     
-      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
-      REAL(wp) ::   ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
-      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
-      !!----------------------------------------------------------------------
-      !
-      IF( kt == nittrc000 )  THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'trc_nxt_off : time stepping'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-         IF( .NOT. ln_linssh ) THEN
-            rfact1 = atfp * rdttrc
-            rfact2 = rfact1 / rau0
-         ENDIF
-        !  
-      ENDIF
-      !
-      DO jn = 1, jptra      
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ze3t_b = e3t(ji,jj,jk,Kbb)
-                  ze3t_n = e3t(ji,jj,jk,Kmm)
-                  ze3t_a = e3t(ji,jj,jk,Krhs)
-                  !                                         ! tracer content at Before, now and after
-                  ztc_b  = tr(ji,jj,jk,jn,Kbb)  * ze3t_b
-                  ztc_n  = tr(ji,jj,jk,jn,Kmm)  * ze3t_n
-                  ztc_a  = tr(ji,jj,jk,jn,Krhs) * ze3t_a
-                  !
-                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
-                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
-                  !
-                  ze3t_f = ze3t_n + atfp * ze3t_d
-                  ztc_f  = ztc_n  + atfp * ztc_d
-                  !
-                  IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN           ! first level 
-                     ze3t_f = ze3t_f - rfact2 * ( emp_b(ji,jj)      - emp(ji,jj)   ) 
-                     ztc_f  = ztc_f  - rfact1 * ( sbc_trc(ji,jj,jn) - sbc_trc_b(ji,jj,jn) )
-                  ENDIF
-
-                  ze3t_f = 1.e0 / ze3t_f
-                  tr(ji,jj,jk,jn,Kbb) = ztc_f * ze3t_f           ! pt(:,:,:,:,Kbb) <-- pt(:,:,:,:,Kmm) filtered
-                  tr(ji,jj,jk,jn,Kmm) = tr(ji,jj,jk,jn,Krhs)     ! pt(:,:,:,:,Kmm) <-- pt(:,:,:,:,Krhs)
-                  !
-               END DO
-            END DO
-         END DO
-         ! 
-      END DO
-      !
-   END SUBROUTINE trc_nxt_off
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option                                         Empty module
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE trc_nxt( kt, Kbb, Kmm, Krhs )  
-      INTEGER, INTENT(in) :: kt
-      INTEGER, INTENT(in) :: Kbb, Kmm, Krhs  ! time level indices
-      WRITE(*,*) 'trc_nxt: You should not have seen this print! error?', kt
-   END SUBROUTINE trc_nxt
-#endif
-   !!======================================================================
-END MODULE trcnxt
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trctrp.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trctrp.F90	(revision 11056)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trctrp.F90	(revision 11057)
@@ -20,5 +20,5 @@
    USE trcadv          ! advection                           (trc_adv routine)
    USE trczdf          ! vertical diffusion                  (trc_zdf routine)
-   USE trcnxt          ! time-stepping                       (trc_nxt routine)
+   USE trcatf          ! time filtering                      (trc_atf routine)
    USE trcrad          ! positivity                          (trc_rad routine)
    USE trcsbc          ! surface boundary condition          (trc_sbc routine)
@@ -78,5 +78,5 @@
 #endif
                                 CALL trc_zdf    ( kt, Kbb, Kmm, Krhs, tr, Kaa  )  ! vert. mixing & after tracer	==> after
-                                CALL trc_nxt    ( kt, Kbb, Kmm, Krhs )            ! tracer fields at next time step     
+                                CALL trc_atf    ( kt, Kbb, Kmm, Kaa , tr )        ! time filtering of "now" tracer fields    
          IF( ln_trcrad )        CALL trc_rad    ( kt, Kbb, Kmm, Krhs, tr       )  ! Correct artificial negative concentrations
          IF( ln_trcdmp_clo )    CALL trc_dmp_clo( kt, Kbb, Kmm )                  ! internal damping trends on closed seas only
@@ -87,5 +87,5 @@
          IF( ln_trcdmp )        CALL trc_dmp( kt, Kbb, Kmm,       tr, Krhs )  ! internal damping trends
                                 CALL trc_zdf( kt, Kbb, Kmm, Krhs, tr, Kaa  )  ! vert. mixing & after tracer	==> after
-                                CALL trc_nxt( kt, Kbb, Kmm, Krhs )            ! tracer fields at next time step     
+                                CALL trc_atf( kt, Kbb, Kmm, Kaa , tr )        ! time filtering of "now" tracer fields
           IF( ln_trcrad )       CALL trc_rad( kt, Kbb, Kmm, Krhs, tr       )  ! Correct artificial negative concentrations
          !
