[[PageOutline]] Last edited [[Timestamp]] '''Author''' : Gurvan Madec & Fabien Roquet '''Ticket ''' : #927 '''Branch v3.4''' : [https://forge.ipsl.jussieu.fr/nemo/browser/branches/2012/dev_r3309_LOCEAN12_Ediag 2012/dev_r3309_LOCEAN12_Ediag] '''The development has been transfer to v3.5''' see [https://forge.ipsl.jussieu.fr/nemo/wiki/ticket/927_Energy_diag_v3.5 v3.5 Ediag wiki page] and Branch v3.5''' : [https://forge.ipsl.jussieu.fr/nemo/browser/branches/2013/dev_r3858_CNRS3_Ediag 2013/dev_r3858_CNRS3_Ediag]''' ---- '''LOCEAN .12 - Energy diagnostics''' ''' Motivation: ''' output 3D trends of tracers, momentum, kinetic energy and potential energy.[[BR]] ''' Status :''' the extraction of trends terms exists, but not the 3D output of the trends [[BR]] ''' Main tasks :''' [[BR]] (1) implement the 3D output of tracers and momentum trends using iom_put [[BR]] (2) compute and output the 3D trends of PE and KE [[BR]] (3) validation + documentation [[BR]] ''' Science Reviewer:''' NOCS guy? [[BR]] ''' System Reviewer:''' NOCS guy? [[BR]] ''' Deadline:''' automn 2013 [[BR]] ''' Priority:''' high [[BR]] ''' Depends on:''' - [[BR]] ''' Principal Investigator : ''' Gurvan Madec and Fabien Roquet (gurvan.madec@locean-ipsl.upmc.fr) [[BR]] [[BR]] ---- ''' Detail of the implementation''' ---- ''' Step I : add the 3D output ''' see [https://forge.ipsl.jussieu.fr/nemo/log/branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC?rev=3316 revision 3316] ''' trdmod_oce''' module and ''' namtrd namelist''' [[BR]] logical flags added in namlist namtrd which now controls what is done with the trends. [[BR]] All the types of treatment of a given trend are available at the same time. The memory requirement will only increase due to the time averaged arrays defined in IOM. {{{ LOGICAL , PUBLIC :: ln_dyn_trd = .FALSE. !: (T) 3D momentum trends or (F) not LOGICAL , PUBLIC :: ln_tra_trd = .FALSE. !: (T) 3D tracer trends or (F) not LOGICAL , PUBLIC :: ln_PE_trd = .FALSE. !: (T) 3D Potential Energy trends or (F) not LOGICAL , PUBLIC :: ln_KE_trd = .FALSE. !: (T) 3D Kinetic Energy trends or (F) not LOGICAL , PUBLIC :: ln_vor_trd = .FALSE. !: (T) 3D barotropic vorticity trends or (F) not LOGICAL , PUBLIC :: ln_glo_trd = .FALSE. !: (T) global domain averaged diag for T, T!^2, KE, and PE LOGICAL , PUBLIC :: ln_dyn_mld = .FALSE. !: (T) 2D tracer trends averaged over the mixed layer LOGICAL , PUBLIC :: ln_tra_mld = .FALSE. !: (T) 2D momentum trends averaged over the mixed layer }}} Add these new logical in the namelist. '''==>>> CAUTION only in the ORCA2_LIM directory''' [[BR]] NB: here is the new name set in revision number 3318[[BR]] ''' trdtra''' module[[BR]] Add a systematic mask of the trend.[[BR]] Change the comments to better describe the purpose of this module. Its purpose is: [[BR]] 'TRA' case: to regroup T & S trends and send them to trd_mod, with, in case of advection, transform the incoming advective fluxes into advctive trend (U.grad[T])[[BR]] 'TRC' case: send trend to ted_mod_trc, with, in case of advection, transform the incoming advective fluxes into advective trend [[BR]] all cases : mask the trend (''' ===>>> PROBABLY add in the module a lbc_lnk so that the trend is defined everywhere''' ) ''' dynadv_cen2 and _ubs''' modules[[BR]] change jpdyn_trd_had into jpdyn_trd_keg. Now in flux form _keg corresponds to the horizontal advection trends and _rvo to the metric terms[[BR]] [[BR]] ''' dynnxt''' module[[BR]] add the output using sum of the total dyn trend (except asselin time filter) ("utrd_tot", "vtrd_tot") and of the asselin time filter trend ("utrd_atf", "vtrd_atf") but with a shift by one time step[[BR]] [[BR]] ''' dynvor, trdvor and trdmod_oce''' modules[[BR]] suppress the call to trd_mod in the jpdyn_trd_dat case (computation of beta.V) add add the calculation of beta.V term in 'trdvor' in jpvor_pvo case. And obviously suppress jpdyn_trd_dat from trdmod_oce [[BR]] Also suppress the jpdyn_trd_had case horizontal advection for the dynamics is 'keg' + 'vor' ; in case of flux form, 'had' is put in 'keg' and the metric terms is put in 'vor' [[BR]] there is now only 10 trends on the dynamics instead of 12[[BR]] [[BR]] ''' trdmod''' module[[BR]] 1- introduce the new logical namelist parameters[[BR]] 2- introduce new subroutines : '''trd_budget''' : computation of the domain averaged T,T^2^, PE, KE trends formerly computes in trd_mod routine)[[BR]] ''' trd_3Diom''': output of the 3D trends using IOM [[BR]] [[BR]] ''' trdicp''' module[[BR]] add in trd_twr routine the computation of the vertical diffusive trend on T & S in case of iso-neutral diffusion (ln_traldf_iso=T). These trends ("ttrd_zdfp", "strd_zdfp") name zdfp for "PURE" vertical diffusion trends are output so that by difference with "zdf" trends we can access to the vertical contribution of the iso-neutral operator [[BR]] [[BR]] ''' iodef.xml''' file[[BR]] add all the trends nick name : '''==>>> CAUTION only in the ORCA2_LIM directory''' [[BR]] [[BR]] ---- ''' Step II : simplification of the structure ''' ''' trdicp and trdicp_oce''' modules[[BR]] see [https://forge.ipsl.jussieu.fr/nemo/log/branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC?rev=3317 revision 3317] 1- suppress trdicp_oce module (put required variables in trdmod_oce parameters)[[BR]] 2- use the jptra_trd_... instead of jpicp_... 3- move trd_budget subroutine from trdmod to trdicp.F90. Suppress trd_icp routine (i.e. trd_2d, trd_3d) as the work is now simply done in trd_budget.[[BR]] 4- rescan all dyn/tra trend indices. add jptra_trd_zdfp for "PURE Kz dissusive trend when ln_traldf_iso=T (see also changes in tranxt where a call to trd_tra is done just before the swap, so that PURE Kz trends can be diagnosed in tra_trd. See also xml file). 5- jptra_trd_nsr and jptra_trd_cdt are 3D trends as their incorporate both surface forcing AND runoff, the later being possibly spread in depth ('''==>>> probably to be changed''')[[BR]] 6- suppress the key_trddyn key trdtra from trdicp[[BR]] ''' trdtra - trddyn - trdmod - trdvor''' modules : impact on ''' TRA''' , ''' DYN''' , ''' TRD''' and almost all ''' TOP_SRC''' modules[[BR]] see [https://forge.ipsl.jussieu.fr/nemo/log/branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC?rev=3318 revision 3318] Change the logic : split DYN and TRA trend diagnostics in separate modules[[BR]] move the call to trd_... from the end of step to the select CASE(jp..._atf) i.e. when the last trend is send to trdtra or trddyn[[BR]] Change the name of the trend indices from jptra_trd_xxx to jptra_xxx and from jpdyn_trd_yyy to jpdyn_yyy. This impact TRA,DYN, TRD but also all TOP_SRC[[BR]] 0- move trdmod_oce into trd_oce[[BR]] 1- trdvor: move the trdvor 1st part from trdmod to trevor module. Suppress the call of trd_vor from step (now called in trddyn for each trend, with the output done for the last trend (_atf). [[BR]] 2- create a trddyn that manage the distribution to iom, global mean, KE, vor, mld [[BR]] 3- trdtra in TRC case pre-process advective trends and then call trdmod_trc ; in TRA case: it regroups T & S trend in one and calla local routine to manage the distribution to iom, global mean, PE, mld[[BR]] 4- trdmod now contains only the trdmod_init routine. Rename all as trdini.F90 and ted_init subroutine 5- suppress all '''key_trd...''' keys (thanks to dynamical allocation) except for the top key: '''key_trdmld_trc'''. ''' ==>>> To be done''' [[BR]] ---- ''' Step III : addition of 3D KE diagnostics ''' ''' trdken''' module[[BR]] see [https://forge.ipsl.jussieu.fr/nemo/log/branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC?rev=3325 revision 3325] 0- bug correction: in trddyn.F90 (ji,jj) ==>> (:,:) in surface stress computation[[BR]] 1- create a trdken.F90 module which compute the KE trends at T-point and out them using IOM TO BE DONE: diagnose the KE local dissipation and associated diffusive fluxes ; diagnose bottom friction in all cases (implicit/explicit) ; in flux form, remove the KE*divh(U) from the H/V advection trends NB: Missing the KE trends in the xml files ; Bug in trdtra: 2 call to ted_tra_mng in case of jptra_zdfp. Both corrected in Step IV release (see below) ---- ''' Step IV : addition of 3D PE diagnostics ''' ''' trdpen''' module[[BR]] see [https://forge.ipsl.jussieu.fr/nemo/log/branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC?rev=3326 revision 3326] and [https://forge.ipsl.jussieu.fr/nemo/log/branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC?rev=3327 revision 3327] 0- correct a bug in trdtra and trdken 1- add xml name for both PE and KE 2- create a trdpen.F90 module which compute the PE trends at T-point and out them using IOM 3- introduce in eosbn2 a subroutine, eos_drau_dtds, that computes the partial derivative off in situ density with respect to T and S (NB: only coded the 2 linear eos, '''TO BE DONE for the UNESCO eos''' (nn_eos=0) 4- correct a few bugs and adopt almost systematic CASE structure in trdtra (see revision 3327) NB: to be verified: the addition of a ssh term in PE for NOT lk_vvl case NB: runoff put as a source term in the divergence computation since v3.3 (see divcur), implication for source term in PE, but also tracer trends ? '''Think about that!''''''''' ---- '''Pending issues''' : atmospheric pressure gradient trend not taken into account (see dynspg.F90 ''' To be done !!!!''' [[BR]] ''' kpp non-local trend put in zdf trends !!! this will not work ! a additional trend term should be add ''' To be done !!!!''' [[BR]] ''' problems to be solved: vvl case for tracer sad trends ; flux form case for had (keg) and zad momentum trends[[BR]] add separate modules for each option ... create the momentum diag over the ML reshape trdtra so that T and S are treated separately in all ted routine (including mld diag...) create umask_i and vmask_i (2D) fields that mirror task_i field but for the velocity points ''''' ''''' ' ---- '''Changes done by Fabien R.''' In oce.F90: * Added 4D variable alpbet: alpbet(:,:,:,jp_tem) is the thermal expansion coef, and alpbet(:,:,:,jp_sal) the haline one. * rhd and rhop are kept for backward compatibility, but they are meant to become obsolete[[BR]](especially rhop which could be easily removed, with minor changes in zdfkpp, zdfmxl and zdfric) * TODO: rn2b might become obsolete (only one rnb variable might be enough). * TODO: check restartibility In step.F90: * Added "IF( .NOT.ln_bfrimp)" before the "CALL dyn_bfr( kstp )" to clarify that dyn_bfr is called only if bottom friction is explicit. * TODO: Diagnostic of BFR should be improved! * added eos calls to set alpha and beta expansion coeffs before bn2 computation. * Consider removing calls to bn2, and compute it only where needed: * rn2b needed to call ldf_slp * also: asmtrj, step_c1d, ldfeiv, zdfddm, zdfevd, zdfgls, zdfkpp, zdfric, zdftke, zdftmx * Also, remove obsolete calls to eos when necessary. * TODO: consider removing calls to eos to get rhop. In trdtra.F90: * removed semi-colon on lines 118 and 267. * replaced ln_glo_trd by ln_PE_trd in line 239. * Moved wrk_alloc and wrk_dealloc in the IF (.NOT.lk_vvl ) structure. In trddyn.F90: * Added wrk_alloc, wrk_dealloc and lbc_lnk in trd_dyn_iom.[[BR]]Remove unused variables ztswu, ztswv. * Following norm in the code, utrd_bfr and vtrd_bfr should be filled only when ln_bfrimp=.FALSE. (non implicit bottom friction). * Implicit bottom friction is calculated at the end of each timestep, in dynzdf_imp.[[BR]]The diagnostic of bottom friction in the implicit case is saved in separate variable names to avoid confusion: utrd_bfri and vtrd_bfri.[[BR]]Similarly, diagnostics of wind stress inputs are provided in 2D fields utrd_tau and vtrd_tau. In tranxt.F90: * Important:[[BR]]If ln_dyn_hpg_imp=.TRUE., Brown and Campana effect is added. The pressure gradient term is computed using a linear combination version of T and S at times b, n and a:[[BR]] Tbc=Tb/4+Tn/2+Ta/4[[BR]] Sbc=Sb/4+Sn/2+Sa/4 In trdglo.F90: * In trd_glo_init, corrected definition for tvolt:[[BR]] tvolt = tvolt + SUM( e1e2t(:,:) * fse3t(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) ) * Treatment of bfr has been clarified: if bfr is implicit (ln_bfrimp=.TRUE.), the bottom friction term remain 0, as bottom friction is implicitly included in the vertical diffusion term. * Correction of outputs of dynamic trends and hke. * What is the test: "pressure gradient u2 = - 1/rau0 u.dz(rhop)" ? Identity is not verified in my configuration! * Wind stress is always zero (although not in the 3d diag provided in trddyn.F90) * In glo_dyn_wri, missing zcof to compute density flux at w-point?! * Conversion is calculated as -g*div(rho*U)/rho0 !! (missing z) * TODO: Once trddyn and trdtra are ready, use them for trdglo computations... In trd_oce.F90 {{{ INTEGER, PUBLIC, PARAMETER :: jptot_dyn = 13 !: Total trend nb: change it when adding/removing one indice below ! =============== ! INTEGER, PUBLIC, PARAMETER :: jpdyn_hpg = 1 !: hydrostatic pressure gradient INTEGER, PUBLIC, PARAMETER :: jpdyn_spg = 2 !: surface pressure gradient INTEGER, PUBLIC, PARAMETER :: jpdyn_keg = 3 !: kinetic energy gradient or horizontal advection INTEGER, PUBLIC, PARAMETER :: jpdyn_rvo = 4 !: relative vorticity or metric term INTEGER, PUBLIC, PARAMETER :: jpdyn_pvo = 5 !: planetary vorticity INTEGER, PUBLIC, PARAMETER :: jpdyn_zad = 6 !: vertical advection INTEGER, PUBLIC, PARAMETER :: jpdyn_ldf = 7 !: horizontal diffusion INTEGER, PUBLIC, PARAMETER :: jpdyn_zdf = 8 !: vertical diffusion INTEGER, PUBLIC, PARAMETER :: jpdyn_bfr = 9 !: bottom stress INTEGER, PUBLIC, PARAMETER :: jpdyn_atf = 10 !: Asselin time filter INTEGER, PUBLIC, PARAMETER :: jpdyn_tau = 11 !: surface stress INTEGER, PUBLIC, PARAMETER :: jpdyn_bfri = 12 !: implicit bottom friction (ln_bfrimp=.TRUE.) INTEGER, PUBLIC, PARAMETER :: jpdyn_ken = 13 !: use for calculation of KE }}} * where jpdyn_tau and jpdyn_bfri are diagnostics of surface and bottom stress, respectively. * jpdyn_ken has been added for the calculation of KE in dynnxt.F90 (done just before the swap) In trd_ken.F90: * l110: DO ji = 2, jpj should be DO ji = 2, jpi * Same in l133, l152 and l178. * In l183, replaced "ketrd_bfr" by "ketrd_bfri" to avoid confusion with "ketrd_bfr". * Important: contribution from atf is calculated with un and vn, after the swap.[[BR]]A small error is made: the contribution of ATF term at time t is recorded with trends taken at t-1. * In trd_ken_init, the 3 variables fse3x_n were replaced by their constant value fse3x for the non-vvl case. * Definition of KE (KE in xml) at time t+1/2: KE(t+1/2)=rau0*u(t)*u(t+1)/2.[[BR]]Add also subroutine ken_p2k to compute conversion rate.[[BR]]KE and ketrd_convP2K both calculated during dynnxt.F90, using ktrd=jpdyn_ken * l. 109: multiplication by rau0 of zke to obtain KE trends in W/m3, and KE in J/m3. * remove r1_2_rau0 variable, which is no more used. In dynspg_flt: * Add a diagnostic of the explicit and implicit (due to filter) contributions to SPG. * ssh_flt: diagnostic of the ssh modification due to filter. In dynnxt.F90 * computation of z1_2dt must be put before the IF( ln_dyn_trd ) block (l. 194 and 195) In eosbn2.F90 * Vallis equation of state added (Vallis 2006, p34-35): simple EOS which accounts for thermobaricity, cabelling and compressibility.[[BR]]Generalization of linear equation of state, so T-linear and TS-linear cases have been removed. [[BR]]Vallis EOS: nn_eos=1. * Add original [wiki:McDougall Jackett and McDougall] (1995) EOS (nn_eos=-1)[[BR]]nn_eos=0 case is a modified version of the Jackett and McDougall (1995) EOS!! Numerically close though... * eos_alpbet modified, which now provides alpha and beta (as a 4d variable) instead of alpha/beta and beta0=0 or 1.[[BR]]ldf_slp_grif in ldfslp.F90 modified accordingly (only place eos_alpbet was used).[[BR]]Now, there is no problem if beta=0, because its inverse is never used. * bn2 modified: * use an exact formulation based on each eos * bn2(ts,pn2): call first eos_alpbet to compute alpha and beta coefficients, then interpolating alpha/beta on w-points and finally calculating pn2. * bn2(ts,alpbet,pn2): alpha and beta coefficients are directly given in arguments, then interpolated to compute pn2. * Consider removing bn2 from eosbn2.F90 * TODO: Modify trabbl.F90 to use alpbet * TODO: Modify tranpc.F90 to use alpbet * TODO: Modify zdfkpp.F90 to use alpbet * eos has now a unified interface: * eos(ts,rhd) gives density anomaly (3D field) * eos(ts,rhd,rhop) should become obsolete (density anomaly + sigma0) * eos(ts,dep,rhop) density anomaly at one zt(pdep) (2D) * eos(ts,alpbet) providing alpha and beta (4D field jpi,jpj,jpk,jpts) * eos(ts,rhd,alpbet) density+alpha/beta (3D + 4D fields) * eos_pen added to compute PE anomaly, and equivalents of alpha and beta for the PE state variable (called alpha_pe and beta_pe, given in alpbet_pe).[[BR]]function is called at each time step trdpen.F90 (PE diagnostic).[[BR]]PE anomaly is defined as: (PE - rau0*gz)/(rau0*gz). For a linear case, PE anomaly is equal to density anomaly (nice isn't it?)[[BR]]When z=0, PE anomaly is defined asymptotically, converging toward the density anomaly value (in the code, if z<0 --> z=1).[[BR]]TODO: Waiting for a new EOS from Trevor McDougall that would allow a more efficient computation of each terms. In trdtra.F90 * Add 'trdt(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)' in CASE jptra_bbc in trd_tra (l. 103) In trdpen.F90 * replace petrd_ldf by petrd_zdf in l. 102 * renamed petrd_for as petrd_nsr for consistency with trd_oce (l. 107) * remove lines on petrd_sad in CASE jptra_atf in trd_pen: petrd_sad cannot be called 2 times. Furthermore, contribution of asselin filter on ssh must be diagnosed in ssh_nxt, not in tra_nxt!! * Remove '( nn_eos == 0 .OR. nn_eos == 3 ) .AND. ' in the condition to call eos_pen in trd_pen (l. 79). These two coefficients are updates at each timestep whatever the eos (very quick anyway for the linear case!) * Add 'CALL iom_put( "PEanom", zpe )' in trd_pen, just after the call to pen_ddt_dds. * Replace ptrdx by ptrdy and remove multiplication by fsde3w to compute pe trends in trd_pen (l. 88-89) Bug during compilation: * make: Circular dynhpg.o <- trddyn.o dependency dropped[[BR]]patched by compiling a first time with the line 'USE trddyn' commented in dynhpg.F90, then recompiling with 'USE trddyn' decommented. Documentation: * In Annex A, * change label name Apdx_A_grad_p in Apdx_A_grad_p2 l. 418 and 430. * sign error for vertical advection in Eq. A.18 (tracer equation) * In Annex C, sign errors corrected in discretization of HPG (p. 293) * In Chap_TRA, * update part 5.8 on equation of state to include reference to Vallis2006+JM95 true formulation * In Chap_DYN, * correction of title 6.1.2, and references * Changes on Asselin description * In Chap_STP, modification on Asselin description General comments: * implement e3t = dk(gdept), and e3w=dk(gdpew). * N2 calculation could be speed up by precalculating ratios: (zw[k]-zt[k])/e3w[k] and (zt[k-1]-zw[k])/e3z[k] * put KE and PE directly in diawri? * Energetics of Brown and Campana effect? * Modify HPG to use alpha/beta + partial step formulation on fixed z-levels ---- 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[[BR]][[BR]]New code doesn't change results when it is switched off[[BR]][[BR]]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?[[BR]][[BR]]!,,!''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 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....... ---- 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..........