[mpas-developers] output bug in trunk

[Michael Duda]

Hi, Folks.

Doug and I have both (independently) uncovered an issue with output in
the trunk code. After some investigation, it appears that the issue is
related to how we track the fields to be written to output files (or
even read from input files, if those input files are read periodically
throughout a model run).

In the new IO layer, we maintain a list of pointers to field types to be
read or written to a stream, and for fields with multiple time levels
(e.g., state variables), we can only maintain a pointer to a specific
time level for the field; for example, when adding xtime to an output
file, we pass a pointer to 

   state % time_levs(1) % state % xtime 

to the MPAS_streamAddField() routine. 

During a model simulation, we shift the time levels for state variables
with a call to mpas_shift_time_levels_state(), which currently swaps
pointers for the field groups:

      sptr => state % time_levs(1) % state
      do i=1,state % nTimeLevels-1
         state % time_levs(i) % state => state % time_levs(i+1) % state
      end do
      state % time_levs(state % nTimeLevels) % state => sptr

After such a call, if the field pointer maintained in the IO layer
pointed to time_levs(1) % state % xtime, it would then point to
time_levs(2) % state % xtime, and writing output at that point would
then write the wrong time level for state fields. So, writing output at
an interval (in time steps) that is evenly divided by the number of time
levels would not reveal the problem; it's only when the output interval
is not evenly divided by the number of time levels that the problem
shows up. Also, since we re-add fields to an output stream when we close
one output file and open a new output file (in the case where, e.g.,
config_frames_per_outfile = 1), we add a new field pointer; so, the
issue doesn't appear when we write a single frame per history file.

To correct the problem, I think we just need to perform the time-level
shift beneath the field level; we can do this by modifying the registry
program to generate code for the mpas_shift_time_levels_state() routine
that swaps array pointers within each field in the state group:

      char0d = state % time_levs(1) % state % xtime % scalar
      do i=1,state % nTimeLevels-1
         state % time_levs(i) % state % xtime % scalar = state % time_levs(i+1) % state % xtime % scalar
      end do
      state % time_levs(state % nTimeLevels) % state % xtime % scalar = char0d

      real2d => state % time_levs(1) % state % u % array
      do i=1,state % nTimeLevels-1
         state % time_levs(i) % state % u % array => state % time_levs(i+1) % state % u % array
      end do
      state % time_levs(state % nTimeLevels) % state % u % array=> real2d

      (and so on for all fields in state...)

This way, the actual data pointed to by the field pointer in the IO
layer are updated with each time level shift, and the change is
transparent to the rest of the model. I've done some testing with the
shallow water model, and for an output interval of an even number of
time steps, the model results match bit-for-bit (as we would expect),
and with odd interval output (and multiple frames per file), the fields
are written from the correct time slice.

The extra work in the revised mpas_shift_time_levels_state() routine is
proportional to the number of fields in the state group, but the time
spent in this routine is very small in an absolute sense: on my desktop,
the routine routine takes about 4e-5 seconds for the shallow-water core,
which has 23 state arrays.

If it sounds reasonable to everyone, I'd like to change the registry
code that generates the contents of the mpas_shift_time_levels_state()
to swap pointers beneath the field level later this afternoon;
otherwise, any comments or suggestions would be welcomed.

Thanks,
Michael