Refactored air temperature and melt functions

[andrewdnolan]

Reference Code:

Mean annual air temperature was formerly a function of day of year ($d$), such that:

T_{\rm s}(d,z) = T_{\rm ma} + \alpha \cos\left( \frac{2 \pi \left(d - T_{\rm p} \right)}{365} \right) + \left(z - z_{\rm ref}\right) \frac{\partial T}{\partial z}, 

In order to convert time in year to day of year, the fractional part of the timestep is parsed and converted to an integer [1,365]. For example if t_i=10.1 then doy_i = int(0.1*365)≈37. As you can imagine there is some (small) temporal resolution lost when converting from the float to an integer on the interval [1,365]. ( In the above example 36.5 != 37). The error is consistent when a constant timestep is used, and doesn't seem to have any effect on the solution. But, when variable timesteps are used there is a bias in the mean annual air temperature (Issue: #13).

Previously, the doy of start (t_i) and end (t_ip1) of the current timestep were calculated and the used to index a vector (of length 365) containing the mean annual air temperature curve.

Changes to Code:

I've re-parametrized the air temperature and melt function to accept time in years ($t$), such that:

T_{\rm s}(t,z) = T_{\rm ma} + \alpha \cos\left( 2 \pi \left( t -T_{\rm p}/365 \right) \right) + \left(z - z_{\rm ref}\right) \frac{\partial T}{\partial z}, 

the units of all the input parameters (e.g. $\alpha, T_{\rm p}, \ldots$) have remain the same, and the code will now convert units where needed to keep things consistent.

The new code evaluates the temperature at 100 evenly spaced points between the start (t_i) and end (t_ip1) of the current timestep. In the case of timesteps less than year, the usual for our thermomechanically coupled experiments, this will serve to increase the resolution of the mean annual airtemp calculation.

[andrewdnolan]

I've rerun the transient initialization for the reference surface forcing ($\hat T_{\rm ma} = -8.5, \Delta \dot{b} = -0.37$) with the new parameterization (orange line). The old implementation is shown in blue. While there is some difference in time, the resulting configuration close enough that I'm comfortable going forward as is.