Index: NEMO/trunk/src/OCE/C1D/step_c1d.F90
===================================================================
--- NEMO/trunk/src/OCE/C1D/step_c1d.F90	(revision 10068)
+++ NEMO/trunk/src/OCE/C1D/step_c1d.F90	(revision 12377)
@@ -6,4 +6,5 @@
    !! History :   2.0  !  2004-04  (C. Ethe)  adapted from step.F90 for C1D
    !!             3.0  !  2008-04  (G. Madec)  redo the adaptation to include SBC
+   !!             4.1  !  2019-08  (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme
    !!----------------------------------------------------------------------
 #if defined key_c1d
@@ -14,9 +15,10 @@
    !!----------------------------------------------------------------------
    USE step_oce        ! time stepping definition modules 
+   USE step, ONLY : Nbb, Nnn, Naa, Nrhs ! time level indices
 #if defined key_top
    USE trcstp          ! passive tracer time-stepping      (trc_stp routine)
 #endif
    USE dyncor_c1d      ! Coriolis term (c1d case)         (dyn_cor_1d     )
-   USE dynnxt          ! time-stepping                    (dyn_nxt routine)
+   USE dynatf          ! time filtering                   (dyn_atf routine)
    USE dyndmp          ! U & V momentum damping           (dyn_dmp routine)
    USE restart         ! restart 
@@ -65,26 +67,26 @@
       ! Update data, open boundaries, surface boundary condition (including sea-ice)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
+                         CALL sbc    ( kstp, Nbb, Nnn )  ! Sea Boundary Condition (including sea-ice)
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Ocean physics update                (ua, va, ta, sa used as workspace)
+      ! Ocean physics update        
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL eos_rab( tsb, rab_b )   ! before local thermal/haline expension ratio at T-points
-                         CALL eos_rab( tsn, rab_n )   ! now    local thermal/haline expension ratio at T-points
-                         CALL bn2( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency
-                         CALL bn2( tsn, rab_n, rn2  ) ! now    Brunt-Vaisala frequency
+                         CALL eos_rab( ts(:,:,:,:,Nbb), rab_b, Nnn )  ! before local thermal/haline expension ratio at T-points
+                         CALL eos_rab( ts(:,:,:,:,Nnn), rab_n, Nnn )  ! now    local thermal/haline expension ratio at T-points
+                         CALL bn2( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency
+                         CALL bn2( ts(:,:,:,:,Nnn), rab_n, rn2 , Nnn ) ! now    Brunt-Vaisala frequency
       
       !  VERTICAL PHYSICS
-                         CALL zdf_phy( kstp )         ! vertical physics update (bfr, avt, avs, avm + MLD)
+                         CALL zdf_phy( kstp, Nbb, Nnn, Nrhs  )    ! vertical physics update (bfr, avt, avs, avm + MLD)
 
-      IF(.NOT.ln_linssh )   CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_hor)
-      IF(.NOT.ln_linssh )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors 
+      IF(.NOT.ln_linssh )   CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )  ! after ssh (includes call to div_hor)
+      IF(.NOT.ln_linssh )   CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn,      Naa )  ! after vertical scale factors 
 
