3.1 Monthly temperature

The impacts of using the appropriate calendar to summarize the data (as opposed to not) are large, often exceeding 1 ^∘C in absolute value (Fig. 11). The effects are spatially variable and are not simple functions of latitude, as might initially be expected, because the effect increases with the amplitude of the annual cycle (which has a substantial longitudinal component) for temperature regimes that are in phase with the annual cycle of insolation. For temperature regimes that are out of phase with insolation, the calendar-adjusted minus unadjusted values would be negative and largest when the temperature variations were exactly out of phase. (If there were no annual cycle, i.e., if a climate variable remained constant over the course of a year, the calendar effect would be zero.) The interaction between the annual cycle and the direct calendar effect on insolation produces patterns of the overall calendar effect that happen to resemble some of the large-scale responses that are frequently found in climate simulations, both past and future, such as high-latitude amplification or damping, continent–ocean contrasts, interhemispheric contrasts, and changes in the seasonality of temperature (see Izumi et al., 2013). Because the month-length calculations use the Northern Hemisphere vernal equinox as a fixed origin for the location of Earth along its orbit, the effects seem to be small during the months surrounding the equinox (i.e., February through April; Fig. 11), and indeed the selection of a different origin would produce different apparent effects (see Joussaume and Braconnot, 1997; Sect. 2.1). However, the selection of a different origin would not change the relative (to present) length of time it would take Earth to transit any particular angular segment of its orbit.

At 6 ka, the largest calendar effects on temperature can be observed over the Northern Hemisphere continents for the months from September through December (Fig. 11), consistent with the earlier beginning of these months (Fig. 2) and the direct calendar effect on insolation at 45^∘ N (Fig. 3). For example, in the interior of the northern continents, as well as North Africa, temperature is in phase with insolation, and therefore the calendar effect on insolation (Fig. 3), which produces strongly positive differences from August through November, is reflected by the calendar effect on temperature. Over the northern oceans, temperature is broadly in phase with insolation but with a lag, which reduces the magnitude of the effect and gives rise to an apparent land–ocean contrast that otherwise might be interpreted in terms of some particular paleoclimatic mechanism. The calendar effect on temperature from January through March produces negative calendar-adjusted minus unadjusted values in the northern continental interiors (Fig. 11), which is also consistent with the calendar effect on insolation. In the Southern Hemisphere at 6 ka, the calendar effects on temperature produce generally negative differences, which is consistent with the calendar effects on mid-month insolation at 45^∘ S (Fig. 5) that produce generally negative differences throughout the year, particularly during the months of August through November. Like the continent–ocean contrast in the Northern Hemisphere, the Northern Hemisphere–Southern Hemisphere contrast in the calendar effect on temperature could also be interpreted in terms of one or another of the mechanisms thought to be responsible for interhemispheric temperature contrasts.

At 127 ka, the calendar effect on temperature is broadly similar to that at 6 ka over the months from September through March, but it differs in sign from April through July and in magnitude in August (Fig. 11). These patterns are also consistent with the direct calendar effects on insolation. At 127 ka, the calendar effect on insolation produces strongly positive differences in the Northern Hemisphere earlier in the northern summer than at 6 ka (Fig. 3), while at 45^∘ S the calendar effect on insolation produces strongly negative differences in July and persists that way through November (Fig. 5). At 116 ka, perihelion occurs in late December in comparison to late June at 127 ka (Figs. 1 and 6), and not surprisingly the calendar effect on temperature is nearly the inverse of that at 127 ka (Fig. 11). This pattern has important implications for paleoclimatic studies because in addition to all of the changes in the forcing and the paleoclimatic responses accompanying the transition out of the last interglacial, the possibility that some of the apparent simulated changes between 127 and 116 ka may be an artifact of data analysis procedures cannot be discounted.

At 97 ka, a time selected to illustrate a different orbital configuration (i.e., one with boreal spring months occurring closer to the June solstice) than the similar (6 and 127 ka) or contrasting (127 and 116 ka) configurations, the calendar effect on temperature in the Northern Hemisphere (Fig. 11) shows a switch from positive differences in the early boreal summer (May and June) to negative in the late summer (August and September). This switch is again consistent with the direct calendar effect on insolation (Fig. 3). Like the other times, these spatial variations in the calendar effect could easily be interpreted in terms of one kind of paleoclimatic mechanism or another.

