Clean up get_grid

[mcgibbon]

There are some points we should talk about in get_grid:

• Should x_name and y_name be inputs? If those are changed, any components in CliMT that use lat or lon dimensions will break (unexpectedly). I feel like if the user wants to change those for diagnostic reasons, they should do it at the time that they output a netCDF file.
• What is the best way to supply top and bottom pressure to get_grid? The two candidates seem to be setting constants in the constants library, or kwargs. I feel like setting constants is a convoluted way to go about this. There's no reason why the surface pressure should be the reference pressure. And I don't currently understand why model top pressure must be read by a component as a Sympl constant - can't the component just read air_pressure_on_interface_levels to determine model top pressure? That way would ensure we don't unintentionally develop in consistencies between the model top pressure in the state and the Sympl constant.

@JoyMonteiro also brought up a point about e.g. ground schemes wanting to be on separate grids from the atmosphere. The way CliMT is currently set up, I would see that being done by get_grid, or by splitting get_grid into multiple functions (get_atmosphere_grid, get_land_grid, get_ocean_grid, etc.). This also requires couplers. What do you think of this, @JoyMonteiro ?

[JoyMonteiro]

I'm fine with getting rid of x/y_name for now. We don't have a use case where this is not applicable currently, and so this is not something we need to support in the near future at least.

model top pressure is decided a priori, and air_pressure_on_interface_levels must describe a pressure
distribution that conforms to this pressure. Hence, this is more fundamental to the model description.

The algorithm to calculate the vertical coordinates requires a reference pressure (not necessarily the actual surface pressure) to calculate the coordinates. You will need the actual surface pressure to calculate air_pressure_on_interface_levels, but that is not a part of the grid.

We will need to put in some thought into the part about getting different grids. Will reply back once we finish submitting the paper ;) I don't think I have the bandwidth to accomplish both in the time remaining.

[mcgibbon]

This is how I think model top pressure should be defined:

• get_grid should be told model top pressure, and create air_pressure arrays that conform to this.
• When GFS is first called (the only component that cares about model top pressure), it remembers the model top pressure, and errors if the model top pressure is not uniform.
• On successive time steps, it errors if the model top pressure has changed. All of this erroring will already happen automatically if the GFS code is passed in the initial model top pressure when it makes its version of the pressure arrays that we check the state pressure arrays against (because those arrays won't match).

This completely avoids the need for a Sympl constant which may or may not actually be in agreement with the state.

[JoyMonteiro]

Nope, any component that modifies the surface pressure should ensure that the other pressure arrays are consistently changed -- DcmipInitialConditions also needs model top pressure. If model top pressure is not available easily to components, then keeping things consistent becomes a mess -- you will need the user to remember which components require model top pressure, which is really unnecessary.

When GFS is first called, for instance, there is no guarantee that the arrays have not been already changed to make it inconsistent. There needs to be a reliable way to access such an important constant.

[mcgibbon]

Look at the modifications I have made to GFS. It does not need the model top or surface pressures until it is called the first time. At that time, it gets that information from the state. The user does not need to keep track of anything. Only get_grid needs the pressures. We can modify Dcmip to work the same way. I can work on this tomorrow if you like, going to bed right now. I'm trying to avoid duplicating information. If we duplicate, we need to do a lot of work to make sure the two places this info is stored don't become inconsistent. I think we can do everything we need while only storing the info in the model state.

…

On Mon, Aug 13, 2018, 12:20 AM JoyMonteiro ***@***.***> wrote: Nope, any component that modifies the surface pressure should ensure that the other pressure arrays are consistently changed -- DcmipInitialConditions also needs model top pressure. If model top pressure is not available easily to components, then keeping things consistent becomes a mess -- you will need the user to remember which components require model top pressure, which is really unnecessary. When GFS is first called, for instance, there is no guarantee that the arrays have not been already changed to make it inconsistent. There needs to be a reliable way to access such an important constant. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#66 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ALvMhWNECZTQqMzFbHydnoxltDzAj43wks5uQSiAgaJpZM4V5QR4> .

[JoyMonteiro]

Yes, I have looked through the GFS code. What I'm pointing out is that if the pressure arrays are modified
by hand before calling GFS and the model top pressure is changed (even by error), then inferring the model top pressure from these arrays is a bad thing to do. If, on the other hand, the model top pressure is available from a reliable source, like the constants dictionary, then GFS can verify that the pressure arrays are consistent with what it calculates.

Sure, we can chat about this once you are awake, but we need to converge on a solution soon.

[mcgibbon]

If the model top pressure is changed before GFS is called and not changed after, then there's no problem. I'll lay out what I'm thinking in concrete terms. This is what I envision:

1. You initialize any components that require model top pressure. Instead of actually running their initialization, they set self._initialized=False. You can't run their initialization code within __init__ because the model top pressure hasn't been defined yet.
2. You initialize the state with get_grid or get_default_state, passing model top pressure as a kwarg. This is the only place where the user sets model top pressure by hand.
3. You run your model the way you normally do.
4. On the first time the aforementioned components are called, they determine model top pressure from air_pressure_on_interface_levels, and store it in self._model_top_pressure. They run their initialization code and set self._initialzed=True before doing computation.
5. On calls after the first time, they determine model top pressure from air_pressure_on_interface_levels and compare it to self._model_top_pressure. If an inconsistency has arisen, the component updates its internal representation of model top pressure if it can, or raises an exception if it can't (such as in GFS where the initialization can't be re-run with a new model top pressure).

