Parallel I/O output module enhancements

[SeanBryan51]

---

📚 Documentation preview 📚: https://cable--655.org.readthedocs.build/en/655/

[SeanBryan51]

Hi @Whyborn, I've opened this PR to highlight the output module specific diffs (includes the aggregators and the list data structure for output variables). The implementation still contains a bunch of TODO's, but it would be awesome if I could get your feedback at this stage.

I've been testing it only in serial (testing in parallel will only be possible once the met data can be read in via parallel I/O) and currently gives the same answers as before for Qh.

However, I did have to change the NetCDF format in the old output module from classic to NetCDF4/HDF5 to easily compare outputs between new and old. This was a hack as the new NetCDF API has the NetCDF4 format hardcoded whenever creating or opening new files. I'll make that configurable so that both classic and NetCDF4 is supported.

[Whyborn]

First pass comments:

• I think the output should be dropped from cable_output_aggregator_t- the aggregator works exactly perfectly fine outside the output module, and can be used any time you need temporal aggregation for output or science.
• I don't think the output variable should handle the aggregation frequency. I think the control should be with the relevant science (I guess particularly the relevant science output module) module as to when an aggregator should be accumulated.
• What's the benefit to having the active argument in the cable_output_add_variable call, over an if statement?
• I think the general approach should be to have the aggregators defined first, and then pass as arguments to the cable_output_add_variable call. Something like:

class(aggregator_t), dimension(:), allocatable :: Qh_aggregators
...
Qh_aggregators = create_aggregators(canopy%fh, ["mean", "max"])
call cable_output_add_variable(
  ...
  aggregator=Qh_aggregators
  ...
)

This way it makes the aggregators available for use where necessary (e.g. the Tmx) variable. Although I'm not sure if you could make the elements of the returned Qh_aggregators targets?

• I don't think the shapes should be pre-defined. I think they should be built by the caller by passing dimensions, something like:

cable_output_shape_land_patch = create_decomposition(["land", "patch"])
cable_output_shape_land_patch_soil = create_decomposition(["land", "patch", "soil"])

where the dimensions are defined earlier in initialisation with something like:

call define_cable_dimension("soil", 6)
call define_cable_dimension("rad", 2)

I think this makes the design more easily extensible? This makes it trivial for someone developing new code which may have new dimensions to output their own stuff.

• In aggregators.f90, what's the point of the aggregator_store_t class? This seems like it could be

type aggregator_store_t
  class(aggregator_t), allocatable :: aggregators(:)
  integer :: num_aggregators
end type aggregator_store_t

where num_aggregators functions in the same way the current num_aggregators functions.

• I do think we need a time module, but we may be able to import that from elsewhere for a full featured timing module.

[SeanBryan51]

Thanks @Whyborn for your thoughts!

• I think the output should be dropped from cable_output_aggregator_t- the aggregator works exactly perfectly fine outside the output module, and can be used any time you need temporal aggregation for output or science.
• I don't think the output variable should handle the aggregation frequency. I think the control should be with the relevant science (I guess particularly the relevant science output module) module as to when an aggregator should be accumulated.

I agree that it would be nice if time related control were handled at a higher level (for example with the relevant science output module as you say). However, it wasn't obvious to me how this approach could allow for driving the intermediate aggregators required for Tmx and Tmn, where the intermediate aggregator accumulates at all frequencies and aggregates at daily frequencies, and the resulting aggregator accumulates at daily frequencies and aggregates at monthly frequencies. The cable_output_aggregator_t, along with the handling of the aggregation and accumulation frequencies at the aggregator level, was introduced so that aggregators can be driven independently at any accumulation or aggregation frequency which seemed easier for supporting the intermediate aggregator case. There might be a better solution here though, so I'm definitely more than happy to discuss further.

• What's the benefit to having the active argument in the cable_output_add_variable call, over an if statement?

For two reasons:

1. I think it's nicer for enabling variables depending on some sort of boolean logic. For example for groups of variables, we could just .or. the variable specific logical with the group logical rather than updating the output inclusion type instance as is done here.
2. The advantage of using active over nesting the call within if (active) is that this allows us to build a list of all possible supported outputs (both active and inactive) which could used to generate documentation on the supported outputs. This was inspired from the CTSM/CLM folks.

• I think the general approach should be to have the aggregators defined first, and then pass as arguments to the cable_output_add_variable call. Something like:

class(aggregator_t), dimension(:), allocatable :: Qh_aggregators
...
Qh_aggregators = create_aggregators(canopy%fh, ["mean", "max"])
call cable_output_add_variable(
  ...
  aggregator=Qh_aggregators
  ...
)

This way it makes the aggregators available for use where necessary (e.g. the Tmx) variable. Although I'm not sure if you could make the elements of the returned Qh_aggregators targets?

