334 rho sigma diagnostic regridding#348
Conversation
|
The main work of the PR is to introduce a new module diag_remap which contains all code to do with vertical remapping of diagnostics. This make diag_mediator a lot cleaner. |
nichannah
force-pushed
the
334-rho-sigma-diagnostic-regridding
branch
from
4e3b25c to
5240529
Compare
- The umo,vmo CMOR diagnostics were using in-place global multiplications which is known to lead to unitialized variable errors within the diag manager. - Registering with conversion=GV%H_to_kg_m2 instead so that the diag_mediator can create a local array with which to do the conversion. - No answer changes.
- register_diag_field now only has a two-way split and (used to be four
for native, CMOR, z, z-CMOR). Call tree is now:
- register_diag_field() calls register_diag_field_expand_cmor() for both
the native and z-remapped field.
- register_diag_field_expand_cmor() calls register_diag_field_expand_axes() for both
the native names and CMOR names. This is where the diag_type is allocated and filled.
- register_diag_field_expand_axes() calls register_diag_field_fms() one of four
ways depending on presence of area interp_method optional arguments.
- Fix: The z-diagnostic cell_methods was incorrectly using the native vertical
dimension name for some fields when the diagnostic was not requested. This
led to the wrong names being provided in available_diags. The re-factor corrected
this unexpectedly.
- No answer changes.
- Removed functions is_u_axes(), is_v_axes(), is_B_axes(), is_layer() by storing logicals in the axes_group type. - Also added needs_remapping logical for non-native 3D axes groups to indicate to do remapping. - No answer changes.
- Interface variables cannot be remapped so our z-coordinate diagnostics have been limited to layer-centered fields. This implements interpolation of interface-centered fields to the interfaces of the target grid (same as used by the diagnostic remapping). - This added a new group of axes for the z-coordinate interfaces. - Note that the (horizontal) B-grid located diagnostics are not implemented (as they are not implemented in remapping code either). - Had to provide a work-around for the current existence of MOM_diag_to_Z which calls register_diag_field for 3D fields that should not be remapped. Todo: defaults for is_native should be flipped. - No answer changes.
- Some layer data is already vertically integrated (e.g. hT, uh, vh, ...) and cannot be remapped as if it were a concentration (intrinsic variable). Setting optional argument v_extrinsic=.true. to register_diag_field() now indicates that the field should be re-gridded as an extrinsic field. - No answer changes.
- 3D diagnostics h, uh, vh have been labeled as vertically extrinsic (meaning integrated over the layer) so that remapping does not treat the fields as concentrations. - Changes availabl_diags. - No answer changes.
| character(len=5), parameter :: REGRIDDING_SIGMA_STRING = "SIGMA" !< Sigma string | ||
| character(len=6), parameter :: REGRIDDING_LAYER_STRING = "LAYER" !< Layer string | ||
| character(len=6), parameter :: REGRIDDING_ZSTAR_STRING = "ZSTAR" !< z* string | ||
| character(len=6), parameter :: REGRIDDING_RHO_STRING = "RHO" !< Rho string |
There was a problem hiding this comment.
I'm curious about why are you changed the declared string lengths to not match the string? I think it's a harmless change but was there a failure somewhere?
There was a problem hiding this comment.
I think I did that because I put these into an array below. My understanding is that Fortran doesn't like having variable length strings in an array. I'll add a comment.
There was a problem hiding this comment.
I see that now. I'm surprised that character(len=*), dimension(n), parameter = works - I would have thought the * is at odds with the parameter...
|
I'm taking a crack at merging our branches. I think it's going to work out. Just giving you the heads up in case we are both working on the merge? |
- New subroutines and functions should be doxygenized at inception.
- The vector diag_remap was a vector of control structures. To help me understand the new code I'm renaming by adding _cs. - Note this code does not compile due to a broken merge of nicjhan-334-rho-sigma-diagnostic-regridding on dev/master.
- For each diag_remap_ctrl element there are now also eight axes_groups just as there are eight native axes_groups. - Not being used yet! - Still doesn't compile after a broken merge of nicjhan-334-rho-sigma-diagnostic-regridding on dev/master.
- The coexistence of axes_groups and axes ids is confusing. I've renamed %axes to $axes_ids to help me understand the code. - Note this code does not compile due to a broken merge of nicjhan-334-rho-sigma-diagnostic-regridding on dev/master.
- This keeps the new name "ZSTAR" but also allows use of the old coordinate name "Z*". - Note this code does not compile due to a broken merge of nicjhan-334-rho-sigma-diagnostic-regridding on dev/master.
…o-sigma-diagnostic-regridding - Moved new inteprolation and reintegration routines from MOM_diag_mediator (added on branch user/aja/diag_interp_and_extrinsic) to MOM_diag_remap. - Added vectors of axes_groups to correspond to each diag_remap_cs. - Added workarounds for keeping _z suffix in module name and input parameter names with ZSTAR. - Compiles and passes regression tests. - A few quick plots indicate remapping, interpolation and reintegration are working as well as they were on branch user/aja/diag_interp_and_extrinsic. # Conflicts: # src/framework/MOM_diag_mediator.F90
- When the source or destination columns were completely vanished the reintegrate_column() in MOM_diag_vkernels was getting stuck in a never ending loop. - Added exit tests for these special cases. - Also added unit tests to cover these scenarios.
- To avoid a conflict with the general coordinate diagnostics, I've renamed the diagnostics module registered within MOM_diag_to_z to ocean_model_zold.
User/aja/diag interp and extrinsic
|
@nicjhan Are you satisfied with this now? Can I pull into dev/master? |
|
@adcroft yes, that would be great. I have an offline testing issue here NOAA-GFDL/MOM6-examples#98 . The only test we currently have is "can all expected diags for each vertical coord be output". If you have ideas for others that you think could be done externally (i.e. by Python code after a run) can you please put them here. I think @angus-g has some ideas too. |
- In MOM_ALE.F90, we read INTERPOLATION_SCHEME which is only used for interpolating the model state to place coordinates when the coordinate is state dependent. We now do not read or log INTERPOLATION_SCHEME unless the coordinate mode use state interpolation.
- The function that returns whether a coordinate is state-dependent now takes either the string name or the integer scheme as dummy argument.
- New function to return the n'th word in a string as an integer - Updated unit tests in MOM_string_functions
- Partial re-factor in preparation for combining coordinate definitions between ALE and diagnostics. - Added WOA09 vertical grid to be used by diagnostics.
- initialize_regridding() now takes param_file as argument and reads most parameters to configure the coordinate. Parameter names are generated from prefix/suffix arguments or default to main ALE parameters (a bit messy but meh). - Moved configure_coordinate() from MOM_ALE into initialize_regridding(). - MOM_ALE no longer does the coordinate configuration inline but calls initialize_regridding(). - Removed allocate_regridding() from MOM_regridding. - Moved dz_function1() from MOM_ALE to MOM_regridding. - MOM_diag_remap now calls initialize_regridding to configure the diagnostic coordinates. - This modifies the parameter documentation.
- When the maximum depth is >3000m it is reasonable to assume the model is global so the default diagnostic z grid is the WOA09 vertical grid.
- We were using the coordinate name to construct the dimension names but with multiple diagnostic coordinates we need to use the label given to the diagnostic coordinate which has to be unique for it to work with the diag_manager.
- UNIFORM can now set number of layers and the resolution, needed for diagnostic coordinates that should not use model nk nor he max_depth parameter. - Added extract_real() to MOM_string_functions for easy extraction of reals from parameter expressions
- An optional argument was no longer needed
- In preparation for adding xy-average diagnostic for all diagnostic coordinates.
- Some interior points were being masked. We only mask z*-based coordinates and so only vanished target cells need be giving missing values.
- Added horizontally_average_diag_field() to MOM_diag_remap.F90 - Added 1D axes group for remapped axes - Removed empty layer egister_diag_field_3d_and_1d() - Also removed unused argument from diag_remap_configure_axes()
|
@nicjhan , here's an update on where this stands: Good news it that the code works. It took longer than expected to get a production run going because we missed a whole class of diagnostic that are used for monitoring the production runs (by dora); in the old z-diagnostics we provided horizontally averaged diagnostic for each field. I've added the same for each new diagnostic-coordinate on my merge of your branch nicjhan-334-rho-sigma-diagnostic-regridding. Bad news is that the production run is not reproducing it's control (based on atmospheric state) which makes evaluating the new diagnostics somewhat difficult (I wanted to run our analysis using the new diagnostics on the same solution). The ocean code is bit-wise passing all our tests but something is very wrong somewhere and now I need to find out what... |
|
@jkrasting Looks like the zonal average analysis isn't working with the new diagnostics because of an assumed name "zw" in plotVars3.py . |
|
@adcroft thank you for the update. It looks like you've made a lot of changes over the past week. I've merged in your branch and am running the ocean reproducibility tests, double checking for any inconsistency there. Is there a way for me to approximate your control run? I'd like to reproduce the problems you're seeing. Note: I needed to make the following change to avoid an assert in horizontally_average_diag_field(): nichannah@58b26a0 |
- The old z-diagnostics had "edge positions" in the output file in addition to the level centers. This was missed in the generalized diagnostic coordinate code.
- Extensive and intensive flags in types and arguments to functions were written as extrinsic and intrinsic.
- @nicjhan correctly disabled the xy-averaging of interface data in 58b26a0 because it wasn't implemented. - Now it is. - Interface data does not use thickness weights. - Layer data that is extensive also does not use thickness weights.
- Adaption of https://github.com/NOAA-GFDL/MOM6/wiki/MOM6-diagnostic-vertical-remapping written by @nicjhan
- The sub-directory mode of doxygen meant we could not predict the url for specific pages (e.g. Diagnostics.html) which makes linking to those pages from outside hard.
|
This branch is ready to merge into dev/master but the associated tools (diag_table, analysis and dora) have the following issue which I needs to fix before we can use it in production. First, all of the automatic analysis that does not use transport seems to work. Woot! Second, for the _z diagnostics we had in the diag_table (u,v,s,t,uh,vh,age) we added ~200s to a run-time of order ~8000s. I think that's acceptable although I've yet to try doing all the CMIP6 diagnostics in _z. Third... The old z-diagnostics had the units for vh reported incorrectly (files said m3/s but was actually kg/s) and we had unwittingly compensated for this in the analysis (converting vh to Sv with a 1e-9 multiplier). The old _z code looks like it is meant to write uh and vh in m3/s for Boussinesq mode and in kg/s for non-Boussinesq mode. Unfortunately it actually always wrote in kg/s even though the reported units were as intended (i.e. m3/s in Boussinesq mode). The analysis scripts thus were able to always use a 1e-9 conversion for vh, independent of reported units. To date the analysis tools have been using uh and vh which are model variables so their units should depend on the Boussinesq mode. The CMIP6 reporting requirements for transport use kg/s regardless of mode (and have redefined 1 Sv = 1e9 kg/s - I wish I'd thought about this before publication but it's a minor oversight). The way out of this is to use umo,vmo for analysis which are always written in kg/s as requested for CMIP6. We already have umo,vmo diagnosed (in the right units) in the code. but we've only had umo,vmo in the diag_able in the native component (not _z component). Updating the diag_table and analysis scripts to use umo,vmo is trivial but the tools on dora are going to need some special workarounds to work on both old and new experiments. Since dora is our primary tool for analysis a.t.t. I need to work with @jkrasting to fix this before making the final merge. All the above explains the "5%" that @nicjhan reported (Nov 3 in this thread). The following figures use vh with either the 1e-9 or 1e-6 conversion. I believe the "new_diags" with 1e-6 conversion are the most consistent (since the 1e-9 did not account for the Boussinesq density not equal to 1000). In other words, MOC has been overestimated by ~2-4% until now (I won't be losing sleep over that...).
|
- umo,vmo, are supposed to include SGS transports but were in fact only the resolved transports. - These diagnostics and their vertical sums have been moved to where uhtr,vhtr are diagnosed (they are the same apart from unit conversions).
nichannah
force-pushed
the
334-rho-sigma-diagnostic-regridding
branch
from
c3d37fa to
c838594
Compare
4 participants
FYI and review only. Presently as conflicts with dev/master and should be merged after:
NOAA-GFDL#347
This branch introduces generalised vertical remapping for diagnostics.