If the user changes model top pressure outside of get_grid, one of two things will happen:

• If it happens before 4), nothing bad happens. There's no issue with changing the model top pressure at that point, because nothing has stored it internally.
  If it happens after 4), the component which has stored the pressure will raise an exception, and the user will fix their code.

I see the above as being pretty air tight against inconsistency. By comparison, this is how I see the constants dictionary approach:

1. You set the model top pressure as a constant.
2. You initialize components. They store the constant model top pressure.
3. You initialize the state, which uses the stored constant model top pressure.
4. You run the model the way you normally do.
5. You could augment this by having the components compare the state model top pressure to its internally stored model top pressure.

This is similarly tight against inconsistency between the state and the components, where you get exceptions only if you change the model top pressure constant between 2) and 3) (as long as we use step 5). However, there's an additional place the model top pressure is stored - the constants dictionary. If you change it any time after step 3, the constants dictionary model top pressure will be inconsistent with both the state model top pressure (which matters) and the component model top pressure (which doesn't matter as much). The user may believe they've modified the state by setting this constant, when they have done no such thing.

It's also uncertain what should actually happen when you change the constants dictionary model top pressure during runtime. Does this mean the defintion of the pressure grid has changed and that the hybrid A and B are actually still valid, so that the state dictionary is inconsistent with itself (assuming it has air pressures stored)? Should GFS not have an internal model top pressure and instead use the updated constant dictionary pressure, and change the pressure arrays when it's passed the state? Or error if the state pressure arrays are inconsistent with the A and B? On top of the potential for inconsistency, I don't like the ambiguities in this way of doing it.

On top of that, I think it's very roundabout to pass in configuration of the grid to the constants dictionary, instead of to get_grid or get_default_state.

[mcgibbon]

In the end it's a matter of separation of concerns. The state is responsible for storing information about model state (including the grid), and model top pressure is information about the state. Since we can write this to run consistently while only storing model top pressure in state, I think we should.

[JoyMonteiro]

I'm coming around to your point of view that storing the information in the constants dictionary may not be the best approach.

Come to think of it, I think this discussion was not necessary -- a more reliable source of the model top pressure is the hybrid coordinate pair a and b -- the model top pressure is given by Eq. 7 in the Eckermann paper when eta = 0, which is simply a[-1] + b[-1]*reference_pressure. For numerical reasons you want b[-1] = 0, so a[-1] gives us the model top pressure. Somebody really needs to search for and modify this variable, so we can assume they know what they are doing.

This gives the advantage that components don't need to store model top pressure at any point in time, and state has all the information to verify its own consistency.

So, shall we do this? use the coordinates to infer model top pressure for calculations with GFS and Dcmip (for now)?

[mcgibbon]

Currently GFS will use the coordinates to infer the full pressure profile, and will check the pressure profile in the state against its inferred pressure profile. Dcmip should probably just use air_pressure from the state. Since it's only modifying initial conditions, it shouldn't care whether the model top pressure is constant (whether b[-1] = 0)

It may still be necessary for GFS to check that b[-1] = 0 and that a[-1] is the same as it was initially. Since GFS checks that the state pressure is consistent with the hybrid coordinates, it doesn't matter whether it's checking against air_pressure_on_interface_levels or a and b.

[JoyMonteiro]

not sure what you mean, but if dcmip uses air_pressure from state to create its output arrays while modifying surface_air_pressure, these arrays will be inconsistent and will error in GFS (which is why those tests were failing).

GFS can check whether the whole of a and b are the same to ensure consistency, but that is a separate issue. Since it uses a grid defined by a, b rather than by pressure, it seems technically more appropriate to use the former two rather than the latter. This would then be an appropriate template when writing components which use height as a coordinate, for example.

At this point, I'm fine either way. If in the future I see an potential issue, we can take it up then.

[mcgibbon]

OK, it sounds like Dcmip needs to modify air_pressure and air_pressure_on_interface_levels when it updates surface_air_pressure, and should use A and B to do so.

You could add an additional test in GFS to ensure that A and B are the same between time steps. I would still keep the pressure check at each time step (it ensures no other components have made air pressure inconsistent).