Oh yep I remember you mentioning this in our discussions. I think this is something we could definitely add in the future - I was hesitant to introduce this now as I want to avoid diverging in functionality from the current output module which assumes a single aggregation method per variable.

As for a list of objects being targets, I got around that problem by working with arrays of type aggregator_handle_t which contain a reference to an aggregator rather than the actual aggregator instance.

• I don't think the shapes should be pre-defined. I think they should be built by the caller by passing dimensions, something like:

cable_output_shape_land_patch = create_decomposition(["land", "patch"])
cable_output_shape_land_patch_soil = create_decomposition(["land", "patch", "soil"])

where the dimensions are defined earlier in initialisation with something like:

call define_cable_dimension("soil", 6)
call define_cable_dimension("rad", 2)

I think this makes the design more easily extensible? This makes it trivial for someone developing new code which may have new dimensions to output their own stuff.

I like this suggestion, I agree it's definitely not that obvious what one needs to do to create a new CABLE_OUTPUT_SHAPE_TYPE_*. I think rather than determining the shape of the decomposition object from the defined dimensions, I would just initialise a new decomposition via the io_decomp_* procedures, or use the io_decomp instance. Removing the SHAPE_TYPE_* variables means I will need to rethink about how each variable will infer the temporary buffer used for grid cell averaging. I'll reflect on that a bit.

• In aggregators.f90, what's the point of the aggregator_store_t class? This seems like it could be

type aggregator_store_t
  class(aggregator_t), allocatable :: aggregators(:)
  integer :: num_aggregators
end type aggregator_store_t

where num_aggregators functions in the same way the current num_aggregators functions.

aggregator_store_t was introduced because allocatable arrays of a polymorphic type are problematic (i.e. class(aggregator_t), allocatable :: aggregators(:)), as the compiler doesn't know in general what the type (and hence size) of each element will be on allocation. In practice, you need a container derived type which itself is not polymorphic, like aggregator_store_t, as the list element. It's similar to how you would declare an array of pointers in fortran. Interestingly the IBM compiler seems to actually allow for allocatable arrays of a polymorphic type: https://www.ibm.com/docs/en/xl-fortran-aix/16.1.0?topic=concepts-array-constructors, but its not standard as far as I know.

• I do think we need a time module, but we may be able to import that from elsewhere for a full featured timing module.

I agree, a timing module for wider use in the code would be great. The time_step_matches function was my hacky attempt to clean up the timing logic code in write_output, definitely something we could introduce properly (independent of these changes).

[Whyborn]

However, it wasn't obvious to me how this approach could allow for driving the intermediate aggregators required for Tmx and Tmn, where the intermediate aggregator accumulates at all frequencies and aggregates at daily frequencies, and the resulting aggregator accumulates at daily frequencies and aggregates at monthly frequencies. The cable_output_aggregator_t, along with the handling of the aggregation and accumulation frequencies at the aggregator level, was introduced so that aggregators can be driven independently at any accumulation or aggregation frequency which seemed easier for supporting the intermediate aggregator case.

Does this allow aggregators to be driven independently? It seems to me like specifically doesn't allow this, because every aggregator which accumulates e.g. on_timestep, daily etc would be accumulated at the same time. This would mess with things that happen at the end of the day, using some aggregation of on_timestep quantities. I think the Tmn and Tmx variables are examples- I'm fairly sure these feed into CASA science. That would mean, to get the correct daily minima and maxima, the order of operations would have to look like:

do_biogeophysics -> do_timestep_accumulation -> do_CASA -> do_daily_accumulation

There has to be a way to trigger accumulation at specific times, for instances like this. This is why I think the aggregators have to be available to work with as standalone objects.

I think it's nicer for enabling variables depending on some sort of boolean logic. For example for groups of variables, we could just .or. the variable specific logical with the group logical rather than updating the output inclusion type instance as is done here.
The advantage of using active over nesting the call within if (active) is that this allows us to build a list of all possible supported outputs (both active and inactive) which could used to generate documentation on the supported outputs. This was inspired from the CTSM/CLM folks.

Yea I can see that, it might be more easily readable to have every call contain a

output%variable .or output%variable_group .or. output%variable_module

instead of a construction like

if (output%variable_module) then
  if (output%variable_group) then
    if (output%variable) then
      ...
    end 
  end
end

Although I don't see why the latter couldn't also be used to create the same documentation in the same way.

I like this suggestion, I agree it's definitely not that obvious what one needs to do to create a new CABLE_OUTPUT_SHAPE_TYPE_. I think rather than determining the shape of the decomposition object from the defined dimensions, I would just initialise a new decomposition via the io_decomp_* procedures, or use the io_decomp instance. Removing the SHAPE_TYPE_ variables means I will need to rethink about how each variable will infer the temporary buffer used for grid cell averaging. I'll reflect on that a bit.

