New URL for NEMO forge!

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
ticket/0663_MLF – NEMO

Version 40 (modified by mlelod, 13 years ago) (diff)


Last edited Timestamp?

Author : Gurvan Madec & Matthieu Leclair

ticket : #663

Branch : DEV_r1837_MLF

useful commands

* check out of NEMO directory:
svn co svn+ssh://

* commit :
svn ci  -m  "ticket: #663  blah_blah"    list_of_file


Introduce the modified Leap-Frog -- Robert-Asselin filter (LF-RA) to fit the one presented in Leclair and Madec Ocean Modelling (2009). The scheme allows a exact conservation of the heat and salt contents when using variable volume thickness (key_vvl). The non-conservation of the standard LF-RA due to the time diffusion of the forcing through the filter. In the modified scheme the exact conservation is achieved by removing the forcing from the filter. Furthermore, the scheme allows to use a much smaller Asselin filter parameter (10-4 instead of 10-1) decreasing the damping effect of the LF-RA scheme, so that it is actually a quasi second order scheme.

The LF-RA flow chart of step is the one introduce in the trunk after the v3.1, to restore the key_vvl option. Its main characteristic is that the ssh is time-stepped first, followed immediately by the computation of the now vertical velocity (see figure below).

step flowchart for the modified LF-RA scheme
Step routine flow chart, form Leclair and Madec (OM 2009).

This update is part of the LOCEAN.8 task (Update of the ocean physics) due by June 25th 2010.

The 3 main modifications concern :

(1) The forcing terms. They are defined as a mean between time step n-1/2 and n+1/2

modules involved : sbc_oce.F90, sbcmod.F90, trasbc.F90, traqsr.F90, dynzdf....F90

step flowchart for the modified LF-RA scheme
Illustration of forcing integration methods. Old formulation (top) and new formulation (bottom).

At each time a forcing term appears in the 3 time evolution equations solved in NEMO/OPA (ssh, dynamics and tracer). The mean of kt+1/2 and kt+1/2 has to be used in order to remove the largest source of divergence of two consecutive time step. Therefore, the surface module (including the sea-ice!) must now provide both kt+1/2 and kt-1/2 forcing fields, and all the module involving a forcing term have to be modified. Note that TKE equation is solved with a two level time stepping scheme : no need of a before field for its forcing term (i.e. taum)

sbc_oce.F90 :

Add "before" forcing fields (denoted by a '_b' sufix)

All ocean forcing components are involved, i.e. utau_b, vtau_b, qns_b, qsr_b, emp_b,and emps_b

6 other fields are introduced: sbc_hc_n and sbc_hc_b are the surface heat content foring at now and before time steps, sbc_sc_n and sbc_sc_b are the surface salt content forcing at now and before time steps and qsr_hc_n and qsr_hc_b are the 3D penetrative solar heat content forcing at now and before time steps.

sbcmod.F90 :

Add the swap of surface fields and the read/write in the restart file, except for sbc_hc_(n/b), sbc_sc_(n/b) ans qsr_hc_(n/b). Note: when writing restart files, the utau_b, vtau_b, qns_b variables etc... contain utau, vtau, qns etc...

Beginning of sbs routine, introduce the swap, but only for kt >= nit000+1. If no restart, the '_b' fields is known only at the end of the sbc calculation. Therefore we have choosen to regroup the setting of '_b' fields at nit000 at the end of sbc routine.

CAUTION : here the swap must be done at each kt, what ever the nn_fsbc value is. Otherwise the ocean forcing will be wrong.

trasbc.F90 : Change in the heat and freshwater/salt fluxes

The thermohalin fluxes (qns, emp, emps) are now applied as the mean of 2 consecutive time-steps (kt+1/2 and kt-1/2). This involes the following changes : (1) incorporate the time mean of the forcing ; (2) put the IF-ELSE-ENDIF outside the DO-loop for vector optimisation ; (3) remove the key_zco optimisation (increase readability) as zco is less and less used ; (4) correct a bug in vvl case, where salt flux associated with freezing/melting were not taken into account.

Swap and restart of sbc_(h/s)c

