Changes between Initial Version and Version 1 of ticket/0927_Energy_diag_v3.5

Timestamp:

2013-04-24T15:52:17+02:00 (11 years ago)

Author:

gm

Comment:

--

ticket/0927_Energy_diag_v3.5

@@ +1 @@
+[[PageOutline]] Last edited [[Timestamp]]
+
+'''Author''' : Gurvan Madec & Fabien Roquet
+
+'''Ticket ''' : #927
+
+'''Branch v3.5''' : [https://forge.ipsl.jussieu.fr/nemo/browser/branches/2013/dev_r3858_CNRS3_Ediag 2013/dev_r3858_CNRS3_Ediag]
+
+----
+'''LOCEAN .3 (2013) - 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. In particular introduce a consistent way of computing PE[[BR]]  (3) add Vallis equation of state in replacement of linear in T and in T-S eos (nn_eos=2)[[BR]]  (4) add the exact Jacket and McDougall   eos (nn_eos=-1) keep for backward compatibility the nn_eos=0 old case[[BR]]  (5) add TEOS10 equation of state[[BR]]  (3) validation + documentation
+
+'''     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'''
+
+'''I. Port from v3.4 to v3.5'''
+
+convert the v3.4 development branch [https://forge.ipsl.jussieu.fr/nemo/browser/branches/2012/dev_r3309_LOCEAN12_Ediag 2012/dev_r3309_LOCEAN12_Ediag] into the following v3.5 development branch [https://forge.ipsl.jussieu.fr/nemo/browser/branches/2013/dev_r3858_CNRS3_Ediag 2013/dev_r3858_CNRS3_Ediag]''' to v3.5'''
+
+see Changeset [https://forge.ipsl.jussieu.fr/nemo/changeset/3874 "[3874]"] ; [https://forge.ipsl.jussieu.fr/nemo/changeset/3876 "[3876]"] and [https://forge.ipsl.jussieu.fr/nemo/changeset/3878 "[3878]"]
+
+'''II. First revision of the eos'''
+
+• introduce Vallis eos as the alternative to UNESCO type eos,
+
+• add systematic calculation of before and now partial derivative of density anomaly with respect to T and s (rab_b and rab_n 4D arrays) and compute rn2 and rrau using rab
+
+To do so, the following changes have been made:
+
+ * Added 4D variable rab_b(:,:,:,jp_ts) and rab_n(:,:,:,jp_ts)  (oce.F90 ; istate.F90, ) and their calculation in eosbn2.F90
+ * Add Vallis eos and suppress the 2 old eos.
+ * introduce the computation of rab_pe, the primitive of rab in eosbn2.F90 module and use it in trdpen.F90
+
+  NB:  rhd and rhop are kept for backward compatibility, but they are meant to become obsolete (especially rhop which could be easily removed, with minor changes in zdfkpp, zdfmxl and zdfric
+
+[[BR]] '''Pending issues''' :
+
+ * diagnostic of trends missing in the following modules :                             '''==>>  To be done !!!!''' '''
+
+  traldf_grif.F90 ;  dynhdf_(iso,lap,bilap and bilapg).F90 ; dynnept.F90 ; dynzdf_(exp,imp).F90 ; dynspg_ts.F90
+
+ * verify the change in sign in transport computation ticket''' '''#1043  (in traldf_iso and griff)
+
+ * 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 trd routine (including mld diag...) create umask_i and vmask_i (2D) fields that mirror task_i field but for the velocity points
+
+----
+'''Fabien addition (copy from the v3.4 page):'''
+
+In oce.F90:
+
+ * 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)
+
+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..........