The generally larger calendar effect on temperature over the continents than over the oceans implicates the amplitude of the seasonal cycle in the size of the effect. This situation suggests that even in model-only intercomparisons (and even in the unlikely case that all models involved in an intercomparison use the same calendar) the calendar effect could be present because the amplitude of the seasonal cycle is dependent on model spatial resolution (and its influence on model orography).

3.2 Mean temperature of the warmest and coldest months

Although the calendar effects on monthly mean temperature show some subcontinental-scale variability, the overall patterns are of relatively large spatial scales and are interpretable in terms of the direct orbital effects on month lengths and insolation. The calendar effects on the mean temperature of the warmest (MTWA) and coldest (MTCO) calendar months (and their differences) are much more spatially variable (Fig. 12). This variability arises in large part because of the way these variables are usually defined (e.g., as the mean temperature of the warmest or coldest conventionally defined month as opposed to the temperature of the warmest or coldest 30 d interval), but also because the calendar adjustment can result in a change in the specific month that is warmest or coldest. These effects are compounded when calculating seasonality (as MTWA minus MTCO). Other definitions of the warmest and coldest month are possible, such as the warmest consecutive 30 d period during the year (e.g., Caley et al., 2014), and such definitions will not be susceptible to the calendar effect. In practice, however, paleoclimatic reconstructions based on calibrations or forward-model simulations routinely use conventional calendar-month definitions of the warmest and coldest months and of seasonality (Bartlein et al., 2011; Harrison et al., 2014), and often only monthly output from paleoclimatic simulations is available, necessitating consistent definitions when summarizing model output.

In the particular set of example times chosen here, the magnitudes of the calendar effects are also smaller than those of individual months because, as it happens, the calendar effects in January and February (typically the coldest months in the Northern Hemisphere) and July and August (typically the warmest months in the Northern Hemisphere) are not large. There are also some surprising patterns. The inverse relationship between the calendar effects at 116 and 127 ka that might be expected from inspection of the monthly effects (Fig. 11) are not present, while the calendar effects on MTCO and MTWA at 97 and 116 ka tend to resemble one another (Fig. 12). Across the four example times, there is an indistinct but still noticeable pattern in reduced seasonality (MTWA minus MTCO) between the adjusted and unadjusted values, which like the other patterns described above could tempt interpretation in terms of some specific climatic mechanisms.

3.3 Monthly precipitation

In contrast to the large spatial-scale patterns of the calendar effect on temperature, the patterns of the calendar effect on precipitation rate are much more complex, showing both continental-scale patterns (like those for temperature) and smaller-scale patterns that are apparently related to precipitation associated with the ITCZ and regional and global monsoons (Fig. 13). The continental-scale patterns are evident in the calendar effects at 6 and 127 ka, particularly in the months from September through November (Fig. 13), when it can also be noted (especially over the midlatitude continents in both hemispheres) that there is a positive association with the calendar effect on temperature. This association is simply related to similarities in the shapes of the annual cycles of those variables and not to some kind of more elaborate thermodynamic constraint. At 116 ka, as for temperature, the large-scale calendar-effect patterns appear to be nearly the inverse of those at 127 ka. The smaller-scale kind of pattern is well illustrated at 127 ka in the tropical North Atlantic, sub-Saharan Africa, and south Asia. There, negative calendar-adjusted minus unadjusted values can be noted for June through August, giving way to positive differences from September through November, and the same transition appears inversely at 116 ka. Another example can be found in the South Pacific Convergence Zone in austral spring and early summer (September through November) at 6 and 127 ka, when generally positive differences between calendar-adjusted and unadjusted values in July and August give way to negative differences from September through December. This second kind of pattern, most evident in the subtropics, is not mirrored by the calendar effects on temperature.

Overall, the magnitude and spatial patterns of the calendar effects on temperature and precipitation (Figs. 11 and 13) resemble those in the paleoclimatic simulations and observations that we attempt to explain in mechanistic terms (Harrison et al., 2016). Depending on the sign of the effect, neglecting to account for the calendar effects could spuriously amplify some signals in long-term mean differences between the experiment and control simulations, while damping others.

3.4 Calendar effects and transient experiments