CAUTION 1 : bug in the salt flux in key_vvl ===> modif in limsbc_2.F90 (emps becomes a salt flux in vvl,limsbc, limsbc_2 and CICE interface must be modified)

CAUTION 2 : bug in with no solar light penetration (ln_traqsr=.false.) ==> modif in trasbc : qns = qns + qsr.

traqsr.F90 : apply the solar heat flux (qsr) as the mean of 2 consecutive time-steps (kt+1/2 and kt-1/2)

The 3D heat content tendency term of penetrative solar flux is stored in "qsr_hc_n" in all cases. To that purpose, the division by fse3t in the computation of etot3 has been removed. It is done in the next step, when heat content tendency is converted to temperature tendency and added to ta.

CAUTION 1 : the case where the light penetration is provided by the bio-model ( i.e. lk_qsr_bio = ln_qsr_bio =.true. is not yet treated. It probably requires changes in p4zopt.F90 and trcopt.F90. Not completely obvious since this depends on the TOP time-stepping strategy.

CAUTION 2 : the thickness used in the penetration should be the mean of the two consecutive time step !

Swap and restart of qsr_hc

dynzdf_exp.F90 : apply the stress (utau, vtau) as the mean of 2 consecutive time-steps (kt+1/2 and kt-1/2)

dynzdf_imp.F90 : apply the stress (utau, vtau) as the mean of 2 consecutive time-steps (kt+1/2 and kt-1/2)

sshwzv.F90 : apply the freshwater flux (emp) as the mean of 2 consecutive time-steps (kt+1/2 and kt-1/2)

Not related to this ticket: surface weighted averaging for the computation of sshf_a and sshf_b in vvl case (as for sshu and sshv)

(2) The Asselin filter remove the forcing from the filter (tranxt and sshwzv)

sshwzv : filter correction when computing sshb in vvl case: remove atfp * rdt * ( 'volume_forcing_n' - 'volume_forcing_b' ) from volume time filter

tranxt.F90 : Changes for modified time filtering - only in the variable volume case

remove atfp * rdt * ( 'heat_content_forcing_n' - 'heat_content_forcing_b' ) from heat content time filter

remove atfp * rdt * ( 'salt_content_forcing_n' - 'salt_content_forcing_b' ) from heat content time filter

remove atfp * rdt * ( 'volume_forcing_n' - 'volume_forcing_b' ) from volume time filter

N.B. : We choose not to remove the momentum forcing from the Asselin time filter.

(3) local compatibility.

Local compatibility between tracer content and volume budgets leads to a full 3D treatment of fse3(t/u/v)_b

the time filter correction on volume fluxes in only added at the first level of fse3t_b. Then fse3(u/v)_b are deduced via a surface averaging (like ssh(u/v) deduced from sshn).

dynnxt.F90 : Computation of filtered and corrected before scale factors

Computation of fse3t_b with filter correction.

Then computation of ze3(u/v)_f the now filtered scale factor used for the thickness weighted time stepping.

Then storage of zse3(u/v)_f into fse3(u/v)_b.

restart.F90 : write and read fse3t_b in the vvl case

domvvl.F90 : initialise fse3(u/v)_b

Some intermediate variables have been introduced so that the computation of ssh(u/v/f) has also been slightly modified

sshwzv.F90 : Computation of the "diffused" vertical scale factor fse3t_d needed for the computation of Robert-Asselin and Brown & Campana filters: fse3t_d = fse3t_b - 2 fse3t_n + fse3t_a

It is used in tranxt and dynnxt. It wouldn't be necessary for these 2 "little" reasons but it will later be used for the computation of horizontal pressure gradient in the semi implicit and vvl case.

dom_oce.F90 : declarations in the vvl case of e3(t/u/v)_b and e3t_d

(4) The semi-implicit hpg.

In tranxt.F90, for the semi-implicit hydrostatic pressure gradient, the coefficients are changed in order to take into account the Asselin filter parameter

time mean T and S for the semi implicit hpg computation was (ta+2*tn+tb) / 4 . It is now replaced by


with rbcp=(1+atfp)*(1+atfp)*(1-atfp)/4, and atfp is the Asselin filter parameter.

The change extend further the stability limit of the time stepping scheme when ln_dynhpg_imp = T

The semi implicit hpg case is included in the loop (tran_nxt_fix and tra_nxt_vvl). It does not break the loop vectorization.

*_d fields contain time laplacian used for both Asselin and Brown & Campana filters (tran_nxt_fix and tra_nxt_vvl)

In variable volume case: Apply B & C averaging on heat and salt contents then devide resulting heat and salt contents by B & C averaging on volume.

(5) Global conservation diagnostic.

A new routine, diatrb, has been added in order to diagnose the global non-conservation level.

N.B. : Three key elements of the LF-RA scheme have already been implemented in former version of NEMO (v3.0 and v3.1). As a remainder, there are:

Forcing evaluation : the momentum and tracers forcing are evaluated at kt+1/2, no-more at kt. This is done in the fldread routine since v3.0

Flow chart : general structure has changed for key_vvl (variable volume or equivalently non-linear free surface) (v3.1). In particular the ssh equation is solved first together followed immediately by the computation of the now vertical velocity, so that wn is available to solve the tracer and the dynamics equations.

The enhanced vertical diffusion (evd) is no more a source of small divergence of two consecutive time step as the static instability is now tested on both before and now fields. See ticket #401.

Pending issues

  • correct the bug in limsbc, limsbc_2 and CICE interface (nned help from Met Office for the later)
  • Do we need to remove the stress component from the filter in dynnxt? NO
  • in vvl case : thickness weighted tracer in the time mean, do we need to compute and use a time mean thickness in the hpg computation? YES
  • The associated changes in TOP are still unclear. A priori, in PISCES p4zopt.F90 should be modified. Nevertheless, the TOP time stepping strategy seems to be a forward in time scheme. If yes, then the evaluation of the forcing at kt+1/2 is an improvement (if not to say a bug fix!). If not, action is required...
  • Check the consistency with the change planned in the emps field (switch to fsalt, a true salt flux).
  • modification required in limsbc, limsbc_2 and CICE interface : conversion of emps into a salt flux (new name 'fsalt' ?)
  • the modifications in trasbc will interfere with the new runoff scheme (see 2009WP/2009Stream3#Substream3.2:fromNEMOteamshort-term)
  • Caution: Vertical mass transport is only the baroclinic mas transport in the vvl case. only the Eulerian part of vertical motions is taken into account.


Testing could consider (where appropriate) other configurations in addition to NVTK].

NVTK Tested'''YES/NO'''
Other model configurations'''YES/NO'''
Processor configurations tested[ Enter processor configs tested here ]
If adding new functionality please confirm that the
New code doesn't change results when it is switched off
and ''works'' when switched on

(Answering UNSURE is likely to generate further questions from reviewers.)

'Please add further summary details here'

  • Processor configurations tested
  • etc

Bit Comparability

Does this change preserve answers in your tested standard configurations (to the last bit) ?'''YES/NO '''
Does this change bit compare across various processor configurations. (1xM, Nx1 and MxN are recommended)'''YES/NO'''
Is this change expected to preserve answers in all possible model configurations?'''YES/NO'''
Is this change expected to preserve all diagnostics?
,,''Preserving answers in model runs does not necessarily imply preserved diagnostics. ''

If you answered '''NO''' to any of the above, please provide further details:

  • Which routine(s) are causing the difference?
  • Why the changes are not protected by a logical switch or new section-version
  • What is needed to achieve regression with the previous model release (e.g. a regression branch, hand-edits etc). If this is not possible, explain why not.
  • What do you expect to see occur in the test harness jobs?
  • Which diagnostics have you altered and why have they changed?Please add details here........

System Changes

Does your change alter namelists?'''YES/NO '''
Does your change require a change in compiler options?'''YES/NO '''

If any of these apply, please document the changes required here.......


''Please ''summarize'' any changes in runtime or memory use caused by this change......''

IPR issues

Has the code been wholly (100%) produced by NEMO developers staff working exclusively on NEMO?'''YES/ NO '''

If No:

  • Identify the collaboration agreement details
  • Ensure the code routine header is in accordance with the agreement, (Copyright/Redistribution? etc).Add further details here if required..........

Attachments (2)

Download all attachments as: .zip