atl_month_secs if off

[pfrater]

I believe that atl_month_secs has two digits rearranged. Currently it is:

mfdbatlantis:::atl_month_secs
[1] 2620800

whereas (year / month) (day sec):

(365 / 12) * (24*3600)
[1] 2628000

The current atlantis_fg_tracer which utilizes this to calculate year, month and day is showing weird spikes in the middle of certain years. Changing atl_month_secs to 2628000 fixes them to the appropriate recruitment months.

Also, I'm guessing that atl_month_days should be 30.41667 instead of 30.33333. The former produces 365 days in a year when multiplied by 12 months whereas the latter produces 364.

Thanks!

[pfrater]

Nevermind. I was not thinking like a programmer (i.e. 0-based indexing). However, there is still an issue with these spikes in numbers. They occur every 7.5 years because the atlantis output from two separate 30-day increments just happens to fall in the same month. But it is an issue because everything then appears twice as high in those particular months (e.g. landings, numbers of fish, etc.). Any thought to an elegant way to avoid this? I could just filter out these twice-in-a-month rows afterwards, too.

[lentinj]

Yes, this is something me and Chris talked about some time ago, Atlantis' calendar is a bit rough and ready: https://github.com/mareframe/mfdb/issues/20#issuecomment-144746321

IIRC after that many not-quite-months in seconds we run out of numerical precision to get the date exactly right. I have recollections of trying to work out a means of reconstructing days/months by assuming that time is going to progress in a uniform manner, but this was met with problems.

My model runs don't go very far into the future. What does the time dimension for your data look like, e.g.:

> nc_out <- ncdf4::nc_open('atlantis-Iceland-2016-05-24/OutCATCH.nc')
> cat(nc_out$dim$t$vals, "\n")

E-mail the values if it's huge.

[pfrater]

Time in Erla's model goes from 0 to 2049840000. She also says that the model is supposed to represent January 1, 1948 to January 1, 2013. The way mfdbatlantis is currently set up, though, it calculates the final month as 3/2013 (i.e. March, 2013). I'm not sure which way is correct. I think it partially depends on whether 86,400 seconds (3600*24) is defined as the end of the current day or the beginning of a new one.

The way it makes sense in my mind is that there are 3600*24 seconds in a day, and 365 of those in a year. The following code uses this line of thinking with integer division to get the year, then skim off the remainder and integer division to get months, then same for days. I think this is a bit different than how it's currently computed in atlantis_time_to_years() and the like.

days.in.year <- 365
sec.in.day <- 3600*24
sec.in.month <- sec.in.day*(days.in.year/12)
sec.in.year <- sec.in.day*days.in.year
td.seconds <- 3600*24*30
atl.final.time <- 2049840000
atl.time <- c(seq(0, 2049840000, by=td.seconds), atl.final.time)
origin.year <- 1948

year <- (atl.time %/% sec.in.year) + origin.year
month <- ((atl.time %% sec.in.year) %/% sec.in.month) + 1
day <- (((atl.time %% sec.in.year) %% sec.in.month) %/% sec.in.day) + 1

time <- data.frame(atl.time, year, month, day)
tail(time)

Changing days.in.year to 364 will show months 1,2,3 in 2013 in tail(time).

Getting the year and month is really the most important. The day doesn't matter, and I also think that, because of the way Atlantis outputs data (e.g. every 30 days), there will just be some months every few years in which 2 different days within the month get printed. Probably filtering these afterwards would be easiest unless you can think of a better way to do it during atlantis_fg_tracer().

[lentinj]

td.seconds <- 36002430 atl.final.time <- 2049840000 atl.time <- c(seq(0, 2049840000, by=td.seconds), atl.final.time)

Does this definitely generate the same sequence as Atlantis?

I think the core difference here is our definition of days in a month. You've got 365 / 12, I've got 30 + 1/3. I did this deliberately since this is what Atlantis seemed to be using (before I also assumed 365 / 12). Unfortunately the evidence of why, I can't find.

Whilst ditching duplicates is possible, the other problem is you need to keep in sync with Atlantis concept of seasons.

[pfrater]

It does generate the same sequence as Atlantis.