Yea that's what I meant, whether the io_decomp is created via dimensions names or just numbers, which could be accessible with something like dim_length(dim_name), is much of the same to me. I do think it's important that the various dimensions can be accessed in this way, since it's something done already throughout the code via use statements, and this would help eliminate that method of variable sharing.

aggregator_store_t was introduced because allocatable arrays of a polymorphic type are problematic (i.e. class(aggregator_t), allocatable :: aggregators(:)), as the compiler doesn't know in general what the type (and hence size) of each element will be on allocation. In practice, you need a container derived type which itself is not polymorphic, like aggregator_store_t, as the list element. It's similar to how you would declare an array of pointers in fortran. Interestingly the IBM compiler seems to actually allow for allocatable arrays of a polymorphic type: https://www.ibm.com/docs/en/xl-fortran-aix/16.1.0?topic=concepts-array-constructors, but its not standard as far as I know.

Are you sure? I thought arrays of polymorphic classes were part of the standard, I used them in some of my aggregator testing. You just need select type blocks to access any of the child components.

I agree, a timing module for wider use in the code would be great. The time_step_matches function was my hacky attempt to clean up the timing logic code in write_output, definitely something we could introduce properly (independent of these changes).

There's the datetime-fortran which we already have a spack package for? I haven't actually looked at it's features yet though.

[SeanBryan51]

Hi @Whyborn, I've made the updates we discussed last week. I think it's looking much better now with your suggestions on how aggregators are organised in the output module, thanks as always for the comments! Please let me know if there is anything else that catches your eye on the design.

I'm going to try add support for specifying non-time varying data in the output and restart files. For this I'm thinking we could introduce a cable_output_add_parameter subroutine (similar to cable_output_add_variable) for specifying the non-time varying parameters in the output, and a similar cable_output_add_restart_variable for the restart variables.

[Whyborn]

Just a couple of comments:

• The cable_output_add_variable routine doesn't have an aggregation_frequency yet- I think we need that and a minimum frequency
• I think it would be good to decouple the output and the grid cell averaging (and reductions in general). A bit like the aggregators- it's a functionality useful for the output, but should be considered a distinct module that the output module uses.
• Along that line, I think it would be nice to define reducers/mapper by name. So you could define the variable like:

if (output%patch)
  decomp = decomp_land_real32
  reducer = "patch"
else
  decomp = decomp_land_patch_real32
  reducer = "none"
end if

call cable_output_add_variable(
  name="Qh",
  ...
  decomp=decomp,
  reducer=reducer
)

And then in the write_variable you could have:

subroutine write_variable(output_variable, output_file, time_index)
...
  if (output_variable%reducer == "patch") then
    call output_file%write_darray(
        var_name=output_variable%name,
        values=compute_patch_average(aggregator%storage)
        decomp=output_variable%decomp,
        frame=time_index
  else if (output_variable%reducer == "none") then
    call output_file%write_darray(
        var_name=output_variable%name,
        values=aggregator%storage
        decomp=output_variable%decomp,
        frame=time_index
  end if

Then you wouldn't have to have all the buffers. The compute_patch_average would be a function interface, with different rank/type module procedures.

• Can you remind me why the all the decompositions for different element types are required to have distinct variables? Aren't the decomps basically just indexing arrays?

[SeanBryan51]

• The cable_output_add_variable routine doesn't have an aggregation_frequency yet- I think we need that and a minimum frequency

I'm definitely happy to add the frequency limit to cable_output_add_variable as that information is no doubt specific to each variable. I removed the aggregation_frequency argument because the aggregation frequency is tied to the output frequency which is set at the "profile" level.

• I think it would be good to decouple the output and the grid cell averaging (and reductions in general). A bit like the aggregators- it's a functionality useful for the output, but should be considered a distinct module that the output module uses.

That sounds good to me, I will rename the module containing the grid cell averaging stuff to be more general (and maybe put this under utils/)?

• Can you remind me why the all the decompositions for different element types are required to have distinct variables? Aren't the decomps basically just indexing arrays?

The distinction of types for decompositions is required by PIO via the basepiotype argument of pio_initdecomp.

• Along that line, I think it would be nice to define reducers/mapper by name.

I'm happy to introduce a string argument for reducer instead of the grid_cell_average logical.

Then you wouldn't have to have all the buffers. The compute_patch_average would be a function interface, with different rank/type module procedures.

I did consider using a function interface instead of a subroutine interface, however I opted for the subroutine approach with the temporary buffer as this avoids introducing a potentially large allocation when computing the grid cell average. I want to limit the number of unnecessary allocations and copy operations in the write procedures where possible. It might be possible to return a preallocated array from a function. I can look into this a bit more.

Thank you again for the feedback on this!

[ccarouge]

If we read a met file during the initialisation (outside the time loop), we could then initialise "CABLE_USER%YearStart = met%year(1)" for the MetType = ''?
It would remove the need for this if condition and it makes sure that cable_user%YearStart is always defined and always has the same meaning.

Not sure that's what you mean in your TODO. It seems to apply more to the code in cable_timing_utils.F90 than here.

[SeanBryan51]

Yep I definitely want to get back to this, I'll test out the solution you propose. I mean to try out the site case (which I believe corresponds to MetType == ' ') to confirm I haven't broken things here for this case.

[SeanBryan51]

Not sure that's what you mean in your TODO. It seems to apply more to the code in cable_timing_utils.F90 than here.

You've interpreted that correctly, it was a reminder to myself to look into an alternative algorithm for monthly timings that doesn't require the start year of the simulation.

[ccarouge]

As an alternate, we could consider an additional "ktau"-like variable that is reset every year. (maybe there is one in the code, I don't know). The time loop tells us when we change year (to check if this is true with the site simulations).

[Whyborn]

What's the reason not to have cable_user%yearstart be something that is unconditionally required in the namelist, and have that be the final truth about the start point of the simulation? Having that able to be overwritten is bound to lead to confusion in the future.

[SeanBryan51]

As an alternate, we could consider an additional "ktau"-like variable that is reset every year. (maybe there is one in the code, I don't know). The time loop tells us when we change year (to check if this is true with the site simulations).

I've seen this pattern used in other models (see #656). I agree that tracking information like this throughout the run makes sense

[SeanBryan51]

What's the reason not to have cable_user%yearstart be something that is unconditionally required in the namelist, and have that be the final truth about the start point of the simulation? Having that able to be overwritten is bound to lead to confusion in the future.

That's a great point, it might make sense to instead error if met%year(1) does not match cable_user%yearstart, and rely on cable_user%yearstart for timing logic.

I don't really understand why met%year(1) is required for this case, I'll see what happens when we rely on cable_user%yearstart alone

[Whyborn]

Isn't ktau already that? The timestep within the current year, with ktau_tot being the total timestep in the simulation?

[SeanBryan51]

I believe ktau can extend beyond an year for site configurations (dependent on kend). For the MetType = 'gswp3' case that I'm testing at the moment, ktau is reset every year. This is probably a symptom of the many met forcing formats / configurations which are supported by the driver.

[Whyborn]

It looks like ktau loops with do ktau = kstart, kend here and kend is set to be the number of steps in a year in almost every instance, except for gswp which determines it based off the length of the time dimension in the met file (this case block)

[SeanBryan51]

Thanks for the comments @ccarouge, I'll add in those suggestions!

[SeanBryan51]

Just an update on this PR! I'm currently adding functionality for restarts to the new output module by adding a restart flag component to cable_output_variable_t. This will roll back the last two commits I made which introduce a separate module for handling restart files: src/offline/cable_serial.F90: write restarts via new restart modules and Add cable_restart_mod and cable_restart_write_mod. There are a few things in the works to get this implemented, namely:

1. Currently the new output module only allows for writing distributed arrays so I'm currently adding support for writing non-distributed (global) arrays as this is required for some restart fields.
2. Quality of life change: instead of specifying the netCDF variable dimensions when declaring output variables, specify the dimensions of the in-memory array and instead infer the netCDF dimensions from these - this is more meaningful as the netCDF dimensions tend to vary based on the configuration, e.g. whether individual patches are written or the x-y vs land grid format (only temporary of course). This also simplifies the API for declaring output variables (see next point).
3. Remove the decomp and temp_buffer_* components from cable_output_variable_t and instead infer these at write time - this removes some data redundancy as the reduction_method and the in-memory array determine decomp and temp_buffer_*. This also simplifies the declaration of output variables in the code as the decomposition also depends on configuration options - accounting for these options makes the code for declaring output variables unnecessarily bloated.

More to come 🙂

[SeanBryan51]

Hi @Whyborn, apologies it has been a while since I've followed up on this! I had to stop myself from tweaking small things 🙃

The output module now has functionality for writing restarts and has gone through a pretty aggressive restructuring. Since a lot of the previous commits on the output module are no longer relevant with the restructuring, I've blown away most of commit history and grouped the commits into these categories:

• Utilities (e.g. enum type, aggregator type)
• NetCDF API enhancements
• Non-output module related bug fixes or comments
• Output module implementation including restart file stuff

Please let me know what you think! (at a time that's convenient of course, it's a pretty dense PR 😅 )

The next thing on the plate after this is would be to introduce working variables for all diagnostics, and then add all output and restart variables and test for (approximate) bitwise reproducibility.