Version 35 (modified by mlelod, 10 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-½ and n+½

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+½ and kt+½ 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+½ and kt-½ 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

sbcmod.F90 : add the swap of surface fields and the read/write in the restart file.

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+½ and kt-½). 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.

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 : applied the solar heat flux (qsr) as the mean of 2 consecutive time-steps (kt+½ and kt-½)

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 !

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

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

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

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

sbc_oce : Declaration of new forcing fields needed for modified forcing and modified filtering

Heat content trend due to SBC non solar fluxes (2D) at now and before time steps stored in sbc_trd_hc_n and sbc_trd_hc_b.

Heat content trend due to SBC fluxes (2D) at now and before time steps stored in sbc_trd_hc_n and sbc_trd_hc_b.

Heat content trends due to penetrative solar radiation (3D) are stored at now and before time steps in qsr_trd_hc_n and qsr_trd_hc_b.

sbcmod : Swap and restart of these before new forcing fields

Add qsr_trd_hc_b, sbc_trd_hc_b and sbc_trd_sc_b in swap and restart read/write.

trasbc : Changes for modified forcing integration

sbc_trd_(h/s)c_n is aded to ta in forward Euler time stepping case

0.5 * ( sbc_trd_(h/s)c_b + sbc_trd_(h/s)c_n ) is aded to "ta" in the leap-frog time stepping case

traqsr : Changes for modified forcing integration

The 3D heat content tendency term of penetrative solar flux is stored in "qsr_trd_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 : 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.

Then only qsr_trd_hc_n variable is aded to ta in the forward Euler time stepping case and 0.5 * ( qsr_trd_hc_b + qsr_trd_hc_n ) is added to ta in the leap-frog time stepping case.

sshwzv : Sea surface elevation time stepping

comment for the forward Euler time stepping case: "In the forward Euler time stepping case, the same formulation as in the leap-frog case can be used because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp"

surface weighted averaging for the computation of sshf_a and sshf_b in vvl case (as for sshu and sshv)

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 * ( 'volume_forcing_n' - 'volume_forcing_b' ) from volume time filter

sshwzv.F90 : a additional term appears in the filter (the last part of the RHS) :

sshn = ssha + atfp * ( sshb -2 sshn + ssha ) - atfp * rdt * ( emp_b - emp ) / rau0

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+½, 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+½ 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