nc_out <- ncdf4::nc_open('OutCATCH.nc')
atl.output.time <- nc_out$dim$t$vals

td.seconds <- 3600*24*30
atl.final.time <- 2049840000
atl.time <- c(seq(0, atl.final.time, by=td.seconds), atl.final.time)

all(atl.output.time == atl.time)
[1] TRUE

Also, I'm not sure if you've noticed this, but the first two time instances output by atllantis_fg_tracer show month 1, day 1 (at atlantis time 0) and then month 2, day 31 (at atlantis time 2592000). After that it seems to get on the right track again with days being every 30 1/3.

[lentinj]

Okay, I'm either out-of-date or this is a bug in Atlantis depending on the output timestep you use. The model output I've got is aiming for every quarter:

> diff(ncdf4::nc_open('atlantis-Iceland-2016-05-24/Out.nc')$dim$t$vals)
 [1] 7862400 7862400 7862400 7862400 7862400 7862400 7862400 7862400 7862400 7862400 7862400 7862400 7862400 7862400 7862400 7862400 7862400 7862400
[19] 7862400 7862400  432000
> 60 * 60 * 24 * (30 + 1/3) * 3
[1] 7862400
> 60 * 60 * 24 * (365 / 12) * 3
[1] 7884000

I'll go back to 365 / 12 in this case.

Also, I'm not sure if you've noticed this, but the first two time instances output ...

Yes I had, but it's not quite as critical. Obviously still needs to be fixed mind :)

[lentinj]

Aha, this is starting to make more sense now, sorry it's taken me a while to remember all this. In the model output I have, Run.xml sets toutinc to "91 day" (~a quarter). This will be where 30 + 1/3 came from. I'm guessing in your model it's "30 day". Either way it's pretty much impossible to have readings that line up with a year/month calendar.

How about we use the following instead?

# Atlantis time units
atl_day_secs <- 60 * 60 * 24  # Seconds-in-day normal
atl_year_days <- 365  # NB: Atlantis years are 365 days, no leap years
atl_year_secs <- atl_year_days * atl_day_secs

# NB: The output we deal with is defined by toutinc, e.g. "30 day" or "91 day", so
# there can never be a neat monthly / yearly output as Atlantis years are 365 days.
# The output from Atlantis will drift, and there's not much that can be done about it.
# Unlike years, we are free to invent our own month definition, so assume whole 30 day months
# with a month 13 remainder of 5 days. Anything in month 13 can be filtered before import.
atl_month_days <- 30
atl_month_secs <- atl_month_days * atl_day_secs

# Convert vector of atl_time seconds to year
atlantis_time_to_years <- function (atl_time) {
    (atl_time %/% atl_year_secs)
}
# month vector from dims$time
atlantis_time_to_months <- function (atl_time) {
    ((atl_time %% atl_year_secs) %/% atl_month_secs) + 1
}
# day-in-month vector from dims$time
atlantis_time_to_days <- function (atl_time) {
    (((atl_time %% atl_year_secs) %% atl_month_secs) / atl_day_secs) + 1
}

atl.final.time <- 2049840000
td.seconds <- 3600*24*30 # i.e. seconds in 30 days
atl.time <- c(seq(0, 2049840000, by=td.seconds), atl.final.time)
rv <- data.frame(
    atl.time,
    m.seq = (seq(0, length(atl.time) - 1) %% 12) + 1,
    y = atlantis_time_to_years(atl.time),
    m = atlantis_time_to_months(atl.time),
    d = atlantis_time_to_days(atl.time))
head(rv, 50)

[pfrater]

This seems like a pretty clean and straightforward solution. Would it also make sense to filter out the 13th month within the functions that use these time conversions? Just thinking about future users wondering where the heck this 13th month came from. Maybe it doesn't matter...just a thought.

[lentinj]

Yes, filtering would be the sensible default. I'll finish up the changes today.

[lentinj]

https://github.com/mareframe/mfdbatlantis/commit/1b088585816d68f4d9daa903bbaf4301e1c99e1a and https://github.com/mareframe/mfdbatlantis/commit/59545181554de464509dbaca5988d93707000f9e should fix this. Let me know if I've slipped up.

[pfrater]

This looks good. Thanks!!