Index: CONFIG/UNIFORM/v6/NEMO_v6.5/SOURCES/p4zbc.F90
===================================================================
--- CONFIG/UNIFORM/v6/NEMO_v6.5/SOURCES/p4zbc.F90	(revision 6877)
+++ CONFIG/UNIFORM/v6/NEMO_v6.5/SOURCES/p4zbc.F90	(revision 6877)
@@ -0,0 +1,388 @@
+MODULE p4zbc
+   !!======================================================================
+   !!                         ***  MODULE p4sbc  ***
+   !! TOP :   PISCES surface boundary conditions of external inputs of nutrients
+   !!======================================================================
+   !! History :   3.5  !  2012-07 (O. Aumont, C. Ethe) Original code
+   !!----------------------------------------------------------------------
+   !!   p4z_bc        :  Read and interpolate time-varying nutrients fluxes
+   !!   p4z_bc_init   :  Initialization of p4z_bc
+   !!----------------------------------------------------------------------
+   USE oce_trc         !  shared variables between ocean and passive tracers
+   USE trc             !  passive tracers common variables 
+   USE sms_pisces      !  PISCES Source Minus Sink variables
+   USE iom             !  I/O manager
+   USE fldread         !  time interpolation
+   USE trcbc
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   p4z_bc
+   PUBLIC   p4z_bc_init   
+
+   LOGICAL , PUBLIC ::   ln_ironsed   !: boolean for Fe input from sediments
+   LOGICAL , PUBLIC ::   ln_hydrofe   !: boolean for Fe input from hydrothermal vents
+   REAL(wp), PUBLIC ::   sedfeinput   !: Coastal release of Iron
+   REAL(wp), PUBLIC ::   icefeinput   !: Iron concentration in sea ice
+   REAL(wp), PUBLIC ::   wdust        !: Sinking speed of the dust 
+   REAL(wp), PUBLIC ::   mfrac        !: Mineral Content of the dust
+   REAL(wp)         ::   hratio       !: Fe:3He ratio assumed for vent iron supply
+   REAL(wp)         ::   distcoast    !: Distance off the coast for Iron from sediments
+   REAL(wp), PUBLIC ::   lgw_rath     !: Weak ligand ratio from hydro sources
+
+   LOGICAL , PUBLIC ::   ll_bc
+   LOGICAL , PUBLIC ::   ll_dust      !: boolean for dust input from the atmosphere
+   LOGICAL , PUBLIC ::   ll_river     !: boolean for river input of nutrients
+   LOGICAL , PUBLIC ::   ll_ndepo     !: boolean for atmospheric deposition of N
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_dust      ! structure of input dust
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_ironsed   ! structure of input iron from sediment
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_hydrofe   ! structure of input iron from sediment
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dust    !: dust fields
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ironsed          !: Coastal supply of iron
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hydrofe          !: Hydrothermal vent supply of iron
+
+   REAL(wp), PUBLIC :: sedsilfrac, sedcalfrac
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
+   !! $Id: p4zbc.F90 10869 2019-04-15 10:34:03Z cetlod $ 
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE p4z_bc( kt, Kbb, Kmm, Krhs )
+      !!----------------------------------------------------------------------
+      !!                  ***  routine p4z_bc  ***
+      !!
+      !! ** purpose :   read and interpolate the external sources of nutrients
+      !!
+      !! ** method  :   read the files and interpolate the appropriate variables
+      !!
+      !! ** input   :   external netcdf files
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt              ! ocean time step
+      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs  ! time level index
+      !
+      INTEGER  ::  ji, jj, jk, jl 
+      REAL(wp) ::  zdep, zwflux, zironice
+      REAL(wp) ::  zcoef, zwdust, zrivdin, zdustdep, zndep
+      !
+      CHARACTER (len=25) :: charout
+      !!---------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('p4z_bc')
+      
+      ! Add the external input of nutrients from dust deposition in the water column
+      ! The inputs at surface have already been added
+      ! ----------------------------------------------------------
+      IF( ll_dust )  THEN
+         !
+         CALL fld_read( kt, 1, sf_dust )
+         dust(:,:) = MAX( rtrn, sf_dust(1)%fnow(:,:,1) )
+         !
+         ! Iron solubilization of particles in the water column
+         ! dust in kg/m2/s ---> 1/55.85 to put in mol/Fe ;  wdust in m/d
+         ! Dust are supposed to sink at wdust sinking speed. 3% of the iron 
+         ! in dust is hypothesized to be soluble at a dissolution rate set to 
+         ! 1/(250 days). The vertical distribution of iron in dust is computed 
+         ! from a steady state assumption. Parameters are very uncertain and 
+         ! are estimated from the literature quoted in Raiswell et al. (2011) 
+         ! ------------------------------------------------------------------- 
+
+         zwdust = 0.03 / ( wdust / rday ) / ( 250. * rday )
+
+         ! Atmospheric input of Iron dissolves in the water column
+         IF ( ln_trc_sbc(jpfer) ) THEN
+            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1 )
+               zdustdep = dust(ji,jj) * zwdust * rfact * EXP( -gdept(ji,jj,jk,Kmm) /( 250. * wdust ) )
+               !
+               tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zdustdep * mfrac / mMass_Fe 
+            END_3D
+
+            IF( lk_iomput ) THEN
+                ! surface downward dust depo of iron
+                jl = n_trc_indsbc(jpfer)
+                CALL iom_put( "Irondep", ( rf_trsfac(jl) * sf_trcsbc(jl)%fnow(:,:,1) / rn_sbc_time ) * 1.e+3 * tmask(:,:,1) )
+
+            ENDIF
+
+         ENDIF
+
+         ! Atmospheric input of PO4 dissolves in the water column
+         IF ( ln_trc_sbc(jppo4) ) THEN
+            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1 )
+               zdustdep = dust(ji,jj) * zwdust * rfact * EXP( -gdept(ji,jj,jk,Kmm) /( 250. * wdust ) )
+               !
+               tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) + zdustdep * 1.e-3 / mMass_P
+            END_3D
+         ENDIF
+
+         ! Atmospheric input of Si dissolves in the water column
+         IF ( ln_trc_sbc(jpsil) ) THEN
+            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1 )
+               zdustdep = dust(ji,jj) * zwdust * rfact * EXP( -gdept(ji,jj,jk,Kmm) /( 250. * wdust ) )
+               !
+               tr(ji,jj,jk,jpsil,Krhs) = tr(ji,jj,jk,jpsil,Krhs) + zdustdep * 0.269 / mMass_Si
+            END_3D
+         ENDIF
+
+         !
+         IF( lk_iomput ) THEN
+             ! dust concentration at surface
+             CALL iom_put( "pdust"  , dust(:,:) / ( wdust / rday ) * tmask(:,:,1) ) ! dust concentration at surface
+         ENDIF
+      ENDIF
+
+      ! -----------------------------------------
+      ! Add the external input of nutrients from river
+      ! ----------------------------------------------------------
+      IF( ll_river ) THEN
+          jl = n_trc_indcbc(jpno3)
+          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+             DO jk = 1, nk_rnf(ji,jj)
+                zcoef = rn_rfact / ( e1e2t(ji,jj) * h_rnf(ji,jj) * rn_cbc_time ) * tmask(ji,jj,1)
+                zrivdin = rf_trcfac(jl) * sf_trccbc(jl)%fnow(ji,jj,1) * zcoef
+                tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) - rno3 * zrivdin * rfact
+            ENDDO
+          END_2D
+      ENDIF
+      
+      ! Add the external input of nutrients from nitrogen deposition
+      ! ----------------------------------------------------------
+      IF( ll_ndepo ) THEN
+         IF( ln_trc_sbc(jpno3) ) THEN
+            jl = n_trc_indsbc(jpno3)
+            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+               zndep = rf_trsfac(jl) * sf_trcsbc(jl)%fnow(ji,jj,1) / e3t(ji,jj,1,Kmm) / rn_sbc_time
+               tr(ji,jj,1,jptal,Krhs) = tr(ji,jj,1,jptal,Krhs) - rno3 * zndep * rfact
+            END_2D
+         ENDIF
+         IF( ln_trc_sbc(jpnh4) ) THEN
+            jl = n_trc_indsbc(jpnh4)
+            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+               zndep = rf_trsfac(jl) * sf_trcsbc(jl)%fnow(ji,jj,1) / e3t(ji,jj,1,Kmm) / rn_sbc_time
+               tr(ji,jj,1,jptal,Krhs) = tr(ji,jj,1,jptal,Krhs) + rno3 * zndep * rfact
+            END_2D
+         ENDIF
+      ENDIF
+      !
+      ! Iron input/uptake due to sea ice : Crude parameterization based on 
+      ! Lancelot et al. Iron concentration in sea-ice is constant and set 
+      ! in the namelist_pisces (icefeinput). ln_ironice is forced to false
+      ! when nn_ice_tr = 1
+      ! ----------------------------------------------------
+      IF( ln_ironice ) THEN
+         !
+         ! Compute the iron flux between sea ice and sea water
+         ! Simple parameterization assuming a fixed constant concentration in
+         ! sea-ice (icefeinput)
+         ! ------------------------------------------------------------------         
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+            zdep     = rfact / e3t(ji,jj,1,Kmm)
+            zwflux   = fmmflx(ji,jj) / 1000._wp
+            zironice =  MAX( -0.99 * tr(ji,jj,1,jpfer,Kbb), -zwflux * icefeinput * zdep )
+            tr(ji,jj,1,jpfer,Krhs) = tr(ji,jj,1,jpfer,Krhs) + zironice
+         END_2D
+         !
+         ! iron flux from ice
+         IF( lk_iomput ) &
+         & CALL iom_put( "Ironice", MAX( -0.99 * tr(:,:,1,jpfer,Kbb), (-1.*fmmflx(:,:)/1000._wp )*icefeinput*1.e+3*tmask(:,:,1)) )
+         !
+      ENDIF
+
+      ! Add the external input of iron from sediment mobilization
+      ! ------------------------------------------------------
+      IF( ln_ironsed .AND. .NOT.lk_sed ) THEN
+          tr(:,:,:,jpfer,Krhs) = tr(:,:,:,jpfer,Krhs) + ironsed(:,:,:) * rfact
+          !
+          IF( lk_iomput )  CALL iom_put( "Ironsed", ironsed(:,:,:) * 1.e+3 * tmask(:,:,:) ) 
+      ENDIF
+
+      ! Add the external input of iron from hydrothermal vents
+      ! ------------------------------------------------------
+      IF( ln_hydrofe ) THEN
+         CALL fld_read( kt, 1, sf_hydrofe )
+         DO jk = 1, jpk
+            hydrofe(:,:,jk) = ( MAX( rtrn, sf_hydrofe(1)%fnow(:,:,jk) ) * hratio ) &
+              &              / ( e1e2t(:,:) * e3t(:,:,jk,Kmm) * ryyss + rtrn ) / 1000._wp &
+              &              * tmask(:,:,jk)
+         ENDDO
+                         tr(:,:,:,jpfer,Krhs) = tr(:,:,:,jpfer,Krhs) + hydrofe(:,:,:) * rfact
+         IF( ln_ligand ) tr(:,:,:,jplgw,Krhs) = tr(:,:,:,jplgw,Krhs) + ( hydrofe(:,:,:) * lgw_rath ) * rfact
+         !
+         IF( lk_iomput ) CALL iom_put( "HYDR", hydrofe(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! hydrothermal iron input
+      ENDIF
+      IF( ln_timing )  CALL timing_stop('p4z_bc')
+      !
+   END SUBROUTINE p4z_bc
+
+
+   SUBROUTINE p4z_bc_init( Kmm ) 
+      !!----------------------------------------------------------------------
+      !!                  ***  routine p4z_bc_init  ***
+      !!
+      !! ** purpose :   initialization of the external sources of nutrients
+      !!
+      !! ** method  :   read the files and compute the budget
+      !!                called at the first timestep (nittrc000)
+      !!
+      !! ** input   :   external netcdf files
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   Kmm  ! time level index
+      INTEGER  :: ji, jj, jk, jm
+      INTEGER  :: ii0, ii1, ij0, ij1
+      INTEGER  :: numiron
+      INTEGER  :: ierr, ierr1, ierr2, ierr3
+      INTEGER  :: ios                 ! Local integer output status for namelist read
+      INTEGER  :: ik50                !  last level where depth less than 50 m
+      REAL(wp) :: zexpide, zdenitide, zmaskt, zsurfc, zsurfp,ze3t, ze3t2, zcslp
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zriver, zcmask
+      !
+      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
+      TYPE(FLD_N) ::   sn_dust, sn_ironsed, sn_hydrofe   ! informations about the fields to be read
+      !!
+      NAMELIST/nampisbc/cn_dir, sn_dust, sn_ironsed, sn_hydrofe, &
+        &                ln_ironsed, ln_ironice, ln_hydrofe,     &
+        &                sedfeinput, distcoast, icefeinput, wdust, mfrac,   &
+        &                hratio, lgw_rath
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'p4z_bc_init : initialization of the external sources of nutrients '
+         WRITE(numout,*) '~~~~~~~~~~~~ '
+      ENDIF
+      !                            !* set file information
+      READ  ( numnatp_ref, nampisbc, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampisbc in reference namelist' )
+      READ  ( numnatp_cfg, nampisbc, IOSTAT = ios, ERR = 902 )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisbc in configuration namelist' )
+      IF(lwm) WRITE ( numonp, nampisbc )
+
+
+      IF(lwp) THEN
+         WRITE(numout,*) '   Namelist : nampissbc '
+         WRITE(numout,*) '      Fe input from sediments                  ln_ironsed  = ', ln_ironsed
+         WRITE(numout,*) '      Fe input from seaice                     ln_ironice  = ', ln_ironice
+         WRITE(numout,*) '      fe input from hydrothermal vents         ln_hydrofe  = ', ln_hydrofe
+         IF( ln_ironsed ) THEN
+            WRITE(numout,*) '      coastal release of iron                  sedfeinput  = ', sedfeinput
+            WRITE(numout,*) '      distance off the coast                   distcoast   = ', distcoast
+         ENDIF
+         IF( ln_ligand ) THEN
+            WRITE(numout,*) '      Weak ligand ratio from sed hydro sources  lgw_rath   = ', lgw_rath
+         ENDIF
+         IF( ln_ironice ) THEN
+            WRITE(numout,*) '      Iron concentration in sea ice            icefeinput  = ', icefeinput
+         ENDIF
+         IF( ln_trc_sbc(jpfer) ) THEN
+            WRITE(numout,*) '      Mineral Fe content of the dust           mfrac       = ', mfrac
+            WRITE(numout,*) '      sinking speed of the dust                wdust       = ', wdust
+         ENDIF
+         IF( ln_hydrofe ) THEN
+            WRITE(numout,*) '      Fe to 3He ratio assumed for vent iron supply hratio  = ', hratio
+         ENDIF
+      END IF
+
+      ll_bc    = ( ln_trcbc .AND. lltrcbc )  .OR. ln_hydrofe .OR. ln_ironsed .OR. ln_ironice .OR. ln_sediment 
+      ll_dust  =  ln_trc_sbc(jpfer) .OR. ln_trc_sbc(jppo4) .OR. ln_trc_sbc(jpsil) .OR. ln_sediment
+      ll_ndepo =  ln_trc_sbc(jpno3) .OR. ln_trc_sbc(jpnh4)   
+      ll_river =  ln_trc_cbc(jpno3)  
+
+      ! dust input from the atmosphere
+      ! ------------------------------
+      IF( ll_dust ) THEN
+         !
+         IF(lwp) WRITE(numout,*) '    initialize dust input from atmosphere '
+         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
+         !
+         ALLOCATE( dust(jpi,jpj) ) 
+         !
+         ALLOCATE( sf_dust(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
+         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_bc_init: unable to allocate sf_dust structure' )
+         !
+         CALL fld_fill( sf_dust, (/ sn_dust /), cn_dir, 'p4z_bc_init', 'Atmospheric dust deposition', 'nampisbc' )
+                                   ALLOCATE( sf_dust(1)%fnow(jpi,jpj,1)   )
+         IF( sn_dust%ln_tint )     ALLOCATE( sf_dust(1)%fdta(jpi,jpj,1,2) )
+         !
+      END IF
+
+      ! coastal and island masks
+      ! ------------------------
+      IF( ln_ironsed ) THEN     
+         !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   ==>>>   ln_ironsed=T , computation of an island mask to enhance coastal supply of iron'
+         !
+         ALLOCATE( ironsed(jpi,jpj,jpk) )    ! allocation
+         !
+         CALL iom_open ( TRIM( sn_ironsed%clname ), numiron )
+         ALLOCATE( zcmask(jpi,jpj,jpk) )
+         CALL iom_get  ( numiron, jpdom_global, TRIM( sn_ironsed%clvar ), zcmask(:,:,:), 1 )
+         CALL iom_close( numiron )
+         !
+         ik50 = 5        !  last level where depth less than 50 m
+         DO jk = jpkm1, 1, -1
+            IF( gdept_1d(jk) > 50. )   ik50 = jk - 1
+         END DO
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) ' Level corresponding to 50m depth ',  ik50,' ', gdept_1d(ik50+1)
+         DO_3D( 0, 0, 0, 0, 1, ik50 )
+            ze3t   = e3t_0(ji,jj,jk)
+            zsurfc =  e1u(ji,jj) * ( 1. - umask(ji  ,jj  ,jk) )   &
+                    + e1u(ji,jj) * ( 1. - umask(ji-1,jj  ,jk) )   &
+                    + e2v(ji,jj) * ( 1. - vmask(ji  ,jj  ,jk) )   &
+                    + e2v(ji,jj) * ( 1. - vmask(ji  ,jj-1,jk) )
+            zsurfp = zsurfc * ze3t / e1e2t(ji,jj)
+            ! estimation of the coastal slope : 5 km off the coast
+            ze3t2 = ze3t * ze3t
+            zcslp = SQRT( ( distcoast*distcoast + ze3t2 ) / ze3t2 )
+            !
+            zcmask(ji,jj,jk) = zcmask(ji,jj,jk) + zcslp * zsurfp
+         END_3D
+         !
+         CALL lbc_lnk( 'p4zbc', zcmask , 'T', 1.0_wp )      ! lateral boundary conditions on cmask   (sign unchanged)
+         !
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
+            zexpide   = MIN( 8.,( gdept_0(ji,jj,jk) / 500. )**(-1.5) )
+            zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2
+            zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 )
+         END_3D
+         ! Coastal supply of iron
+         ! -------------------------
+         ironsed(:,:,jpk) = 0._wp
+         DO jk = 1, jpkm1
+            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( e3t_0(:,:,jk) * rday )
+         END DO
+         DEALLOCATE( zcmask)
+      ENDIF
+      !
+      ! Iron from Hydrothermal vents
+      ! ------------------------
+      IF( ln_hydrofe ) THEN
+         !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   ==>>>   ln_hydrofe=T , Input of iron from hydrothermal vents'
+         !
+         ALLOCATE( hydrofe(jpi,jpj,jpk) )    ! allocation
+         !
+         ALLOCATE( sf_hydrofe(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
+         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_bc_init: unable to allocate sf_hydro structure' )
+         !
+         CALL fld_fill( sf_hydrofe, (/ sn_hydrofe /), cn_dir, 'p4z_bc_init', 'Input of iron from hydrothermal vents', 'nampisbc' )
+                                   ALLOCATE( sf_hydrofe(1)%fnow(jpi,jpj,jpk)   )
+         IF( sn_hydrofe%ln_tint )    ALLOCATE( sf_hydrofe(1)%fdta(jpi,jpj,jpk,2) )
+         !
+      ENDIF
+      !
+   END SUBROUTINE p4z_bc_init
+
+   !!======================================================================
+END MODULE p4zbc
Index: CONFIG/UNIFORM/v6/NEMO_v6.5/SOURCES/stpctl.F90
===================================================================
--- CONFIG/UNIFORM/v6/NEMO_v6.5/SOURCES/stpctl.F90	(revision 6877)
+++ CONFIG/UNIFORM/v6/NEMO_v6.5/SOURCES/stpctl.F90	(revision 6877)
@@ -0,0 +1,311 @@
+MODULE stpctl
+   !!======================================================================
+   !!                       ***  MODULE  stpctl  ***
+   !! Ocean run control :  gross check of the ocean time stepping
+   !!======================================================================
+   !! History :  OPA  ! 1991-03  (G. Madec) Original code
+   !!            6.0  ! 1992-06  (M. Imbard)
+   !!            8.0  ! 1997-06  (A.M. Treguier)
+   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
+   !!            2.0  ! 2009-07  (G. Madec)  Add statistic for time-spliting
+   !!            3.7  ! 2016-09  (G. Madec)  Remove solver
+   !!            4.0  ! 2017-04  (G. Madec)  regroup global communications
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   stp_ctl      : Control the run
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain variables 
+   USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables
+   USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy
+   !  
+   USE diawri          ! Standard run outputs       (dia_wri_state routine)
+   USE in_out_manager  ! I/O manager
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp         ! distributed memory computing
+   USE eosbn2, ONLY: ln_SEOS, rn_b0
+   !
+   USE netcdf          ! NetCDF library
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC stp_ctl           ! routine called by step.F90
+
+   INTEGER, PARAMETER         ::   jpvar = 8
+   INTEGER                    ::   nrunid   ! netcdf file id
+   INTEGER, DIMENSION(jpvar)  ::   nvarid   ! netcdf variable id
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: stpctl.F90 15023 2021-06-18 14:35:25Z gsamson $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+#  include "single_precision_substitute.h90"
+#  include "single_precision_substitute.h90"
+CONTAINS
+
+   SUBROUTINE stp_ctl( kt, Kmm )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE stp_ctl  ***
+      !!
+      !! ** Purpose :   Control the run
+      !!
+      !! ** Method  : - Save the time step in numstp
+      !!              - Stop the run IF problem encountered by setting nstop > 0
+      !!                Problems checked: |ssh| maximum larger than 10 m
+      !!                                  |U|   maximum larger than 10 m/s 
+      !!                                  negative sea surface salinity
+      !!
+      !! ** Actions :   "time.step" file = last ocean time-step
+      !!                "run.stat"  file = run statistics
+      !!                 nstop indicator sheared among all local domain
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
+      INTEGER, INTENT(in   ) ::   Kmm      ! ocean time level index
+      !!
+      INTEGER, PARAMETER              ::   jptst = 4
+      INTEGER                         ::   ji                                    ! dummy loop indices
+      INTEGER                         ::   idtime, istatus
+      INTEGER , DIMENSION(jptst)      ::   iareasum, iareamin, iareamax
+      INTEGER , DIMENSION(3,jptst)    ::   iloc                                  ! min/max loc indices
+      REAL(dp)                        ::   zzz, zminsal, zmaxsal                 ! local real
+      REAL(dp), DIMENSION(jpvar+1)    ::   zmax
+      REAL(dp), DIMENSION(jptst)      ::   zmaxlocal
+      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns, ll_0oce
+      LOGICAL, DIMENSION(jpi,jpj,jpk) ::   llmsk
+      CHARACTER(len=20)               ::   clname
+      !!----------------------------------------------------------------------
+      IF( nstop > 0 .AND. ngrdstop > -1 )   RETURN   !   stpctl was already called by a child grid
+      !
+      ll_wrtstp  = ( MOD( kt-nit000, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend )
+      ll_colruns = sn_cfctl%l_runstat .AND. ll_wrtstp .AND. jpnij > 1
+      ll_wrtruns = sn_cfctl%l_runstat .AND. ll_wrtstp .AND. lwm
+      !
+      IF( kt == nit000 ) THEN
+         !
+         IF( lwp ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'stp_ctl : time-stepping control'
+            WRITE(numout,*) '~~~~~~~'
+         ENDIF
+         !                                ! open time.step    ascii file, done only by 1st subdomain
+         IF( lwm )   CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+         !
+         IF( ll_wrtruns ) THEN
+            !                             ! open run.stat     ascii file, done only by 1st subdomain
+            CALL ctl_opn( numrun, 'run.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+            !                             ! open run.stat.nc netcdf file, done only by 1st subdomain
+         ENDIF
+         !
+      ENDIF
+      !
+      !                                   !==              write current time step              ==!
+      !                                   !==  done only by 1st subdomain at writting timestep  ==!
+      IF( lwm .AND. ll_wrtstp ) THEN
+         WRITE ( numstp, '(1x, i8)' )   kt
+         REWIND( numstp )
+      ENDIF
+      !                                   !==            test of local extrema           ==!
+      !                                   !==  done by all processes at every time step  ==!
+      !
+      llmsk(     1:nn_hls,:,:) = .FALSE.                                          ! exclude halos from the checked region
+      llmsk(Nie0+1:   jpi,:,:) = .FALSE.
+      llmsk(:,     1:nn_hls,:) = .FALSE.
+      llmsk(:,Nje0+1:   jpj,:) = .FALSE.
+      !
+      llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0) == 1._wp         ! define only the inner domain
+      !
+      ll_0oce = .NOT. ANY( llmsk(:,:,1) )                                         ! no ocean point in the inner domain?
+      !
+      IF( ll_wd ) THEN
+         zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) )   ! ssh max
+      ELSE
+         zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm)           ), mask = llmsk(:,:,1) )   ! ssh max
+      ENDIF
+      llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
+      zmax(2) = MAXVAL(  ABS( uu(:,:,:,Kmm) ), mask = llmsk )                     ! velocity max (zonal only)
+      llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
+      zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     ! minus salinity max
+      zmax(4) = MAXVAL(  ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     !       salinity max
+      IF( ll_colruns .OR. jpnij == 1 ) THEN     ! following variables are used only in the netcdf file
+         zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  ! minus temperature max
+         zmax(6) = MAXVAL(  ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  !       temperature max
+         IF( ln_zad_Aimp ) THEN
+            zmax(7) = MAXVAL(   Cu_adv(:,:,:)   , mask = llmsk )                  ! partitioning coeff. max
+            llmsk(:,:,:) = wmask(:,:,:) == 1._wp
+            zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = llmsk )                  ! implicit vertical vel. max
+         ELSE
+            zmax(7:8) = 0._wp
+         ENDIF
+      ELSE
+         zmax(5:8) = 0._wp
+      ENDIF
+      zmax(jpvar+1) = REAL( nstop, wp )                                           ! stop indicator
+      !
+      !                                   !==               get global extrema             ==!
+      !                                   !==  done by all processes if writting run.stat  ==!
+      IF( ll_colruns ) THEN
+         zmaxlocal(:) = zmax(1:jptst)
+         CALL mpp_max( "stpctl", zmax )          ! max over the global domain: ok even of ll_0oce = .true. 
+         nstop = NINT( zmax(jpvar+1) )           ! update nstop indicator (now sheared among all local domains)
+      ELSE
+         ! if no ocean point: MAXVAL returns -HUGE => we must overwrite this value to avoid error handling bellow.
+         IF( ll_0oce )   zmax(1:jptst) = (/ 0._wp, 0._wp, -1._wp, 1._wp /)   ! default "valid" values...
+      ENDIF
+      !
+      zmax(3) = -zmax(3)                              ! move back from max(-zz) to min(zz) : easier to manage! 
+      zmax(5) = -zmax(5)                              ! move back from max(-zz) to min(zz) : easier to manage!
+      IF( ll_colruns ) zmaxlocal(3) = -zmaxlocal(3)   ! move back from max(-zz) to min(zz) : easier to manage!
+      !
+      !                                   !==              write "run.stat" files              ==!
+      !                                   !==  done only by 1st subdomain at writting timestep  ==!
+      IF( ll_wrtruns ) THEN
+         WRITE(numrun,9500) kt, zmax(1:jptst)
+         IF( jpnij == 1 ) CALL FLUSH(numrun)
+      END IF
+      !                                   !==               error handling               ==!
+      !                                   !==  done by all processes at every time step  ==!
+      !
+      IF ( ln_SEOS.AND.(rn_b0==0._wp) ) THEN             ! Discard checks on salinity
+         zmaxsal =  HUGE(1._wp)                               ! if not used in eos
+         zminsal = -HUGE(1._wp)
+      ELSE
+         zmaxsal = 100._wp
+         zminsal =   0._wp
+      ENDIF 
+      ! 
+      IF(  zmax(1) >   20._wp .OR.   &                        ! too large sea surface height ( > 20 m )
+         & zmax(2) >   10._wp .OR.   &                        ! too large velocity ( > 10 m/s)
+         & zmax(3) <= zminsal .OR.   &                        ! negative or zero sea surface salinity
+         & zmax(4) >= zmaxsal .OR.   &                        ! too large sea surface salinity ( > 100 )
+         & zmax(4) <  zminsal .OR.   &                        ! too large sea surface salinity (keep this line for sea-ice)
+         & ISNAN( SUM(zmax(1:jptst)) ) .OR.   &               ! NaN encounter in the tests
+         & ABS(   SUM(zmax(1:jptst)) ) > HUGE(1._wp) ) THEN   ! Infinity encounter in the tests
+         !
+         iloc(:,:) = 0
+         IF( ll_colruns ) THEN   ! zmax is global, so it is the same on all subdomains -> no dead lock with mpp_maxloc
+            ! first: close the netcdf file, so we can read it
+            ! get global loc on the min/max
+            llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0 ) == 1._wp         ! define only the inner domain
+            CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,         Kmm)), llmsk(:,:,1), zzz, iloc(1:2,1) )   ! mpp_maxloc ok if mask = F 
+            llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
+ CALL mpp_maxloc( 'stpctl', CASTDP(ABS( uu(:,:,:, Kmm))), llmsk(:,:,:), zzz, iloc(1:3,2) )
+            llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
+            CALL mpp_minloc( 'stpctl',      ts(:,:,:,jp_sal,Kmm) , llmsk(:,:,:), zzz, iloc(1:3,3) )
+            CALL mpp_maxloc( 'stpctl',      CASTDP(ts(:,:,:,jp_sal,Kmm)) , llmsk(:,:,:), zzz, iloc(1:3,4) )
+            ! find which subdomain has the max.
+            iareamin(:) = jpnij+1   ;   iareamax(:) = 0   ;   iareasum(:) = 0
+            DO ji = 1, jptst
+               IF( zmaxlocal(ji) == zmax(ji) ) THEN
+                  iareamin(ji) = narea   ;   iareamax(ji) = narea   ;   iareasum(ji) = 1
+               ENDIF
+            END DO
+            CALL mpp_min( "stpctl", iareamin )         ! min over the global domain
+            CALL mpp_max( "stpctl", iareamax )         ! max over the global domain
+            CALL mpp_sum( "stpctl", iareasum )         ! sum over the global domain
+         ELSE                    ! find local min and max locations:
+            ! if we are here, this means that the subdomain contains some oce points -> no need to test the mask used in maxloc
+            llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0 ) == 1._wp        ! define only the inner domain
+            iloc(1:2,1) = MAXLOC( ABS( ssh(:,:,         Kmm)), mask = llmsk(:,:,1) )
+            llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
+            iloc(1:3,2) = MAXLOC( ABS(  uu(:,:,:,       Kmm)), mask = llmsk(:,:,:) )
+            llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
+            iloc(1:3,3) = MINLOC(       ts(:,:,:,jp_sal,Kmm) , mask = llmsk(:,:,:) )
+            iloc(1:3,4) = MAXLOC(       ts(:,:,:,jp_sal,Kmm) , mask = llmsk(:,:,:) )
+            DO ji = 1, jptst   ! local domain indices ==> global domain indices, excluding halos
+               iloc(1:2,ji) = (/ mig0(iloc(1,ji)), mjg0(iloc(2,ji)) /)
+            END DO
+            iareamin(:) = narea   ;   iareamax(:) = narea   ;   iareasum(:) = 1         ! this is local information
+         ENDIF
+         !
+         WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests'
+         CALL wrt_line( ctmp2, kt, '|ssh| max', CASTSP(zmax(1)), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) )
+         CALL wrt_line( ctmp3, kt, '|U|   max', CASTSP(zmax(2)), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) )
+         CALL wrt_line( ctmp4, kt, 'Sal   min', CASTSP(zmax(3)), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) )
+         CALL wrt_line( ctmp5, kt, 'Sal   max', CASTSP(zmax(4)), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) )
+         IF( Agrif_Root() ) THEN
+            WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort* files'
+         ELSE
+            WRITE(ctmp6,*) '      ===> output of last computed fields in '//TRIM(Agrif_CFixed())//'_output.abort* files'
+         ENDIF
+         !
+         CALL dia_wri_state( Kmm, 'output.abort' )     ! create an output.abort file
+         !
+         IF( ll_colruns .OR. jpnij == 1 ) THEN   ! all processes synchronized -> use lwp to print in opened ocean.output files
+            IF(lwp) THEN   ;   CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 )
+            ELSE           ;   nstop = MAX(1, nstop)   ! make sure nstop > 0 (automatically done when calling ctl_stop)
+            ENDIF
+         ELSE                                    ! only mpi subdomains with errors are here -> STOP now
+            CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 )
+         ENDIF
+         !
+      ENDIF
+      !
+      IF( nstop > 0 ) THEN                                                  ! an error was detected and we did not abort yet...
+         ngrdstop = Agrif_Fixed()                                           ! store which grid got this error
+         IF( .NOT. ll_colruns .AND. jpnij > 1 )   CALL ctl_stop( 'STOP' )   ! we must abort here to avoid MPI deadlock
+      ENDIF
+      !
+9500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16)
+      !
+   END SUBROUTINE stp_ctl
+
+
+   SUBROUTINE wrt_line( cdline, kt, cdprefix, pval, kloc, ksum, kmin, kmax )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE wrt_line  ***
+      !!
+      !! ** Purpose :   write information line
+      !!
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*),      INTENT(  out) ::   cdline
+      CHARACTER(len=*),      INTENT(in   ) ::   cdprefix
+      REAL(wp),              INTENT(in   ) ::   pval
+      INTEGER, DIMENSION(3), INTENT(in   ) ::   kloc
+      INTEGER,               INTENT(in   ) ::   kt, ksum, kmin, kmax
+      !
+      CHARACTER(len=80) ::   clsuff
+      CHARACTER(len=9 ) ::   clkt, clsum, clmin, clmax
+      CHARACTER(len=9 ) ::   cli, clj, clk
+      CHARACTER(len=1 ) ::   clfmt
+      CHARACTER(len=4 ) ::   cl4   ! needed to be able to compile with Agrif, I don't know why
+      INTEGER           ::   ifmtk
+      !!----------------------------------------------------------------------
+      WRITE(clkt , '(i9)') kt
+      
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpnij  ,wp))) + 1     ! how many digits to we need to write ? (we decide max = 9)
+      !!! WRITE(clsum, '(i'//clfmt//')') ksum                   ! this is creating a compilation error with AGRIF
+      cl4 = '(i'//clfmt//')'   ;   WRITE(clsum, cl4) ksum
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1    ! how many digits to we need to write ? (we decide max = 9)
+      cl4 = '(i'//clfmt//')'   ;   WRITE(clmin, cl4) kmin-1
+                                   WRITE(clmax, cl4) kmax-1
+      !
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpiglo,wp))) + 1      ! how many digits to we need to write jpiglo? (we decide max = 9)
+      cl4 = '(i'//clfmt//')'   ;   WRITE(cli, cl4) kloc(1)      ! this is ok with AGRIF
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpjglo,wp))) + 1      ! how many digits to we need to write jpjglo? (we decide max = 9)
+      cl4 = '(i'//clfmt//')'   ;   WRITE(clj, cl4) kloc(2)      ! this is ok with AGRIF
+      !
+      IF( ksum == 1 ) THEN   ;   WRITE(clsuff,9100) TRIM(clmin)
+      ELSE                   ;   WRITE(clsuff,9200) TRIM(clsum), TRIM(clmin), TRIM(clmax)
+      ENDIF
+      IF(kloc(3) == 0) THEN
+         ifmtk = INT(LOG10(REAL(jpk,wp))) + 1                   ! how many digits to we need to write jpk? (we decide max = 9)
+         clk = REPEAT(' ', ifmtk)                               ! create the equivalent in blank string
+         WRITE(cdline,9300) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj), clk(1:ifmtk), TRIM(clsuff)
+      ELSE
+         WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpk,wp))) + 1      ! how many digits to we need to write jpk? (we decide max = 9)
+         !!! WRITE(clk, '(i'//clfmt//')') kloc(3)               ! this is creating a compilation error with AGRIF
+         cl4 = '(i'//clfmt//')'   ;   WRITE(clk, cl4) kloc(3)   ! this is ok with AGRIF
+         WRITE(cdline,9400) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj),    TRIM(clk), TRIM(clsuff)
+      ENDIF
+      !
+9100  FORMAT('MPI rank ', a)
+9200  FORMAT('found in ', a, ' MPI tasks, spread out among ranks ', a, ' to ', a)
+9300  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j   ', a, ' ', a, ' ', a, ' ', a)
+9400  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j k ', a, ' ', a, ' ', a, ' ', a)
+      !
+   END SUBROUTINE wrt_line
+
+
+   !!======================================================================
+END MODULE stpctl
Index: CONFIG/UNIFORM/v6/NEMO_v6.5/SOURCES/trcini.F90
===================================================================
--- CONFIG/UNIFORM/v6/NEMO_v6.5/SOURCES/trcini.F90	(revision 6877)
+++ CONFIG/UNIFORM/v6/NEMO_v6.5/SOURCES/trcini.F90	(revision 6877)
@@ -0,0 +1,313 @@
+MODULE trcini
+   !!======================================================================
+   !!                         ***  MODULE trcini  ***
+   !! TOP :   Manage the passive tracer initialization
+   !!======================================================================
+   !! History :   -   ! 1991-03 (O. Marti)  original code
+   !!            1.0  ! 2005-03 (O. Aumont, A. El Moussaoui) F90
+   !!            2.0  ! 2005-10 (C. Ethe, G. Madec) revised architecture
+   !!            4.0  ! 2011-01 (A. R. Porter, STFC Daresbury) dynamical allocation
+   !!----------------------------------------------------------------------
+#if defined key_top
+   !!----------------------------------------------------------------------
+   !!   'key_top'                                                TOP models
+   !!----------------------------------------------------------------------
+   !!   trc_init  :   Initialization for passive tracer
+   !!   top_alloc :   allocate the TOP arrays
+   !!----------------------------------------------------------------------
+   USE par_trc         ! need jptra, number of passive tracers
+   USE oce_trc         ! shared variables between ocean and passive tracers
+   USE trc             ! passive tracers common variables
+   USE trcnam          ! Namelist read
+   USE daymod          ! calendar manager
+   USE prtctl          ! Print control passive tracers (prt_ctl_init routine)
+   USE trcrst
+   USE lib_mpp         ! distribued memory computing library
+   USE trcice          ! tracers in sea ice
+   USE trcbc           ! generalized Boundary Conditions
+   USE trcais          ! tracers from Antartic Ice Sheet
+   USE trcbdy          ! passive-tracer open boundary conditions
+ 
+   IMPLICIT NONE
+   PRIVATE
+   
+   PUBLIC   trc_init   ! called by opa
+
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
+   !! $Id: trcini.F90 15446 2021-10-26 14:34:38Z cetlod $ 
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+   
+   SUBROUTINE trc_init( Kbb, Kmm, Kaa )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_init  ***
+      !!
+      !! ** Purpose :   Initialization of the passive tracer fields 
+      !!
+      !! ** Method  : - read namelist
+      !!              - control the consistancy 
+      !!              - compute specific initialisations
+      !!              - set initial tracer fields (either read restart 
+      !!                or read data or analytical formulation
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa   ! time level indices
+      !
+      IF( ln_timing )   CALL timing_start('trc_init')
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'trc_init : initial set up of the passive tracers'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~'
+      !
+      CALL trc_nam       ! read passive tracers namelists
+      CALL top_alloc()   ! allocate TOP arrays
+
+      !
+      IF(.NOT.ln_trcdta )   ln_trc_ini(:) = .FALSE.
+      !
+      IF(lwp) WRITE(numout,*)
+      IF( ln_rsttr .AND. .NOT. l_offline ) CALL trc_rst_cal( nit000, 'READ' )   ! calendar
+      IF(lwp) WRITE(numout,*)
+      !
+      CALL trc_ini_sms( Kmm )   ! SMS
+      CALL trc_ini_trp          ! passive tracers transport
+      CALL trc_ice_ini          ! Tracers in sea ice
+      !
+      IF( lwm .AND. sn_cfctl%l_trcstat ) THEN
+         CALL ctl_opn( numstr, 'tracer.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp , narea )
+      ENDIF
+      !
+      CALL trc_ini_state( Kbb, Kmm, Kaa )  !  passive tracers initialisation : from a restart or from clim
+      !
+      CALL trc_ini_inv( Kmm )              ! Inventories
+      !
+      IF( ln_timing )   CALL timing_stop('trc_init')
+      !
+   END SUBROUTINE trc_init
+
+
+   SUBROUTINE trc_ini_inv( Kmm )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_ini_stat  ***
+      !! ** Purpose :      passive tracers inventories at initialsation phase
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   Kmm    ! time level index
+      INTEGER             ::  jk, jn  ! dummy loop indices
+      CHARACTER (len=25) :: charout
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra) :: zzmsk
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'trc_ini_inv : initial passive tracers inventories'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      !
+      !                          ! masked grid volume
+      DO jk = 1, jpk
+         cvol(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
+      END DO
+      !                          ! total volume of the ocean 
+      areatot = glob_sum( 'trcini', cvol(:,:,:) )
+      !
+      trai(:) = 0._wp            ! initial content of all tracers
+      DO jn = 1, jptra
+         trai(jn) = trai(jn) + glob_sum( 'trcini', tr(:,:,:,jn,Kmm) * cvol(:,:,:)   )
+      END DO
+
+      IF(lwp) THEN               ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) '   ==>>>   Total number of passive tracer jptra = ', jptra
+         WRITE(numout,*) '           Total volume of ocean                = ', areatot
+         WRITE(numout,*) '           Total inital content of all tracers '
+         WRITE(numout,*)
+         DO jn = 1, jptra
+            WRITE(numout,9000) jn, TRIM( ctrcnm(jn) ), trai(jn)
+         ENDDO
+         WRITE(numout,*)
+      ENDIF
+      IF(lwp) WRITE(numout,*)
+      IF(sn_cfctl%l_prttrc) THEN            ! print mean trends (used for debugging)
+         CALL prt_ctl_init( 'top', jptra )
+         WRITE(charout, FMT="('ini ')")
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl( tab4d_1=tr(:,:,:,:,Kmm), mask1=tmask, clinfo=ctrcnm )
+      ENDIF
+9000  FORMAT('      tracer nb : ',i2,'      name :',a10,'      initial content :',e18.10)
+      !
+   END SUBROUTINE trc_ini_inv
+
+
+   SUBROUTINE trc_ini_sms( Kmm )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_ini_sms  ***
+      !! ** Purpose :   SMS initialisation
+      !!----------------------------------------------------------------------
+      USE trcini_pisces  ! PISCES   initialisation
+      USE trcini_cfc     ! CFC      initialisation
+      USE trcini_c14     ! C14  initialisation
+      USE trcini_age     ! age initialisation
+      USE trcini_my_trc  ! MY_TRC   initialisation
+      !
+      INTEGER, INTENT(in) ::   Kmm ! time level indices
+      INTEGER :: jn
+      !!----------------------------------------------------------------------
+      !
+      ! Pass sn_tracer fields to specialized arrays 
+      DO jn = 1, jp_bgc
+         ctrcnm    (jn) = TRIM( sn_tracer(jn)%clsname )
+         ctrcln    (jn) = TRIM( sn_tracer(jn)%cllname )
+         ctrcun    (jn) = TRIM( sn_tracer(jn)%clunit  )
+         ln_trc_ini(jn) =       sn_tracer(jn)%llinit
+         ln_trc_sbc(jn) =       sn_tracer(jn)%llsbc
+         ln_trc_cbc(jn) =       sn_tracer(jn)%llcbc
+         ln_trc_obc(jn) =       sn_tracer(jn)%llobc
+         ln_trc_ais(jn) =       sn_tracer(jn)%llais
+      END DO
+      !
+      IF( .NOT.ln_trcbc ) THEN
+         DO jn = 1, jp_bgc
+            ln_trc_sbc(jn) = .FALSE.
+            ln_trc_cbc(jn) = .FALSE.
+            ln_trc_obc(jn) = .FALSE.
+         END DO
+      ENDIF
+     
+      lltrcbc = ( COUNT(ln_trc_sbc) + COUNT(ln_trc_obc) + COUNT(ln_trc_cbc) ) > 0 
+      !    
+      IF( ln_pisces      )   CALL trc_ini_pisces( Kmm )     !  PISCES model
+      IF( ln_my_trc      )   CALL trc_ini_my_trc( Kmm )     !  MY_TRC model
+      IF( ll_cfc         )   CALL trc_ini_cfc   ( Kmm )     !  CFC's
+      IF( ln_c14         )   CALL trc_ini_c14   ( Kmm )     !  C14 model
+      IF( ln_age         )   CALL trc_ini_age   ( Kmm )     !  AGE
+      !
+      IF(lwp) THEN                   ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'trc_init_sms : Summary for selected passive tracers'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '    ID     NAME     INI  SBC  CBC  OBC  AIS'
+         DO jn = 1, jptra
+            WRITE(numout,9001) jn, TRIM(ctrcnm(jn)), ln_trc_ini(jn),ln_trc_sbc(jn),ln_trc_cbc(jn),ln_trc_obc(jn),ln_trc_ais(jn)
+         END DO
+      ENDIF
+      IF( lwp .AND. ln_trcbc .AND. lltrcbc ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) ' Applying tracer boundary conditions '
+      ENDIF
+      !
+      IF( lwp .AND. ln_trcais ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) ' Applying tracer from Antarctic Ice Sheet '
+      ENDIF
+     
+9001  FORMAT(3x,i3,1x,a10,3x,l2,3x,l2,3x,l2,3x,l2,3x,l2)
+      !
+   END SUBROUTINE trc_ini_sms
+
+
+   SUBROUTINE trc_ini_trp
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_ini_trp  ***
+      !!
+      !! ** Purpose :   Allocate all the dynamic arrays of the OCE modules
+      !!----------------------------------------------------------------------
+      USE trcdmp , ONLY:  trc_dmp_ini
+      USE trcadv , ONLY:  trc_adv_ini
+      USE trcldf , ONLY:  trc_ldf_ini
+      USE trcrad , ONLY:  trc_rad_ini
+      USE trcsink, ONLY:  trc_sink_ini
+      !
+      INTEGER :: ierr
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_trcdmp )  CALL  trc_dmp_ini          ! damping
+                       CALL  trc_adv_ini          ! advection
+                       CALL  trc_ldf_ini          ! lateral diffusion
+                       !                          ! vertical diffusion: always implicit time stepping scheme
+                       CALL  trc_rad_ini          ! positivity of passive tracers 
+                       CALL  trc_sink_ini         ! Vertical sedimentation of particles
+      !
+   END SUBROUTINE trc_ini_trp
+
+
+   SUBROUTINE trc_ini_state( Kbb, Kmm, Kaa )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_ini_state ***
+      !! ** Purpose :          Initialisation of passive tracer concentration 
+      !!----------------------------------------------------------------------
+      USE trcrst          ! passive tracers restart
+      USE trcdta          ! initialisation from files
+      !
+      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa   ! time level index
+      INTEGER             :: jn, jl          ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_trcdta )   CALL trc_dta_ini( jptra )           ! set initial tracers values
+      !
+      IF( ln_rsttr ) THEN              ! restart from a file
+        !
+        CALL trc_rst_read( Kbb, Kmm )
+        !
+      ELSE                             ! Initialisation of tracer from a file that may also be used for damping
+        IF( ln_trcdta .AND. nb_trcdta > 0 ) THEN
+            ! update passive tracers arrays with input data read from file
+            DO jn = 1, jptra
+               IF( ln_trc_ini(jn) ) THEN
+                  jl = n_trc_index(jn) 
+                  CALL trc_dta( nit000, jl, tr(:,:,:,jn,Kmm) )
+               ENDIF
+            END DO
+            !
+        ENDIF
+        !
+        tr(:,:,:,:,Kbb) = tr(:,:,:,:,Kmm)
+        ! 
+      ENDIF
+      !
+      tr(:,:,:,:,Kaa) = 0._wp
+      !
+      IF( ln_bdy ) CALL trc_bdy_ini( jptra )
+      IF( ln_trcbc .AND. lltrcbc ) CALL trc_bc_ini( jptra, Kmm  )            ! set tracers Boundary Conditions
+      !
+      IF( ln_trcais ) CALL trc_ais_ini   ! set tracers from Antarctic Ice Sheet
+      !                                                         ! Partial top/bottom cell: GRADh(tr(Kmm))
+   END SUBROUTINE trc_ini_state
+
+
+   SUBROUTINE top_alloc
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE top_alloc  ***
+      !!
+      !! ** Purpose :   Allocate all the dynamic arrays of the OCE modules
+      !!----------------------------------------------------------------------
+      USE trc           , ONLY:   trc_alloc
+      USE trdtrc_oce    , ONLY:   trd_trc_oce_alloc
+#if defined key_trdmxl_trc 
+      USE trdmxl_trc    , ONLY:   trd_mxl_trc_alloc
+#endif
+      !
+      INTEGER ::   ierr   ! local integer
+      !!----------------------------------------------------------------------
+      !
+      ierr =        trc_alloc()
+      ierr = ierr + trd_trc_oce_alloc()
+#if defined key_trdmxl_trc 
+      ierr = ierr + trd_mxl_trc_alloc()
+#endif
+      !
+      CALL mpp_sum( 'trcini', ierr )
+      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'top_alloc : unable to allocate standard ocean arrays' )
+      !
+   END SUBROUTINE top_alloc
+
+#else
+   !!----------------------------------------------------------------------
+   !!  Empty module :                                     No passive tracer
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_init                      ! Dummy routine   
+   END SUBROUTINE trc_init
+#endif
+
+   !!======================================================================
+END MODULE trcini
Index: CONFIG/UNIFORM/v6/NEMO_v6.5/SOURCES/trcstp.F90
===================================================================
--- CONFIG/UNIFORM/v6/NEMO_v6.5/SOURCES/trcstp.F90	(revision 6877)
+++ CONFIG/UNIFORM/v6/NEMO_v6.5/SOURCES/trcstp.F90	(revision 6877)
@@ -0,0 +1,291 @@
+MODULE trcstp
+   !!======================================================================
+   !!                       ***  MODULE trcstp  ***
+   !! Time-stepping    : time loop of opa for passive tracer
+   !!======================================================================
+   !! History :  1.0  !  2004-03  (C. Ethe)  Original
+   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme
+   !!----------------------------------------------------------------------
+#if defined key_top
+   !!----------------------------------------------------------------------
+   !!   trc_stp       : passive tracer system time-stepping
+   !!----------------------------------------------------------------------
+   USE par_trc        ! need jptra, number of passive tracers
+   USE oce_trc        ! ocean dynamics and active tracers variables
+   USE sbc_oce
+   USE trc
+   USE trcbc           ! Tracers boundary condtions          ( trc_bc routine)
+   USE trcais          ! Antarctic Ice Sheet tracers         (trc_ais routine)
+   USE trctrp         ! passive tracers transport
+   USE trcsms         ! passive tracers sources and sinks
+   USE trcwri
+   USE trcrst
+   USE trdtrc_oce
+   USE trdmxl_trc
+   USE sms_pisces,  ONLY : ln_check_mass
+   !
+   USE prtctl         ! Print control for debbuging
+   USE iom            !
+   USE in_out_manager !
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   trc_stp    ! called by step
+
+   LOGICAL  ::   llnew                   ! ???
+   REAL(wp) ::   rdt_sampl               ! ???
+   INTEGER  ::   nb_rec_per_day, ktdcy   ! ???
+   REAL(wp) ::   rsecfst, rseclast       ! ???
+   REAL(wp), DIMENSION(:,:,:), SAVE, ALLOCATABLE ::   qsr_arr   ! save qsr during TOP time-step
+
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
+   !! $Id: trcstp.F90 15446 2021-10-26 14:34:38Z cetlod $ 
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE trc_stp( kt, Kbb, Kmm, Krhs, Kaa )
+      !!-------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_stp  ***
+      !!                      
+      !! ** Purpose :   Time loop of opa for passive tracer
+      !! 
+      !! ** Method  :   Compute the passive tracers trends 
+      !!                Update the passive tracers
+      !!-------------------------------------------------------------------
+      INTEGER, INTENT( in ) :: kt                  ! ocean time-step index
+      INTEGER, INTENT( in ) :: Kbb, Kmm, Krhs, Kaa ! time level indices
+      !
+      INTEGER ::   jk, jn   ! dummy loop indices
+      INTEGER ::   ibb      ! local time-level index
+      REAL(wp)::   ztrai    ! local scalar
+      LOGICAL ::   ll_trcstat, ll_trcpis ! local logical
+      CHARACTER (len=25) ::   charout   !
+      !!-------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('trc_stp')
+      !
+      ibb = Kbb                     ! default "before" time-level index
+      IF( l_1st_euler .OR. ln_top_euler ) THEN     ! at nittrc000
+         rDt_trc =  rn_Dt           ! = rn_Dt (use or restarting with Euler time stepping)
+         ibb = Kmm                  ! time-level index used to substitute the "before" with the "now" time level
+      ELSEIF( kt <= nittrc000 + 1 ) THEN                                     ! at nittrc000 or nittrc000+1 
+         rDt_trc = 2. * rn_Dt       ! = 2 rn_Dt (leapfrog) 
+      ENDIF
+      !
+      ll_trcstat  = ( sn_cfctl%l_trcstat ) .AND. &
+     &              ( ( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend ) )
+
+      IF( kt == nittrc000 )                      CALL trc_stp_ctl   ! control 
+      IF( kt == nittrc000 .AND. lk_trdmxl_trc )  CALL trd_mxl_trc_init    ! trends: Mixed-layer
+      !
+      IF( .NOT.ln_linssh ) THEN                                           ! update ocean volume due to ssh temporal evolution
+         DO jk = 1, jpk
+            cvol(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
+         END DO
+         IF( ln_pisces )  THEN
+            IF ( iom_use( "pno3tot" ) .OR. iom_use( "ppo4tot" ) .OR. iom_use( "psiltot" ) &
+               &  .OR. iom_use( "palktot" ) .OR. iom_use( "pfertot" ) ) &
+               & ll_trcpis = .TRUE.
+         ELSE
+            ll_trcpis = .FALSE.
+         ENDIF
+         IF ( ll_trcstat .OR. kt == nitrst .OR. ( ln_check_mass .AND. kt == nitend ) .OR. ll_trcpis ) &
+            &     areatot = glob_sum( 'trcstp', cvol(:,:,:) )
+      ENDIF
+      !
+      IF( l_trcdm2dc )   CALL trc_mean_qsr( kt )
+      !    
+      !
+      IF(sn_cfctl%l_prttrc) THEN
+         WRITE(charout,FMT="('kt =', I4,'  d/m/y =',I2,I2,I4)") kt, nday, nmonth, nyear
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+      ENDIF
+      !
+      tr(:,:,:,:,Krhs) = 0._wp
+      !
+      CALL trc_rst_opn  ( kt )                            ! Open tracer restart file 
+      IF( lrst_trc )  CALL trc_rst_cal  ( kt, 'WRITE' )   ! calendar
+      CALL trc_wri      ( kt,      Kmm            )       ! output of passive tracers with iom I/O manager
+      IF( ln_trcbc .AND. lltrcbc )  CALL trc_bc ( kt, Kmm, tr, Krhs )      ! tracers: surface and lateral Boundary Conditions
+      IF( ln_trcais )               CALL trc_ais( kt, Kmm, tr, Krhs )      ! tracers from Antarctic Ice Sheet (icb, isf)
+      CALL trc_sms      ( kt, ibb, Kmm, Krhs      )       ! tracers: sinks and sources
+#if ! defined key_sed_off
+      CALL trc_trp      ( kt, ibb, Kmm, Krhs, Kaa )       ! transport of passive tracers
+#endif
+           !
+           ! Note passive tracers have been time-filtered in trc_trp but the time level
+           ! indices will not be swapped until after tra_atf/dyn_atf/ssh_atf in stp. Subsequent calls here
+           ! anticipate this update which will be: Nrhs= Nbb ; Nbb = Nnn ; Nnn = Naa ; Naa = Nrhs
+           ! and use the filtered levels explicitly.
+           !
+      IF( kt == nittrc000 ) THEN
+         CALL iom_close( numrtr )                         ! close input tracer restart file
+         IF(lrxios) CALL iom_context_finalize(      cr_toprst_cxt          )
+         IF(lwm) CALL FLUSH( numont )                     ! flush namelist output
+      ENDIF
+      IF( lk_trdmxl_trc  )      CALL trd_mxl_trc  ( kt,      Kaa       )       ! trends: Mixed-layer
+      !
+      IF( ln_top_euler ) THEN 
+         ! For Euler timestepping for TOP we need to copy the "after" to the "now" fields 
+         ! here then after the (leapfrog) swapping of the time-level indices in OCE/step.F90 we have 
+         ! "before" fields = "now" fields.
+         tr(:,:,:,:,Kmm) = tr(:,:,:,:,Kaa)
+      ENDIF
+      !
+      IF( lrst_trc ) CALL trc_rst_wri( kt, Kmm, Kaa, ibb )   ! write tracer restart file
+      !
+      IF (ll_trcstat) THEN
+         ztrai = 0._wp                                                   !  content of all tracers
+         DO jn = 1, jptra
+            ztrai = ztrai + glob_sum( 'trcstp', tr(:,:,:,jn,Kaa) * cvol(:,:,:)   )
+         END DO
+         IF( lwm ) WRITE(numstr,9300) kt,  ztrai / areatot
+      ENDIF
+9300  FORMAT(i10,D23.16)
+      !
+      IF( ln_timing )   CALL timing_stop('trc_stp')
+      !
+   END SUBROUTINE trc_stp
+
+
+   SUBROUTINE trc_stp_ctl
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_stp_ctl  ***
+      !!----------------------------------------------------------------------
+      !
+      ! Define logical parameter ton control dirunal cycle in TOP
+      l_trcdm2dc = ( ln_trcdc2dm .AND. .NOT. ln_dm2dc  ) 
+      !
+      IF( l_trcdm2dc .AND. lwp )   CALL ctl_warn( 'Coupling with passive tracers and used of diurnal cycle.',   &
+         &                           'Computation of a daily mean shortwave for some biogeochemical models ' )
+      !
+   END SUBROUTINE trc_stp_ctl
+
+
+   SUBROUTINE trc_mean_qsr( kt )
+      !!----------------------------------------------------------------------
+      !!             ***  ROUTINE trc_mean_qsr  ***
+      !!
+      !! ** Purpose :  Compute daily mean qsr for biogeochemical model in case
+      !!               of diurnal cycle
+      !!
+      !! ** Method  : store in TOP the qsr every hour ( or every time-step if the latter 
+      !!              is greater than 1 hour ) and then, compute the  mean with 
+      !!              a moving average over 24 hours. 
+      !!              In coupled mode, the sampling is done at every coupling frequency 
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      !
+      INTEGER  ::   jn   ! dummy loop indices
+      REAL(wp) ::   zkt, zrec     ! local scalars
+      CHARACTER(len=1) ::   cl1   ! 1 character
+      CHARACTER(len=2) ::   cl2   ! 2 characters
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('trc_mean_qsr')
+      !
+      IF( kt == nittrc000 ) THEN
+         !
+         rdt_sampl = REAL( ncpl_qsr_freq )
+         nb_rec_per_day = INT( rday / ncpl_qsr_freq )
+         !
+         IF(lwp) THEN
+            WRITE(numout,*) 
+            WRITE(numout,*) ' Sampling frequency dt = ', rdt_sampl, 's','   Number of sampling per day  nrec = ', nb_rec_per_day
+            WRITE(numout,*) 
+         ENDIF
+         !
+         ALLOCATE( qsr_arr(jpi,jpj,nb_rec_per_day ) )
+         !
+         !                                            !* Restart: read in restart file
+         IF( ln_rsttr .AND. nn_rsttr /= 0 .AND. iom_varid( numrtr, 'qsr_mean' , ldstop = .FALSE. ) > 0  &
+           &                              .AND. iom_varid( numrtr, 'qsr_arr_1', ldstop = .FALSE. ) > 0  &
+           &                              .AND. iom_varid( numrtr, 'ktdcy'    , ldstop = .FALSE. ) > 0  &
+           &                              .AND. iom_varid( numrtr, 'nrdcy'    , ldstop = .FALSE. ) > 0  ) THEN
+            CALL iom_get( numrtr, 'ktdcy', zkt )  
+            rsecfst = INT( zkt ) * rn_Dt
+            IF(lwp) WRITE(numout,*) 'trc_qsr_mean:   qsr_mean read in the restart file at time-step rsecfst =', rsecfst, ' s '
+            CALL iom_get( numrtr, jpdom_auto, 'qsr_mean', qsr_mean )   !  A mean of qsr
+            CALL iom_get( numrtr, 'nrdcy', zrec )   !  Number of record per days
+            IF( INT( zrec ) == nb_rec_per_day ) THEN
+               DO jn = 1, nb_rec_per_day 
+                  IF( jn <= 9 )  THEN
+                    WRITE(cl1,'(i1)') jn
+                    CALL iom_get( numrtr, jpdom_auto, 'qsr_arr_'//cl1, qsr_arr(:,:,jn) )   !  A mean of qsr
+                  ELSE
+                    WRITE(cl2,'(i2.2)') jn
+                    CALL iom_get( numrtr, jpdom_auto, 'qsr_arr_'//cl2, qsr_arr(:,:,jn) )   !  A mean of qsr
+                  ENDIF
+              END DO
+            ELSE
+               DO jn = 1, nb_rec_per_day
+                  qsr_arr(:,:,jn) = qsr_mean(:,:)
+               ENDDO
+            ENDIF
+         ELSE                                         !* no restart: set from nit000 values
+            IF(lwp) WRITE(numout,*) 'trc_qsr_mean:   qsr_mean set to nit000 values'
+            rsecfst  = kt * rn_Dt
+            !
+            qsr_mean(:,:) = qsr(:,:)
+            DO jn = 1, nb_rec_per_day
+               qsr_arr(:,:,jn) = qsr_mean(:,:)
+            END DO
+         ENDIF
+         !
+      ENDIF
+      !
+      rseclast = kt * rn_Dt
+      !
+      llnew   = ( rseclast - rsecfst ) .ge.  rdt_sampl    !   new shortwave to store
+      IF( llnew ) THEN
+          ktdcy = kt
+          IF( lwp .AND. kt < nittrc000 + 100 ) WRITE(numout,*) ' New shortwave to sample for TOP at time kt = ', ktdcy, &
+             &                      ' time = ', rseclast/3600.,'hours '
+          rsecfst = rseclast
+          DO jn = 1, nb_rec_per_day - 1
+             qsr_arr(:,:,jn) = qsr_arr(:,:,jn+1)
+          ENDDO
+          qsr_arr (:,:,nb_rec_per_day) = qsr(:,:)
+          qsr_mean(:,:                ) = SUM( qsr_arr(:,:,:), 3 ) / nb_rec_per_day
+      ENDIF
+      !
+      IF( lrst_trc ) THEN    !* Write the mean of qsr in restart file 
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'trc_mean_qsr : write qsr_mean in restart file  kt =', kt
+         IF(lwp) WRITE(numout,*) '~~~~~~~'
+         zkt  = REAL( ktdcy, wp )
+         zrec = REAL( nb_rec_per_day, wp )
+         CALL iom_rstput( kt, nitrst, numrtw, 'ktdcy', zkt  )
+         CALL iom_rstput( kt, nitrst, numrtw, 'nrdcy', zrec )
+          DO jn = 1, nb_rec_per_day 
+             IF( jn <= 9 )  THEN
+               WRITE(cl1,'(i1)') jn
+               CALL iom_rstput( kt, nitrst, numrtw, 'qsr_arr_'//cl1, qsr_arr(:,:,jn) )
+             ELSE
+               WRITE(cl2,'(i2.2)') jn
+               CALL iom_rstput( kt, nitrst, numrtw, 'qsr_arr_'//cl2, qsr_arr(:,:,jn) )
+             ENDIF
+         END DO
+         CALL iom_rstput( kt, nitrst, numrtw, 'qsr_mean', qsr_mean(:,:) )
+      ENDIF
+      !
+      IF( ln_timing )   CALL timing_stop('trc_mean_qsr')
+      !
+   END SUBROUTINE trc_mean_qsr
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default key                                     NO passive tracers
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_stp( kt )        ! Empty routine
+      WRITE(*,*) 'trc_stp: You should not have seen this print! error?', kt
+   END SUBROUTINE trc_stp
+#endif
+
+   !!======================================================================
+END MODULE trcstp
Index: CONFIG/UNIFORM/v6/NEMO_v6.5/SOURCES/trctrp.F90
===================================================================
--- CONFIG/UNIFORM/v6/NEMO_v6.5/SOURCES/trctrp.F90	(revision 6877)
+++ CONFIG/UNIFORM/v6/NEMO_v6.5/SOURCES/trctrp.F90	(revision 6877)
@@ -0,0 +1,122 @@
+MODULE trctrp
+   !!======================================================================
+   !!                       ***  MODULE trctrp  ***
+   !! Ocean Physics    : manage the passive tracer transport
+   !!======================================================================
+   !! History :   1.0  !  2004-03 (C. Ethe) Original code
+   !!             3.3  !  2010-07 (C. Ethe) Merge TRA-TRC
+   !!----------------------------------------------------------------------
+#if defined key_top
+   !!----------------------------------------------------------------------
+   !!   'key_top'                                                TOP models
+   !!----------------------------------------------------------------------
+   !!   trc_trp        : passive tracer transport
+   !!----------------------------------------------------------------------
+   USE par_trc         ! need jptra, number of passive tracers
+   USE oce_trc         ! ocean dynamics and active tracers variables
+   USE trc             ! ocean passive tracers variables 
+   USE trcbbl          ! bottom boundary layer               (trc_bbl routine)
+   USE trcdmp          ! internal damping                    (trc_dmp routine)
+   USE trcldf          ! lateral mixing                      (trc_ldf routine)
+   USE trcadv          ! advection                           (trc_adv routine)
+   USE trczdf          ! vertical diffusion                  (trc_zdf routine)
+   USE trcatf          ! time filtering                      (trc_atf routine)
+   USE trcrad          ! positivity                          (trc_rad routine)
+   USE trcsbc          ! surface boundary condition          (trc_sbc routine)
+   USE trcbc           ! Tracers boundary condtions          ( trc_bc routine)
+   USE trcais          ! Antarctic Ice Sheet tracers         (trc_ais routine)
+   USE zpshde          ! partial step: hor. derivative       (zps_hde routine)
+   USE bdy_oce   , ONLY: ln_bdy
+   USE trcbdy          ! BDY open boundaries
+   USE in_out_manager
+
+#if defined key_agrif
+   USE agrif_top_sponge ! tracers sponges
+#endif
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   trc_trp    ! called by trc_stp
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
+   !! $Id: trctrp.F90 15023 2021-06-18 14:35:25Z gsamson $ 
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE trc_trp( kt, Kbb, Kmm, Krhs, Kaa )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_trp  ***
+      !!                      
+      !! ** Purpose :   Management of passive tracers transport
+      !! 
+      !! ** Method  : - Compute the passive tracers trends 
+      !!              - Update the passive tracers
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) :: kt                  ! ocean time-step index
+      INTEGER, INTENT( in ) :: Kbb, Kmm, Krhs, Kaa ! time level indices (not swapped in this routine)
+      !! ---------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('trc_trp')
+      !
+      IF( .NOT. ln_c1d ) THEN
+         !
+         !                                                         ! Partial top/bottom cell: GRADh( trb )  
+         IF( ln_zps ) THEN
+            IF( ln_isfcav ) THEN ; CALL zps_hde_isf( kt, jptra, tr(:,:,:,:,Kbb), pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )  ! both top & bottom
+            ELSE                 ; CALL zps_hde    ( kt, jptra, tr(:,:,:,:,Kbb), gtru, gtrv )                                      !  only bottom
+            ENDIF
+         ENDIF
+         !
+                                CALL trc_sbc    ( kt,      Kmm, tr, Krhs )      ! surface boundary condition
+         IF( ln_trabbl )        CALL trc_bbl    ( kt, Kbb, Kmm, tr, Krhs )      ! advective (and/or diffusive) bottom boundary layer scheme
+         IF( ln_trcdmp )        CALL trc_dmp    ( kt, Kbb, Kmm, tr, Krhs )      ! internal damping trends
+         IF( ln_bdy )           CALL trc_bdy_dmp( kt, Kbb,      Krhs )      ! BDY damping trends
+#if defined key_agrif
+         IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_trc       ! tracers sponge
+#endif
+                                CALL trc_adv    ( kt, Kbb, Kmm, tr, Krhs )      ! horizontal & vertical advection 
+                                CALL trc_ldf    ( kt, Kbb, Kmm,       tr, Krhs )  ! lateral mixing
+                                CALL trc_zdf    ( kt, Kbb, Kmm, Krhs, tr, Kaa  )  ! vert. mixing & after tracer	==> after
+                                CALL trc_atf    ( kt, Kbb, Kmm, Kaa , tr )        ! time filtering of "now" tracer fields    
+         !
+         ! Subsequent calls use the filtered values: Kmm and Kaa 
+         ! These are used explicitly here since time levels will not be swapped until after tra_atf/dyn_atf/ssh_atf in stp
+         !
+         IF( ln_trcrad )        CALL trc_rad    ( kt, Kmm, Kaa, tr       )    ! Correct artificial negative concentrations
+         IF( ln_trcdmp_clo )    CALL trc_dmp_clo( kt, Kmm, Kaa )              ! internal damping trends on closed seas only
+
+         !
+      ELSE                                               ! 1D vertical configuration
+                                CALL trc_sbc( kt,      Kmm,       tr, Krhs )  ! surface boundary condition
+         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_atf( kt, Kbb, Kmm, Kaa , tr )        ! time filtering of "now" tracer fields
+         !
+         ! Subsequent calls use the filtered values: Kmm and Kaa 
+         ! These are used explicitly here since time levels will not be swapped until after tra_atf/dyn_atf/ssh_atf in stp
+         !
+         IF( ln_trcrad )       CALL trc_rad( kt, Kmm, Kaa, tr       )  ! Correct artificial negative concentrations
+         !
+      END IF
+      !
+      IF( ln_timing )   CALL timing_stop('trc_trp')
+      !
+   END SUBROUTINE trc_trp
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module :                                        No TOP models
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_trp( kt )              ! Empty routine
+      INTEGER, INTENT(in) ::   kt
+      WRITE(*,*) 'trc_trp: You should not have seen this print! error?', kt
+   END SUBROUTINE trc_trp
+#endif
+   
+   !!======================================================================
+END MODULE trctrp
