Version 16 (modified by mlelod, 14 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://gm@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/DEV_r1837_MLF/NEMO * commit : svn ci m "ticket: #663 blah_blah" list_of_file
Description
Introduce the modified LeapFrog  RobertAsselin filter (LFRA) 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 nonconservation of the standard LFRA 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 LFRA scheme, so that it is actually a quasi second order scheme.
The LFRA 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 timestepped first, followed immediately by the computation of the now vertical velocity (see figure below).
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 n1/2 and n+1/2
modules involved : sbc_oce.F90, sbcmod.F90, trasbc.F90, traqsr.F90, dynzdf....F90
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 seaice!) must now provide both kt+1/2 and kt1/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
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 timesteps (kt+1/2 and kt1/2). This involes the following changes : (1) incorporate the time mean of the forcing ; (2) put the IFELSEENDIF outside the DOloop 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 timesteps (kt+1/2 and kt1/2)
CAUTION 1 : the case where the light penetration is provided by the biomodel ( 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 timestepping 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 timesteps (kt+1/2 and kt1/2)
dynzdf_imp.F90 : applied the stress (utau, vtau) as the mean of 2 consecutive timesteps (kt+1/2 and kt1/2)
sshwzv.F90 : applied the freshwater flux (emp) as the mean of 2 consecutive timesteps (kt+1/2 and kt1/2)
(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 leapfrog 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 biomodel ( 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 timestepping 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 leapfrog 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 leapfrog 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. Therefore dynnxt.F90 is not changed.
(3) The semiimplicit hpg.
In tranxt.F90, for the semiimplicit 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
tn+rbcp*(ta2*tn+tb)
with rbcp=(1+atfp)*(1+atfp)*(1atfp)/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.
*_m 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
N.B. : Three key elements of the LFRA 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, nomore 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 nonlinear 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:fromNEMOteamshortterm)
 Add 'real 3D' treatment of fse3(t/u/v)_b:
fse3t_b is not the repartition of sshb. Local conservation implies that the filter correction done on ssh is concentrated in the first level and not distributed on the vertical. => full 3D field is needed. fse3(u/v)_b are deduced from fse3t_b by a surface averaging (as ssh(u/v) from ssh) Add a non corrected before ssh field, sshbnc, in the restart file in order to compute fse3ub. This variable won't be used in the code, just for the restart.
Testing
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  '''YES/NO/NA''' 
(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. ''  '''YES/NO''' 
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 sectionversion
 What is needed to achieve regression with the previous model release (e.g. a regression branch, handedits 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.......
Resources
''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)

NEW_LF_RA.jpg
(53.6 KB) 
added by gm 14 years ago.
step flowchart for the modified LFRA scheme

NEW_LF_RA_FlowCart.jpg
(79.1 KB) 
added by gm 14 years ago.
step flowchart for the modified LFRA scheme
Download all attachments as: .zip