-      IF(.NOT.ln_linssh )   CALL wzv           ( kstp )  ! now cross-level velocity 
+      IF(.NOT.ln_linssh )   CALL wzv           ( kstp, Nbb, Nnn, ww,  Naa )  ! now cross-level velocity 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! diagnostics and outputs             (ua, va, ta, sa used as workspace)
+      ! diagnostics and outputs       
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL dia_wri( kstp )       ! ocean model: outputs
-      IF( lk_diahth  )   CALL dia_hth( kstp )       ! Thermocline depth (20°C)
+                         CALL dia_wri( kstp, Nnn )  ! ocean model: outputs
+      IF( lk_diahth  )   CALL dia_hth( kstp, Nnn )  ! Thermocline depth (20°C)
 
 
@@ -93,21 +95,21 @@
       ! Passive Tracer Model
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                        CALL trc_stp( kstp )       ! time-stepping
+                        CALL trc_stp( kstp, Nbb, Nnn, Nrhs, Naa  )   ! time-stepping
 #endif
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Active tracers                              (ua, va used as workspace)
+      ! Active tracers                              (uu(:,:,:,Nrhs), vv(:,:,:,Nrhs) used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                        tsa(:,:,:,:) = 0._wp       ! set tracer trends to zero
+                        ts(:,:,:,:,Nrhs) = 0._wp       ! set tracer trends to zero
 
-                        CALL tra_sbc( kstp )       ! surface boundary condition
-      IF( ln_traqsr )   CALL tra_qsr( kstp )       ! penetrative solar radiation qsr
-      IF( ln_tradmp )   CALL tra_dmp( kstp )       ! internal damping trends- tracers
-      IF(.NOT.ln_linssh)CALL tra_adv( kstp )       ! horizontal & vertical advection
-      IF( ln_zdfosm  )  CALL tra_osm( kstp )       ! OSMOSIS non-local tracer fluxes
-                        CALL tra_zdf( kstp )       ! vertical mixing
-                        CALL eos( tsn, rhd, rhop, gdept_0(:,:,:) )   ! now potential density for zdfmxl
-      IF( ln_zdfnpc )   CALL tra_npc( kstp )       ! applied non penetrative convective adjustment on (t,s)
-                        CALL tra_nxt( kstp )       ! tracer fields at next time step
+                        CALL tra_sbc( kstp,      Nnn, ts, Nrhs  )  ! surface boundary condition
+      IF( ln_traqsr )   CALL tra_qsr( kstp,      Nnn, ts, Nrhs  )  ! penetrative solar radiation qsr
+      IF( ln_tradmp )   CALL tra_dmp( kstp, Nbb, Nnn, ts, Nrhs  )  ! internal damping trends- tracers
+      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 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_atf( kstp, Nbb, Nnn, Nrhs,     Naa, ts   )     ! time filtering of "now" tracer fields
 
 
@@ -116,22 +118,28 @@
       ! Dynamics                                    (ta, sa used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                        ua(:,:,:) = 0._wp          ! set dynamics trends to zero
-                        va(:,:,:) = 0._wp
+                        uu(:,:,:,Nrhs) = 0._wp          ! set dynamics trends to zero
+                        vv(:,:,:,Nrhs) = 0._wp
 
-      IF( ln_dyndmp )   CALL dyn_dmp    ( kstp )   ! internal damping trends- momentum
-                        CALL dyn_cor_c1d( kstp )   ! vorticity term including Coriolis
-      IF( ln_zdfosm  )  CALL dyn_osm    ( kstp )   ! OSMOSIS non-local velocity fluxes
-                        CALL dyn_zdf    ( kstp )   ! vertical diffusion
-                        CALL dyn_nxt    ( kstp )   ! lateral velocity at next time step
-      IF(.NOT.ln_linssh)CALL ssh_swp    ( kstp )   ! swap of sea surface height
-
-      IF(.NOT.ln_linssh)CALL dom_vvl_sf_swp( kstp )! swap of vertical scale factors
+      IF( ln_dyndmp )   CALL dyn_dmp    ( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! internal damping trends- momentum
+                        CALL dyn_cor_c1d( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity term including Coriolis
+      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_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_update( kstp, Nbb, Nnn, Naa )                    ! update of vertical scale factors
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Control and restarts
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                             CALL stp_ctl( kstp, indic )
-      IF( kstp == nit000 )   CALL iom_close( numror )      ! close input  ocean restart file
-      IF( lrst_oce       )   CALL rst_write( kstp )        ! write output ocean restart file
+                             CALL stp_ctl( kstp, Nnn, indic )
+      IF( kstp == nit000 )   CALL iom_close( numror )          ! close input  ocean restart file
+      IF( lrst_oce       )   CALL rst_write( kstp, Nbb, Nnn )  ! write output ocean restart file
       !
 #if defined key_iomput