Calendar effects must also be considered in the analysis of transient climate-model simulations (even if those data are available on the daily time step). This can be illustrated for a variety of variables and regions using data from the TraCE-21ka transient simulations (Liu et al., 2009; https://www.earthsystemgrid.org/project/trace.html, last access: 11 August 2019). The series plotted in Fig. 14 are area averages for individual months on a yearly time step, with 100-year (window half-width) locally weighted regression curves added to emphasize century-timescale variations. The original yearly time-step data were aggregated using a perpetual noleap (365 d) calendar (using the present-day month lengths for all years). The gray and black curves in Fig. 14 show these unadjusted “original” values, while the colored curves show month-length adjusted values (i.e., pseudo-daily interpolated values, reaggregated using the appropriate paleo fixed-angular calendar). Area averages were calculated for ice-free land points.

Figure 14Time series of original and month-length adjusted annual area-weighted averages of TraCE-21ka data (Liu et al., 2009), expressed as differences from the 1961–1989 long-term mean for (a–c) 2 m air temperature, (d) precipitation rate, and (e–f) precipitation minus evaporation (P−E). The original or unadjusted data are plotted in gray and black, and the adjusted data are in color. The area averages are grid-cell area-weighted values for land grid points in each region, and the smoother curves are locally weighted regression curves with a window half-width of 100 years. The regions are defined as (a) 15 to 75^∘ N and −170 to 60^∘ E, (b) 10 to 50^∘ S and 110 to 160^∘ E, (c) global ice-free land area, (d) 0 to 30^∘ N and −30 to 120^∘ E, (e) 5 to 17^∘ N and −5 to 30^∘ E, and (f) 31 to 43^∘ N and −5 to 30^∘ E.

Download

Figure 14a shows area-weighted averages for 2 m air temperature for a region that spans 15 to 75^∘ N and −170 to 60^∘ E, the region used by Marsicek et al. (2018) to discuss Holocene temperature trends in simulations and reconstructions. The largest differences between month-length adjusted values and unadjusted values occur in October between 14 and 6 ka, when perihelion occurred during the northern summer months. October month lengths during this interval were generally within 1 d of those at present (Fig. 1), but the generally shorter months from April through September resulted in Octobers beginning up to 10 d earlier in the calendar than at present, i.e., closer in time to the boreal summer solstice (Fig. 2). The calendar-effect-adjusted October values therefore average up to 4 ^∘C higher than the unadjusted values during this interval (Fig. 14a), consistent with the direct calendar effects on insolation at 45^∘ N (Fig. 3). The calendar effect also changes the shape of the temporal trends in the data, particularly during the Holocene. October temperatures in the unadjusted data showed a generally increasing trend over the Holocene (i.e., since 11.7 ka), reaching a maximum around 3 ka, comparable with present-day values, while the adjusted data reached levels consistently above present-day values by 7.5 ka. The unadjusted October temperature data could be described as reaching a “Holocene thermal maximum” only in the late Holocene (i.e., after 4 ka), while the adjusted data display more of a mid-Holocene maximum. As is the case with the mapped assessments of the pure calendar effect, the differences between unadjusted and adjusted time series are of the kind that could be interpreted in terms of various hypothetical mechanisms. For example, the calendar-effect adjustment advances the time of the occurrence of a Holocene thermal maximum in October by about 3 kyr for North America and Europe.

As in North America and Europe, the adjusted temperature trends in Australia (10 to 50^∘ S and 110 to 160^∘ E) (Fig. 14b) are consistent with the direct calendar effects on insolation (i.e., for 45^∘ S; Fig. 5). The difference between adjusted and unadjusted values are again largest in October between 14 and 6 ka, but the difference is the inverse of that for the North America and Europe region because the annual cycle of temperature for Australia is inversely related to the annual cycle of the insolation anomalies (Fig. 9) and therefore to the direct calendar effects on insolation (Fig. 5). Again, the shapes of the Holocene trends in the adjusted and unadjusted data are noticeably different. In the Australia (Fig. 14b) and North America and Europe (Fig. 14a) examples, relatively large areas are being averaged, and the calendar effect becomes more apparent as the size of the area decreases. Notably, the effect does not completely disappear at the largest scales, i.e., for area-weighted averages for the globe (for ice-free land grid cells) (Fig. 14c). The differences are smaller but still discernible.

In the Northern Hemisphere (African–Asian) monsoon region (0 to 30^∘ N and −30 to 120^∘ E), the calendar effects on precipitation rate are similar to those on temperature in the midlatitudes because the annual cycle of precipitation is roughly in phase with that of insolation (Fig. 7). There is little effect in the winter and spring but a substantial effect in summer and autumn over the interval from 17 to about 3 ka (Fig. 14d). The calendar effect reverses sign between July and August (when the month-length adjusted precipitation rate values are less than the unadjusted ones) and September and October (when the adjusted values are greater than the unadjusted ones). In July, the timing of relative maxima and minima in the two data sets is similar, while in October, in particular, the Holocene precipitation maximum is several thousand years earlier in the adjusted data than in the unadjusted data.

The time-series expression of the latitudinally reversing calendar effect on precipitation rate evident in Fig. 13 (e.g., July vs. October at 127 ka) can be illustrated by comparing precipitation or precipitation minus evaporation (P−E) for the North African (sub-Saharan) monsoon region (5 to 17^∘ N and −5 to 30^∘ E) with the Mediterranean region (31 to 43^∘ N and −5 to 30^∘ E) (Fig. 14e and f). The differences between the adjusted and unadjusted data in the North African region (Fig. 14e) parallel that of the larger monsoon region (Fig. 14d). The Mediterranean region, which is characteristically moister in winter and drier in summer, shows the reverse pattern: when the calendar-adjusted minus unadjusted P−E difference is positive in the monsoon region, it is negative in the Mediterranean region. Dipoles are frequently observed in climatic data, both present-day and paleo, and are usually interpreted in terms of broad-scale circulation changes in the atmosphere or ocean. This example illustrates that they could also be artifacts of the calendar effect. Such changes in the timing of extrema could also influence the interpretation of phase relationships among simulated time series and time series of potential forcing (Joussaume and Braconnot, 1997; Timm et al., 2008; Chen et al., 2011).

There are other interesting patterns in the monthly time series from the transient simulations, some of which are amplified by the calendar effect and others damped. The monthly time series suggest that the traditional meteorological seasons (i.e., December–February, March–May, June–August, September–November) are not necessarily the optimal way to aggregate data – the September time series in Fig. 14 often look like they are more similar to, and should be grouped with, July and August rather than with October and November, the traditional other (northern) autumnal months. Figure 14a (North America and Europe), for example, suggests that the July through November time series are similar in their overall trends, and even more so for the adjusted data (in pink and red). Similarly, months that appear highly correlated over some intervals (e.g., July and June global temperatures from the Last Glacial Maximum to the Holocene) become decoupled at other times. The impacts of the calendar effect on temporal trends in transient simulations (Fig. 14), when compounded by the spatial effects (Figs. 11–13), make it even more likely that spurious climatic mechanisms could be inferred in analyzing transient simulations rather than in the simpler time-slice simulations.

3.5 Summary

Several observations can be made about the calendar effect and its potential role in the interpretation of paleoclimatic simulations and comparisons with observations.

• The variations in eccentricity and perihelion over time are large enough to produce differences in the length of (fixed-angular) months that are as large as 4 or 5 d, as well as differences in the beginning and ending times of months on the order of 10 d or more (Fig. 1).

• These month-length and beginning and ending date differences are large enough to have noticeable impacts on the location in time of a fixed-length month relative to the solstices and hence on the insolation receipt during that interval (Figs. 2 through 5). The average insolation (and its difference from the present) during a fixed-length month will thus include the effects of the orbital variations on insolation and the changing month length.

• However, such insolation effects are not offset by the changing insolation itself but instead can be reinforced or damped (Figs. 7 through 9). (In other words, orbitally related variations in insolation do not take care of the calendar definition issue.)

• The pure calendar effects on temperature and precipitation (illustrated by comparing adjusted and non-adjusted data assuming no climate change; Figs. 11–13) are large, spatially variable, and could easily be mistaken for real paleoclimatic differences (from the present).

• The impact of the calendar effect on transient simulations is also large (Fig. 14), affecting the timing and phasing of maxima and minima, which, when combined with the spatial impacts of the calendar effect, makes transient simulations even more prone to misinterpretation.

4.1 Pseudo-daily interpolation

The first step in adjusting monthly time-step model output to reflect the calendar effect is to interpolate the monthly data (either long-term means or time-series data) to pseudo-daily values (a step that is not required if the data are daily time-step values). It turns out that the most common way of producing pseudo-daily values, linear interpolation between monthly means, is not mean preserving; the monthly (or annual) means of the interpolated daily values will generally not match the original monthly values. An alternative approach, and the one we use here, is the mean-preserving “harmonic” interpolation method of Epstein (1991), which is easy to implement and performs the same function as the parabolic-spline interpolation method of Pollard and Reusch (2002). As is also the case with Pollard and Reusch's (2002) method, Epstein's (1991) approach can occasionally produce overshoots that are physically impossible, as can happen in the application of the method to variables like precipitation, which may have monthly values that alternate between zero and nonzero values. For practical reasons, variables like precipitation are therefore “clamped” at zero, which can introduce small differences between the annual and monthly means of the original and interpolated data, and we illustrate a pathological case of this below.

Figure 15Pseudo-daily interpolated temperature (a, b) and precipitation (c, d) for some representative locations: (a, c) Madison, Wisconsin, USA, (b) Australia, and (d) the Indian Ocean. The original monthly mean data are shown by the black dots and stepped curves (black lines), daily values linearly interpolated between the monthly mean values are shown in blue, and daily values using the mean-preserving approach of Epstein (1991) are shown in red. The annual interpolation error (or the difference between the annual average calculated using the original data and the pseudo-daily interpolated data) is given for the mean-preserving approach in each case. The interpolated data for this figure were generated using the program demo_01_pseudo_daily_interp.f90.

Download

The linear and mean-preserving interpolation methods can be compared using the Climate Forecast System Reanalysis (CFSR) near-surface air temperature and CPC Merged Analysis of Precipitation (CMAP) 1981–2010 long-term mean data (Fig. 15). A typical example for temperature appears in Fig. 15a for a grid point near Madison, Wisconsin (USA). The difference between the annual mean values of the interpolated data for the two approaches is small and similar (ca. $\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{6}}$), but the difference between the original monthly means and the monthly mean of the linearly interpolated daily values can exceed 0.8 ^∘C in some months (e.g., December). (The differences from the original monthly means for the mean-preserving interpolation method are less than $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{3}}$ ^∘C for every month in Fig. 15a.) Figure 15b shows an example for a grid point in Australia, where again the difference between the original monthly means and the monthly means of the linearly interpolated daily values is not negligible (i.e., 0.4 ^∘C). Similar results hold for precipitation (Fig. 15c), whereby the difference can exceed 0.1 mm d⁻¹. Like other harmonic-based approaches, the Epstein (1991) approach can create interpolated curves that are wavy (see Pollard and Reusch, 2002, for a discussion), but these effects are small enough to not be practically important in nearly all cases. The pathological case for precipitation is shown in Fig. 15d at a grid point in the Indian Ocean. Here, the difference between an original monthly mean value and one calculated using the mean-preserving interpolation method reaches −0.12 mm d⁻¹ in March and April, but the differences between the original monthly means and the monthly means of the linearly interpolated daily values are nearly 3 times larger.

