Preparation of the NEMO code for 2LTS schemes: changes to NEMO 2019

1. Overview and motivation

The current NEMO code (version 4 and earlier) is coded to use a three-level leapfrog timestepping scheme. As part of the IMMERSE project a new two-level timestepping scheme is being implemented. In order to lay the groundwork for this, a major refactoring of the NEMO code has been undertaken to be included in the 2019 merge. This document gives an overview of the changes.

2. New prognostic variables

The main point of the change is to recast all the main prognostic variables to have a time dimension rather than separate arrays for before, now, and after time levels. The following table gives a non-exhaustive, but representative list of the changes to the main prognostic variables in NEMO. Here (ji,jj,jk) are the three spatial dimensions, jn is the index on tracers, and jt is the new time dimension.

OLD	NEW
tsb/tsn/tsa (ji, jj, jk, jn)	ts(ji, jj, jk, jn, jt)
ub/un/ua (ji, jj, jk)	uu(ji, jj, jk, jt)
wn(ji,jj, jk)	ww(ji, jj, jk)
sshb/sshn/ssha (ji,jj)	ssh(ji, jj, jt)
ub_b/un_b/ua_b (ji, jj)	uu_b(ji, jj, jt)
e3t_0/e3t_b/e3t_n/e3t_a (ji, jj, jk)	e3t_0(ji, jj, jk) / e3t(ji, jj, jk, jt)
hu_0/hu_b/hu_n/hu_a (ji, jj)	hu_0(ji, jj) / hu(ji, jj, jt)
ht_0/ht_n (ji, jj)	ht_0(ji, jj) / ht(ji, jj)

Notes:

• Where only one time level of a field is stored in the old code (eg. wn) we simply rename the variable (wn -> ww) to reflect the new schema.
• The barotropic variables defined on the baroclinic timesteps (ssh, uu_b) are changed but the variables used within the barotropic subcycling (ub_e, un_e etc) are not changed.
• The names of the variables that NEMO writes to and reads from a restart are currently unchanged (eg. tb, tn, ub, un).

3. Time level indices

Time level indices are used to index the time dimension of prognostic variables. The values of these indices can be updated to achieve the updating of the prognostic variables (eg. now -> before) without moving data in memory. The indices are defined in step.F90 and passed through the code in subroutine arguments. The names of the indices are different in step.F90 and in other modules as per the DOCTOR coding norm. They are as follows:

time level	name in step.F90	name in other modules
before	Nbb	Kbb
now / intermediate estimate	Nnn	Kmm
after	Naa	Kaa
right-hand side = trend	Nrhs	Krhs

Notes:

• For leapfrog timestepping Naa=Nrhs always because the same area of memory is used to store the trend and the updated field at different stages. However two different index variables are used to indicate whether that part of the array is currently storing the trend or the updated field.
• The option to use a different timestep for TOP than is used for OCE has been deprecated and TOP uses the same time index variables as OCE. The ln_top_euler option to allow TOP to do a forward Euler timestep throughout the integration is still available but likely to be deprecated when the new NEMO timestepping scheme is implemented.

4. Subroutine interfaces

Time level indices are passed through the code in subroutine argument lists rather than via USE statements. In addition, for the routines to calculate dynamics and tracer trends called from the "stp" routines, the prognostic fields to be modified are passed as arguments. The last time level argument indicates whether the routine returns the trend or the updated "after" field. For example the call to dyn_adv from stp is modified as follows.

old code:

      CALL dyn_adv       ( kstp )

new code:

      CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )

where the final Nrhs argument indicates that this routine is returning a modified trend. The call to the vertical advection routine in the new code is

      CALL dyn_zdf( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )

where the final Naa argument indicates that this routine is returning the "after" field. The prognostic fields passed through argument lists have a leading "p" in the low-level routines as per the DOCTOR coding norm, so for example in dyn_adv, uu and vv are referred to as puu and pvv.

Notes:

• Exceptions had to be made to the passing of time-levels through subroutine arguments for AGRIF routines. The time-levels need to be known in some of the routines in the NST directory but these routines are called by several layers of generated interfaces and passing arguments all the way down the calling chain proved problematic. The less-invasive solution was to allow agrif_oce.F90 to store its own copies of the time-level variables within it's own module. These copies are used within AGRIF routines and need to be set prior to calling the AGRIF routines.

     agrif_oce.F90:
     INTEGER , PUBLIC,  SAVE  ::  Kbb_a, Kmm_a, Krhs_a
                              !: AGRIF module-specific copies of time-level indices
  

These are typically set in nemogcm and step:

   Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs

Or can be set explicitly before a call to an AGRIF routine:

   ./DYN/sshwzv.F90:
   Kbb_a = Kbb; Kmm_a = Kmm; Krhs_a = Kaa; CALL agrif_ssh( kt )

5. Reorganisation of time filtering and time level swapping

In the old version of the code, the Asselin time filter is applied and the time levels swapped in the "nxt" routines: dyn_nxt, tra_nxt, trc_nxt. In the new code, the time level swapping is achieved by swapping the values of the time index variables in stp:

      ! Swap time levels
	      Nrhs = Nbb
	      Nbb = Nnn
	      Nnn = Naa
	      Naa = Nrhs

The "nxt" modules and routines are renamed to reflect the fact that their main function is now just to do the Asselin time filtering, so dynnxt.F90 ? dynatf.F90, dyn_nxt -> dyn_atf etc. Note that other functionality apart from the filtering and the time-level swapping has been retained in the "atf" routines; for example the code to adjust the barotropic part of the solution in dyn_nxt.