Tracer conservation for configurations with forward Euler passive-tracer time stepping and time-varying grid-cell volumes

Context

The reference configuration ORCA2_ICE_PISCES of NEMO version 4.2 does not conserve tracer inventories. By selecting the modified leapfrog passive-tracer time-stepping option, however, tracer conservation can be restored. Further, the built-in computation of basic passive-tracer statistics (module trcrst, subroutine trc_rst_stat) is inaccurate when grid-cell volumes vary in time, which hampers rapid validation of tracer-inventory conservation.

Analysis

The lack of tracer conservation found in the ORCA2_ICE_PISCES reference configuration can be demonstrated in a simpler configuration by enabling the nonlinear free surface in the GYRE_PISCES reference configuration,

$ ./makenemo -r GYRE_PISCES -n TEST del_key key_linssh
$ cd cfgs/TEST/EXP00/
$ sed -i -e 's/ln_linssh.*/ln_linssh=.false./' -e 's/ln_hpg_zco.*/ln_hpg_sco=.true./' namelist_cfg

by replacing the PISCES option with the transport of a passive tracer with initial uniform distribution and without external sources and sinks,

$ sed -i -e 's/jp_bgc.*/jp_bgc=1/' -e 's/ln_pisces.*/ln_pisces=.false./' -e 's/ln_my_trc.*/ln_my_trc=.true./' namelist_top_cfg

and by activating forward Euler time stepping for passive tracers,

$ sed -i -e '/&namtrc_run/aln_top_euler=.true./' namelist_top_cfg

The output from the corresponding model run reports a substantial drift in the tracer inventory:

$ grep -e 'DET.*drift' ocean.output | cut -d ':' -f 7
 -0.4331065867E-07 %

When selecting the modified leapfrog time-stepping scheme instead,

$ sed -i -e 's/ln_top_euler.*/ln_top_euler=.false./' namelist_top_cfg

the tracer-inventory drift reduces to

$ grep -e 'DET.*drift' ocean.output | cut -d ':' -f 7
 -0.1484324163E-09 %

and, after adjusting the grid-cell volumes used for the computation of tracer statistics at the end of the model run,

--- a/src/TOP/trcrst.F90
+++ b/src/TOP/trcrst.F90
@@ -356,7 +356,7 @@ CONTAINS
       ENDIF
       !
       DO jk = 1, jpk
-         zvol(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Krhs) * tmask(:,:,jk)
+         zvol(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
       END DO
       !
       DO jn = 1, jptra

and subsequent compilation, this drift reduces further to

$ grep -e 'DET.*drift' ocean.output | cut -d ':' -f 7
 -0.2386728137E-11 %

which appears to be within the uncertainty that is expected from using double-precision arithmetic operations. After switching back to forward Euler time stepping,

$ sed -i -e 's/ln_top_euler.*/ln_top_euler=.true./' namelist_top_cfg

a new model run reveals the actual tracer-inventory drift of

$ grep -e 'DET.*drift' ocean.output | cut -d ':' -f 7
 -0.4322051485E-07 %

for this configuration.

Further analysis suggests that the replication of the "now" tracer concentrations at the "before" time level effectively changes the local tracer inventories (due to a difference between the grid-cell volumes at the "before" and "now" time levels) that are used as the tracer state before the trends of the current time step are added.

Fix

It is recommended to implement the forward Euler time stepping independent of the tracer concentrations and vertical scale factors for the "before" time step, which could be achieved by supplying the "now" time-level index as both the "before" and "now" time-level indices to the passive-tracer transport and source subroutines (subroutines trc_trp and trc_sms) in this case:

--- a/src/TOP/trcstp.F90
+++ b/src/TOP/trcstp.F90
@@ -58,6 +58,7 @@ CONTAINS
       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 ! local logical
       CHARACTER (len=25) ::   charout   !
@@ -65,8 +66,10 @@ CONTAINS
       !
       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
@@ -100,9 +103,9 @@ CONTAINS
       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
-      CALL trc_sms      ( kt, Kbb, Kmm, Krhs      )       ! tracers: sinks and sources
+      CALL trc_sms      ( kt, ibb, Kmm, Krhs      )       ! tracers: sinks and sources
 #if ! defined key_sed_off
-      CALL trc_trp      ( kt, Kbb, Kmm, Krhs, Kaa )       ! transport of passive tracers
+      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
@@ -115,7 +118,7 @@ CONTAINS
          IF(lrxios) CALL iom_context_finalize(      cr_toprst_cxt          )
          IF(lwm) CALL FLUSH( numont )                     ! flush namelist output
       ENDIF
-      IF( lrst_trc )            CALL trc_rst_wri  ( kt, Kmm, Kaa, Kbb  )       ! write tracer restart file
+      IF( lrst_trc )            CALL trc_rst_wri  ( kt, Kmm, Kaa, ibb  )       ! write tracer restart file
       IF( lk_trdmxl_trc  )      CALL trd_mxl_trc  ( kt,      Kaa       )       ! trends: Mixed-layer
       !
       IF( ln_top_euler ) THEN 

In addition, the adjustment of the built-in computation of passive-tracer statistics,

--- a/src/TOP/trcrst.F90
+++ b/src/TOP/trcrst.F90
@@ -356,7 +356,7 @@ CONTAINS
       ENDIF
       !
       DO jk = 1, jpk
-         zvol(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Krhs) * tmask(:,:,jk)
+         zvol(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
       END DO
       !
       DO jn = 1, jptra

is also recommended.