The map patterns of the interpolation errors (the monthly mean values recalculated using the linear or mean-preserving pseudo-daily interpolated values minus the original values) appear in Fig. 16. (Note the differing scales for the linear interpolation errors and the mean-preserving interpolation errors.) The linear interpolation errors are quite large, with absolute values exceeding 1 ^∘C and 1 mm d⁻¹, and have distinct seasonal and spatial patterns: underpredictions of Northern Hemisphere temperature in summer (and overpredictions in winter), underpredictions of precipitation in the wet season (e.g., southern Asia in July), and overpredictions in the dry season (southern Asia in May). The magnitude and patterns of these effects again rival those we attempt to infer or interpret in the paleo record. The mean-preserving interpolation errors for temperature are very small and show only vague spatial patterns (note the differing scales). The errors for precipitation are also quite small but can be locally larger, as in the pathological case illustrated above. However, the map patterns of the interpolation errors strongly suggest that those cases are not practically important.

Figure 16Pseudo-daily interpolation errors for CFSR near-surface air temperature (a, c) and CMAP precipitation rate (b, d). (a, b) The interpolation errors, or the differences between the original monthly mean values and the monthly mean values recalculated from linear interpolation of pseudo-daily values. (c, d) The interpolation errors for mean-preserving (Epstein, 1991) interpolation. The errors for linear interpolation of the temperature data (^∘C) range from −1.20851 to 1.29904, with a mean of 0.05664 and standard deviation of 0.16129 (over all months and grid points), while those for mean-preserving interpolation range from −0.00002 to 0.00050, with a mean of −0.0061 and standard deviation of 0.00007. The errors for linear interpolation of the precipitation data (mm d⁻¹) range from −1.10617 to 1.40968, with a mean of 0.00087 and standard deviation of 0.11851, while those for mean-preserving interpolation range from −0.00002 to 0.00383, with a mean of 0.00001 and standard deviation of 0.00163.

