MAM4xx: Fixing the ne120 case for the microphysics interface. #3095
Conversation
|
Tagging @susburrows |
There was a problem hiding this comment.
I think this should be reworked.
- Why don't you consider fixing the files so that they have the correct dimensions needed?
- Did you check if the units align? I don't think LEV and altitude would have the same units ... but I could be wrong
- If you want to impl this in code, I think you'd want to potentially pass some sort of a (old_name, new_name) map/dict for these renaming ops instead of overloading the functions
- If you want to impl this in code, in the added logic, I don't think you should differentiate between remapping or not (you are not seeing its effect because the no-remapping branch is no-op, but fundamentally the setup is wrong and so the code should be applied to all)
- it seems this had another bug, which makes me think, you should be adding a simple test that does the remapping (say from ne2 to ne4 or something)
Because we are relying on the vertical interpolation from EAM that takes altitude instead of levels. See here, which is invoked here.
Because this part relies on horizontal interpolation, neither levels nor altitude should be involved in the computation. Note, the third dimension (for columns, I will use two dimensions even though in the code we only have one) is involved, but I do not think units will affect this step. We followed the SPA code, where horizontal interpolation is done first. Then, vertical interpolation is computed on the resulting finer grid.
I use overloading because we only use this part of the code for vertical emissions. Linoz and oxidant reads do not need to be involved. I am open to changing this code. Could you please provide more details on what you are proposing?
I've created two code blocks for remapping and no remapping to avoid cloning the source grid in the no remapping case, where renaming the LEV tag is necessary. This remapper object requires a target grid and a map file as input. In cases where no map files are needed, the target and source grid are the same object. We rename the LEV tag for the target grid here. For this reason, you do not see the renaming of the LEV tag in the source grid.
Yes, we can add another test. @singhbalwinder , could you please tell me where I can find a map file for the horizontal interpolation from ne2 to ne4? |
I meant: if I understand your code, all your code is doing is renaming the dims. Why not rename the dims in the input files? Then, most of this PR won't be needed. |
Because altitude and lev represent different quantities, as you mentioned, their units differ: altitude is measured in km, while lev is in 'hPa'. Thus, using lev with units of km would appear strange. |
There was a problem hiding this comment.
I will defer to Luca
I think this PR requires two things before completion:
- you should add a test for this, given that the initial impl was wrong and buggy; a test will catch stuff
- i would simplify the code to make it easier to maintain in the long run, i don't like the unnecessary overloading (but that's likely a personal opinion) and even if you chose to keep overloading, I think you have to do it more logically than by branching arbitrarily due to the presence of a map or not
| const auto io_grid = horiz_remapper->get_src_grid(); | ||
|
|
||
| // NOTE: If we are using a vertical emission NC file with altitude instead of levels, | ||
| // we must rename this tag. This is only necessary when a map file is used. | ||
| if(tracer_data.file_type == VERT_EMISSION && extfrc_map_file != ""){ | ||
| auto horiz_interp_src_grid = | ||
| io_grid->clone("tracer_horiz_interp_src_grid", true); | ||
| horiz_interp_src_grid->reset_field_tag_name(LEV, "altitude"); | ||
| horiz_interp_src_grid->reset_field_tag_name(ILEV, "altitude_int"); | ||
| return std::make_shared<AtmosphereInput>(tracer_data_file, horiz_interp_src_grid, io_fields, | ||
| true); | ||
| } else { | ||
| return std::make_shared<AtmosphereInput>(tracer_data_file, io_grid, io_fields, | ||
| true); | ||
| } | ||
|
|
There was a problem hiding this comment.
i still don't think you should overload, see my comments previously --- just add this stuff to the other function with some extra args
if you still want to, consider this change
| const auto io_grid = horiz_remapper->get_src_grid(); | |
| // NOTE: If we are using a vertical emission NC file with altitude instead of levels, | |
| // we must rename this tag. This is only necessary when a map file is used. | |
| if(tracer_data.file_type == VERT_EMISSION && extfrc_map_file != ""){ | |
| auto horiz_interp_src_grid = | |
| io_grid->clone("tracer_horiz_interp_src_grid", true); | |
| horiz_interp_src_grid->reset_field_tag_name(LEV, "altitude"); | |
| horiz_interp_src_grid->reset_field_tag_name(ILEV, "altitude_int"); | |
| return std::make_shared<AtmosphereInput>(tracer_data_file, horiz_interp_src_grid, io_fields, | |
| true); | |
| } else { | |
| return std::make_shared<AtmosphereInput>(tracer_data_file, io_grid, io_fields, | |
| true); | |
| } | |
| auto io_grid = horiz_remapper->get_src_grid(); | |
| // NOTE: we must rename this tag to be consistent | |
| io_grid->reset_field_tag_name(LEV, "altitude"); | |
| io_grid->reset_field_tag_name(ILEV, "altitude_int"); | |
| return std::make_shared<AtmosphereInput>(tracer_data_file, io_grid, io_fields, | |
| true); | |
Even better, if you could just add this stuff in the first call as a map
inline std::shared_ptr<AtmosphereInput> create_tracer_data_reader(
const std::shared_ptr<AbstractRemapper> &horiz_remapper,
const std::string &tracer_data_file
std::map<FieldTag, std::string>* some_map = nullptr) {
std::vector<Field> io_fields;
for(int i = 0; i < horiz_remapper->get_num_fields(); ++i) {
io_fields.push_back(horiz_remapper->get_src_field(i));
}
auto io_grid = horiz_remapper->get_src_grid();
if (some_map) {
for (const auto& [key, value] : *some_map) {
io_grid->reset_field_tag_name(key, value);
}
}
return std::make_shared<AtmosphereInput>(tracer_data_file, io_grid, io_fields,
true);
} // create_tracer_data_readerI am assuming this stuff isn't happening on device, and you have to double check the code
There was a problem hiding this comment.
He can't change the remapper src grid, since it's const. He'd have to (shallow) clone it first.
There was a problem hiding this comment.
Good point, forgot that ... in the above:
- auto io_grid = horiz_remapper->get_src_grid();
+ auto io_grid = horiz_remapper->get_src_grid().clone(...);There was a problem hiding this comment.
I committed a cleaner version of this part of the code. I will clone only the source grid if the file type corresponds to a file with altitude, i.e., the VERT_EMISSION file type. I did not use a map because it would require additional changes in the interface code (I need to create a map). I deleted the overload of create_tracer_data_reader by passing the file type as input and setting the file type default to NONE.
| @@ -306,7 +308,7 @@ void MAMMicrophysics::set_grids( | |||
| // I am assuming the order of species in extfrc_lst_. | |||
| // Indexing in mam4xx is fortran. | |||
| forcings_[i].frc_ndx = i + 1; | |||
| const auto io_grid_emis = VertEmissionsHorizInterp_[i]->get_src_grid(); | |||
| const auto io_grid_emis = VertEmissionsHorizInterp_[i]->get_tgt_grid(); | |||
There was a problem hiding this comment.
To create a test, I will need to have a map file. Is it possible to include this test in a future PR? Note that the evaluation team is running additional cases that allow us to catch issues.
There was a problem hiding this comment.
I think this is fine. In the devs/IO group we are working on unifying the way we read data from file, so I would not spend too much time in polishing how your read data. So long as it works, and the details/gotchas are clear (e.g., inline comments), then I'm fine.
| const auto io_grid = horiz_remapper->get_src_grid(); | ||
|
|
||
| // NOTE: If we are using a vertical emission NC file with altitude instead of levels, | ||
| // we must rename this tag. This is only necessary when a map file is used. | ||
| if(tracer_data.file_type == VERT_EMISSION && extfrc_map_file != ""){ | ||
| auto horiz_interp_src_grid = | ||
| io_grid->clone("tracer_horiz_interp_src_grid", true); | ||
| horiz_interp_src_grid->reset_field_tag_name(LEV, "altitude"); | ||
| horiz_interp_src_grid->reset_field_tag_name(ILEV, "altitude_int"); | ||
| return std::make_shared<AtmosphereInput>(tracer_data_file, horiz_interp_src_grid, io_fields, | ||
| true); | ||
| } else { | ||
| return std::make_shared<AtmosphereInput>(tracer_data_file, io_grid, io_fields, | ||
| true); | ||
| } | ||
|
|
There was a problem hiding this comment.
He can't change the remapper src grid, since it's const. He'd have to (shallow) clone it first.
…o change tag Lev of vertical emissions.
There was a problem hiding this comment.
I am changing my review to "Comment" so that you can feel free to override and merge. However, I am not a fan of this PR, because I think it's covering sloppy code by more sloppy code. That's not a sustainable strategy. You asked me to review it, and I am giving you my reasons why I don't like it. I think my reasons are well-intentioned and logical; and I am not convinced my concerns have been addressed.
Maybe there's miscommunication, so here's another way to request simple clarifications related to my review:
- Could you please document in this PR the differences between the different file types in the TracerFileType enum? This should be done with a call like
ncdump -h <file>where<file>will be one example of each type (FORMULA_PS, ZONAL, VERT_EMISSION, NONE). It's not clear to me why you're branching the logic, but maybe I am missing something in the file types. The ncdump details will clarify that. - I still don't understand why you're resisting adding a test. I think you should. You have three options (from least desirable to more desirable):
- Document clearly in this PR a way to reproduce the failure before the PR, and how the fixes in this PR actually address the failure. Ideally, this should be something like a script from your testing.
- Add a testmod where the ne120 case (in the title of the PR) can be run with a simple cime call. The full test name could then be something like SMS_D_Ln5.ne120pg2_ne120pg2.F2010-SCREAMv1.$machine_$compiler.scream-mam4xx-micro-remap. This test won't be added to tests we run regularly.
- Add a testmod where a simpler ne2pg2/ne4pg2 can be run. This could be something like: SMS_D_Ln5.ne4pg2_oQU480.F2010-SCREAMv1-MPASSI.chrysalis_intel.scream-mam4xx-micro-remap.
An easy approach is to simply find a remap file that can remap ne30pg2 files to ne4pg2 and then you have everything ready.I don't think we allow coarsening for these, so the test will require some more work, but it's not difficult.
|
Thanks, Naser, for re-explaining everything in detail. We will address your concerns before merging it. Here is how we are going to address all your concerns, please let us know if you need anything else.
|
mahf708
dismissed
their stale review
so that authors can override request for edits if desired
Merge ProtectionsYour pull request matches the following merge protections and will not be merged until they are valid. 🟠 Enforce checks passingWaiting checks:
|
There was a problem hiding this comment.
Thanks for providing the header info from the files.
I wrote a lengthy comment about the confusion. You're encouraged to edit the misleading comments and names.
I also think you should add a test of your choosing so that when we rewrite all this IO stuff we won't break your setup.
I won't have time to re-review this PR again, so feel free to do what you want with my feedback. I am approving and there's no need for my review/approval again.
| true); | ||
| } else{ | ||
| // We do not need to rename tags in or clone io_grid for other types of files. |
There was a problem hiding this comment.
This is partly the confusion and that's why I was unhappy about this PR. Your comments here don't add up to what's going on. What you should actually say is the following: "we must not rename the tags". If you rename the tags here for the zonal or formula_ps files, you will end up with errors. (You will be to rename the tags for the surface/none type though.) So, your branching logic is justified, but that wasn't clear to me (because the notes here made it sound it was just a preference; if it were just a preference, it'd simply be bad code imo --- we should only branch if we have to)
Based on the info you provided in the main post, I believe the following is the best way to describe this messy situation. Let's take your enum for illustration:
enum TracerFileType {
- // file with PS ncol, lev, and time
+ // file with ncol, lev, ilev, time and has P0 and PS as variables
+ // example: oxidants
FORMULA_PS,
- // nc zonal file from ncremap
+ // file with ncol, lev, ilev, time
+ // example: linoz
ZONAL,
- // vertical emission files
+ // file with ncol, altitude, altitude_int, time
+ // example: elevated (at a height) emissions of aerosols and precursors
+ // NOTE: we must rename the default vertical tags when horiz remapping
+ // NOTE: we vert remap in a different routine in mam4xx
VERT_EMISSION,
- //
+ // file with ncol, time
+ // example: surface emissions of aerosols and precursors
NONE
};Note I referred to the "VERT_EMISSION" files as they are actually named "elevated" (that's in the name of the files). I think it would have been good to keep that terminology from EAM. I believe the "NONE" type includes the surface emission. In the messy namelist world of EAM, we roughly call these "ext_frc_" (VERT_) and "srf_emis_" (NONE). So, I think you should actually change the "NONE" type to something more meaningful. I would encourage you to check with people like Balwinder, Kai, and Mingxuan.
There was a problem hiding this comment.
@mingxuanwupnnl explained very well why to use the term "external forcings" for emissions (as they are applied like forcing in the equations). I think NONE here means a tracer file with no vertical dimension (e.g., surface emissions). We can rename NONE to SURFACE, but NONE should also be fine as it would mean a tracer file with a NONE vertical extent.
There was a problem hiding this comment.
We can rename the enums. Note that we are not using the TracerFileType enum for surface emissions. The reader for surface emissions was implemented by @singhbalwinder. We worked in parallel on this, so we ended up with duplicate code. Because @bartgol is working on a general reader, it is better to just keep NONE, and I do not think we should add a reference for surface emissions. I will update the comments in this part of the code, as recommended by @mahf708.
@singhbalwinder, is this change okay: renaming VERT_EMISSION to ELEVATED_EMISSION? Or do we want to keep it as EXT_FRC?
There was a problem hiding this comment.
Let's change it to ELEVATED_EMISSIONS. I would like to keep a reference to "emissions" to be clear that these are actually emissions.
|
I have added a test mod that must be run with an NE30 grid. This test forces the NE30 grid to use NE4 emission files. A remapping file is provided that remaps NE4 emission files to NE30. We should probably revisit it later. I do not like restricting test mods to model resolution; they should be independent. |
Thank you for your detailed review and comments.
The single process test for microphysics and its baseline should catch DIFFs for the tracer reader. These tests were added in PR #3013 |
We are fixing a couple of issues encountered when running with ne120 for the vertical emission NetCDF files. Note that for resolutions finer than ne30, the case requires a remap file to perform horizontal interpolation. Our standalone test and CIME case did not catch this issue because, in these instances, we are not employing a remap file, so there is no horizontal interpolation.
The main issue was that the vertical emission files mam4xx is using have dimension names
altitudeandaltitude_intinstead oflevandilev. Thus, we encountered an IO error due to a mismatch between layouts.To fix this issue, we modified the source grid tag for LEV and ILEV using:
Details
@mahf708
The distinction among the file types FORMULA_PS, ZONAL, and VERT_EMISSION lies in their approach to vertical interpolation. The fourth enum, NONE, was introduced as a 'none' type. In EAM, these files utilize latitude and longitude to represent the horizontal dimensions. However, we opted to use ncremap to convert latitude/longitude files to the ne grid. More details can be found here
Therefore, while FORMULA_PS, VERT_EMISSION, and ZONAL types share the same method for horizontal interpolation, they employ distinct vertical interpolations. Although I examined how SPA performs vertical interpolation, I chose to directly port the vertical interpolation routines from EAM because these routines are compact and straightforward to transfer. Moreover, adopting the same vertical interpolation techniques in both EAM and EAMxx for mam4xx will simplify the code evaluation process.
FORMULA_PS
The files corresponding to FORMULA_PS include invariant tracers such as O3, OH, NO3, and HO2. For this type of file, the source pressure is computed using the following equation.
Vertical interpolation is performed using the ported vert_interp routine, available here.
The logic to detect this file type checks if the PS variable is present in the nc file; see here
This file is located here:
${DIN_LOC_ROOT}/atm/scream/mam4xx/invariants/ne4pg2/oxid_1.9x2.5_L26_1850-2015_ne4pg2_c20241009.nc
ncdump -h oxid_ne2np4_L26_1850-2015_c20240827.ncVERT_EMISSION
Vertical interpolation is performed by the ported 'rebin' routine: here
To select this file, the code looks for the 'altitude' dimension in the nc file; see here
This type of file corresponds to vertical emission files, where a file for each species is provided. For example:
${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne30pg2/elevated/cmip6_mam4_so2_elev_1x1_2010_clim_ne30pg2_c20241008.nc
ncdump -h cmip6_mam4_so2_elev_ne2np4_2010_clim_c20240726.ncZONAL
Vertical interpolation is performed using the ported vert_interp routine, available here.
The key difference between the FORMULA_PS and ZONAL types in our implementation is that the former does not utilize the PS variable to compute source pressure; instead, the source pressure is provided directly in the nc file.
I have named this file type 'zonal,' reflecting its usage in the EAM code as 'zonal_ave'. see here
This file is located here:
${DIN_LOC_ROOT}/atm/scream/mam4xx/linoz/ne4pg2/linoz1850-2015_2010JPL_CMIP6_10deg_58km_ne4pg2_c20240724.nc
ncdump -h linoz1850-2015_2010JPL_CMIP6_10deg_58km_ne2np4_c20240724.nc