The mean-preserving interpolation method is implemented in the Fortran 90 module named pseudo_daily_interp_subs.f90. The subroutine hdaily(…) manages the interpolation, first getting the harmonic coefficients (Eq. 6 of Epstein, 1991) using the subroutine named harmonic_coeffs(…) and then applying these coefficients in the subroutine xdhat(…) to get the interpolated values.

4.2 Month-length calculations

The calculation of the length and the beginning, middle, and ending (real number or fractional) days of each month at a particular time is based on an approach for calculating orbital position as a function of time using Kepler's equation:

where M is the angular position along a circular orbit (referred to by astronomers as the mean anomaly), ε is eccentricity, and E is the eccentric anomaly (Curtis, 2014; Eq. 3.14). Given the angular position of the orbiting body (Earth) along the elliptical orbit, θ (the “true anomaly”), E can be found using the following expression (Curtis, 2014; Eq. 3.13b):

Substituting E into Eq. (1) gives us M, and then the time since perihelion is given by

where T is the orbital period (i.e., the length of the year) (Curtis, 2014; Eq. 3.15).

This expression can be used to determine the traverse time or time of flight of individual days or of segments of the orbit equivalent to the fixed-angular definition of months or seasons. Doing so involves determining the traverse times between the vernal equinox and perihelion, between the vernal equinox and 1 January (set at the appropriate number of degrees prior to the vernal equinox for a particular calendar), and the angle between perihelion and 1 January and using these values to translate time since perihelion to time since 1 January. The true anomaly angles along the elliptical orbit (θ) are determined using the present-day (e.g., 1950 CE) definitions of the months in different calendars (e.g., January is defined as having 30, 31, and 31 d in calendars with a 360, 365, or 366 d year, respectively). For example, January in a 365 d year is defined as the arc or “month angle” between 0.0 and 31.0 d × (360.0^∘ ∕ 365.0 d). Note that when perihelion is in the Northern Hemisphere winter, the arc may begin after 1 January as a consequence of the occurrence of shorter winter months, and when perihelion is in the Northern Hemisphere summer, the arc may begin before 1 January as a consequence of longer winter months (Fig. 1).

We also implemented the approximation approach described by Kutzbach and Gallimore (1988; Appendix A) for calculating month lengths. There were no practical differences between their approach and our implementation of Kepler's equation based on the Curtis (2014) approach.

Application of this algorithm requires as input eccentricity and the longitude of perihelion (^∘) relative to the vernal equinox, and the generalization of the approach to other calendars, such as the proleptic Gregorian calendar (that includes leap years; http://cfconventions.org), also requires the (real number or fractional) day of the vernal equinox. To calculate the orbital parameters using the Berger (1978) solution and the timing of the (northern) vernal equinox (as well as insolation itself), we adapted a set of programs provided by the National Aeronautics and Space Administration (NASA) Goddard Institute for Space Studies (GISS) (now available at https://web.archive.org/web/20150920211936/http://data.giss.nasa.gov/ar5/solar.html, last access: 11 August 2019).

4.3 Simulation ages and simulation years

Inspection shows that different climate models employ different starting dates in their output files for both present-day (piControl) and paleo (e.g., midHolocene) simulations (https://esgf-node.llnl.gov/projects/cmip5/). For models that use a noleap (constant 365 d year) calendar, such as CCSM4 (Otto-Bliesner, 2014), the starting date is not an issue, but for MPI-ESM-P (Jungclaus et al., 2012), which uses a proleptic Gregorian calendar, or CNRM-CM5 (Sénési et al., 2014), which uses a “standard” (i.e., mixed Julian–Gregorian) calendar, as examples, the specific starting date influences the date of the vernal equinox through the occurrence of individual leap years. For example, in the CMIP5–PMIP4 midHolocene simulations, output from MPI-ESM-P starts in 1850 CE and that from CNRM-CM5 in 2050 CE (and it can be verified that leap years in those output files occur in a fashion consistent with the “modern” calendar). Consequently, we need to make a distinction between two notions of time here: (1) the simulation age, expressed in (negative) years BP (i.e., before 1950 CE), and (2) the simulation year, expressed in years CE. The simulation age controls the orbital parameter values, while the simulation year, along with the specification of the CF-compliant calendar attribute (http://cfconventions.org), controls the date and time of the vernal equinox.

4.4 Month-length programs and subprograms

Month lengths are calculated in the subroutine get_month_lengths(…) (contained in the Fortran 90 module named month_length_subs.f90), which in turn calls the subroutine monlen(…) to get real-type month lengths for a particular simulation age and year. (The subroutine get_month_lengths(…) can be exercised to produce tables of month lengths, beginning, middle, and ending days of the kind used to produce Figs. 1–5 and 7–9 using a driver program named month_length.f90.) The subroutine get_month_lengths(…) uses two other modules, GISS_orbpar_subs.f90 and GISS_srevents_subs.f90 (based on programs originally downloaded from GISS; now available at https://web.archive.org/web/20150920211936/http://data.giss.nasa.gov/ar5/solar.html), to get the orbital parameters and vernal equinox dates.

The specific tasks involved in the calculation of either a single year's set of month lengths, or a series of month lengths for multiple years, include the following steps, implemented in get_month_lengths(…):

• 1.

  generate a set of “target” dates based on the simulation ages and simulation years;

• 2.

  obtain the orbital parameters for 0 ka (1950 CE), which will be used to adjust the calculated month-length values to the conventional definition of months for 1950 CE as the reference year; and

• 3.

  obtain the present-day (i.e., 1950 CE) month lengths (along with the beginning, middle, and ending days relative to 1 January) for the appropriate calendar using the subroutine monlen(…).

Then loop over the simulation ages and simulation years, and for each combination do the following:

• 4.

  obtain the orbital parameters for each simulation age using the subroutine GISS_orbpars(…);

• 5.

  calculate real-type month lengths (along with the beginning, middle, and ending days relative to 1 January) for the appropriate calendar using monlen(…);

• 6.

  adjust (using the subroutine adjust_to_ref_length(…)) those month-length values to the reference year (e.g., 1950 CE) and its conventional set of month-length definitions so that, for example, January will have 31 d, February 28 or 29 d, etc., in that reference year;

• 7.

  further adjust the month-length values to ensure that the individual monthly values will sum exactly to the year length in days using the subroutine adjust_to_yeartot(…);

• 8.

  convert real-type month lengths to integers using the subroutine integer_monlen(…) (these integer values are not used anywhere but may be useful in conceptualizing the pattern of month-length variations over time);

• 9.

  get integer-valued beginning, middle, and ending days for each month; and

• 10.

  determine the mid-March day using the subroutine GISS_srevents(…) to get the vernal equinox date for calendars in which it varies.

4.5 Month-length tables and time series

Tables and time series of month lengths, beginning, middle, and ending days, and dates of the vernal equinox can be calculated using the program month_length.f90. This program reads an “info file” (month_length_info.csv) consisting of an identifying output file name prefix, the calendar type, the beginning and ending simulation age (in years BP), the age step, the beginning simulation year (in years CE), and the number of simulation years. Note that in the approach described above, orbital parameters are calculated once per year (step 4 in Sect. 4.4) and are assumed to apply for the whole year. This assumption can lead to small differences (ranging from −0.000863 to 0.000787 d over the past 22 kyr with a mean of −0.00000389 d) in the ending day of one year and the beginning day of the next.

5.1 Interpolation and (re)aggregation

The pseudo-daily interpolation and (re)aggregation are done using two subroutines, mon_to_day_ts(…) and day_to_mon_ts(…), in the module calendar_effects_subs.f90. The pseudo-daily interpolation is done a year at a time, creating slight discontinuities between one year and the next in the case of transient or multiyear “snapshot” simulations. The subroutine mon_to_day_ts(…) has options for smoothing those discontinuities and restoring the long-term mean of the interpolated daily data to that of the original monthly data.

The (re)aggregation of the daily data is also done a year at a time by collecting the daily data for a particular year and “padding” it at the beginning and end with data from the previous and following year if available, as in transient or multiyear simulations (to accommodate the fact that under some orbital configurations the first day of the current year may occur in the previous year or the last day in the following year; Fig. 1). For example, at 6 ka, the changes in the shape of the orbit and the consequently longer months from January through March (32.5, 29.5, and 32.4 d, respectively) displace the beginning of January 4 d into the previous year, with the last day of December consequently falling just before day 361 in a 365 d year. In the case of long-term mean “climatological” data (Aclim data; see Sect. 5.2), the padding is done with the ending and beginning days of the single year of pseudo-daily data.

The calculation of monthly means is done by calculating weighted averages of the days that overlap a particular month as defined by the (real number or fractional) beginning and ending days of that month (from the subroutine get_month_lengths(…)). Each whole day in that interval gets a weight of 1.0, and each partial day gets a weight proportional to its part of a whole day. It should be noted that in transient simulations, annual averages, constructed either by averaging actual or pseudo-daily data (or by month-length weighted averages), will differ from the unadjusted data.

5.2 Processing individual netCDF files

The cal_adjust_PMIP.f90 program reads an info file that provides file and variable details and can handle CMIP6–PMIP4-formatted files (https://pcmdi.llnl.gov/CMIP6/Guide/modelers.html#5-model-output-requirements, last access: 11 August 2019) as they become available. The fields in the info file include (for each netCDF file) the activity (PMIP3 or PMIP4), the variable (e.g., tas, pr), the realm-plus-time-frequency type (e.g., Amon, Aclim), the model name, the experiment name (e.g., midHolocene), the ensemble member (e.g., r1i1p1), the grid label (for PMIP4 files), and the simulation year beginning date and ending date (as a YYYYMM or YYYYMMDD string). An input file name suffix field is also read (which is usually blank but is “-clim” for Aclim-type files), as is an output file name suffix field (e.g., _cal_adj), which is added to the output file name to indicate that it has been modified from the original. The info file also contains the simulation age beginning and end (in years BP), the increment between simulation ages (usually 1 in the application here), the beginning simulation year (years CE), the number of simulation years, and the paths to the source and adjusted files. This information could also be obtained by parsing the netCDF file names and reading the calendar attribute and time-coordinate variables, but that would add to the complexity of the program.

The output netCDF files have the string _cal_adj appended to the end of the file name. In the case of monthly time series (e.g., Amon) or long-term means (e.g., Aclim) the file names are otherwise the same as the input data. In the case of the daily input data, with “day” as the realm-plus-time-frequency string, that string is changed to Amon2.

The adjustment of a file using cal_adjust_PMIP.f90 includes the following steps:

1. read the info file, construct various file names, and allocate month-length variables;

2. generate month lengths using the subroutine get_month_lengths(…); and

3. open input and output netCDF files. For each file,

4. redefine the time-coordinate variable as appropriate using the subroutines new_time_day(…) and new_time_month(…) in the module CMIP_netCDF_subs.f90;

5. create the new netCDF file, copy the dimension and global attributes from the input file using the subroutine copy_dims_and_glatts(…), and define the output variable using the subroutine define_outvar(…);

6. get the input variable to be adjusted;

7. for each model grid point, get calendar-adjusted values as described above using the subroutines mon_to_day_ts(…) and day_to_mon_ts(…); and

8. write out the adjusted data, and close the output file.

5.3 Further examples

Five other main programs that serve as “drivers” for some of the subroutines or that demonstrate particular aspects of procedures used here are included in the GitHub repository for the programs (https://github.com/pjbartlein/PaleoCalAdjust, last access: 11 August 2019).

• GISS_orbpar_driver.f90 and GISS_srevents_driver.f90 are the main programs that call the subroutines GISS_orbpars(…) and GISS_srevents(…) to produce tables of orbital parameters and “solar events” like the dates of equinoxes, solstices, and perihelion and aphelion.

• demo_01_pseudo_daily_interp.f90 is the main program that demonstrates linear and mean-preserving pseudo-daily interpolation.

• demo_02_adjust_1yr.f90 is the main program that demonstrates the paleo calendar adjustment of a single year's data.

• demo_03_adjust_TraCE_ts.f90 is the main program that demonstrates the adjustment of a 22 040-year-long time series of monthly TraCE-21ka data.