Skip to content

Surface Energy Budget

gamma_linear_approx(time, temp, lambda_const, lambda_time_lag=None, lambda_nl=None, temp_anom_nl=None)

OUTDATED FUNCTION - USED FOR get_temp_fourier but now use get_temp_fourier_numerical

This approximates \(\Gamma^{\uparrow} = LW^{\uparrow} - LW^{\downarrow} + LH^{\uparrow} + SH^{\uparrow}\) as:

\[\Gamma^{\uparrow} \approx \lambda_0 + \sum_{i=1}^{N_{\lambda}}\lambda_i T(t-\Lambda_i) + \sum_{j=2}^{j_{max}}\lambda_{nl_j}\left(T'^{j}(t) - \overline{T'^j}\right)\]

where:

  • \(T' = T(t) - \overline{T}\) is the surface temperature anomaly
  • \(LW^{\uparrow}\) is upward longwave radiation at the surface i.e. \(\sigma T^4\) or lwup_sfc output by Isca. Units: \(Wm^{-2}\).
  • \(LW^{\downarrow}\) is downward longwave radiation at the surface i.e. lwdn_sfc output by Isca. Units: \(Wm^{-2}\).
  • \(LH^{\uparrow}\) is the upward latent heat radiation at the surface i.e. flux_lhe output by Isca. Units: \(Wm^{-2}\).
  • \(SH^{\uparrow}\) is the upward sensible heat radiation at the surface i.e. flux_t output by Isca. Units: \(Wm^{-2}\).

Parameters:

Name Type Description Default
time ndarray

float [n_time]
Time in days (assumes periodic e.g. annual mean, so time = np.arange(360)).

 required
temp ndarray

float [n_time]
Surface temperature, \(T\) for each day in time. Assumes periodic so temp[0] is the temperature at time time[-1]+1.

 required
lambda_const ndarray

float [n_lambda+1]
The constants \(\lambda_i\) used in the approximation.
lambda_const[0] is \(\lambda_0\) and lambda_const[i] is \(\lambda_{i}\) for \(i>0\).

 required
lambda_time_lag Optional[ndarray]

float [n_lambda]
The constants \(\Lambda_i\) used in the approximation.
lambda_time_lag[0] is \(\Lambda_1\) and lambda_time_lag[i] is \(\Lambda_{i+1}\) for \(i>0\).

 None
lambda_nl Optional[Union[float, ndarray]]

float [n_lambda_nl] The constants \(\lambda_{nl}\) used in the approximation. [0] is squared contribution, [1] is cubed, ...

 None
temp_anom_nl Optional[ndarray]

float [n_time]
The value of \(T'(t)\) to use in the calculation of non-linear part: \(\sum_{j=2}^{j_{max}}\lambda_{nl_j}\left(T'^{j}(t) - \overline{T'^j}\right)\). If None, will compute from temp.

 None

Returns:

Type Description
ndarray

float [n_time]
The approximation \(\Gamma^{\uparrow} \approx \lambda_0 + \sum_{i=1}^{N_{\lambda}}\lambda_i T(t-\Lambda_i) + \lambda_{sq}\left(T(t) - \overline{T}\right)^2\) with units of \(Wm^{-2}\).
Source code in isca_tools/thesis/surface_energy_budget.py 
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
def gamma_linear_approx(time: np.ndarray, temp: np.ndarray,
                        lambda_const: np.ndarray, lambda_time_lag: Optional[np.ndarray] = None,
                        lambda_nl: Optional[Union[float, np.ndarray]] = None,
                        temp_anom_nl: Optional[np.ndarray] = None) -> np.ndarray:
    """
    OUTDATED FUNCTION - USED FOR `get_temp_fourier` but now use `get_temp_fourier_numerical`

    This approximates $\Gamma^{\\uparrow} = LW^{\\uparrow} - LW^{\\downarrow} + LH^{\\uparrow} + SH^{\\uparrow}$ as:

    $$\Gamma^{\\uparrow} \\approx \lambda_0 + \sum_{i=1}^{N_{\lambda}}\lambda_i T(t-\Lambda_i) +
    \sum_{j=2}^{j_{max}}\lambda_{nl_j}\\left(T'^{j}(t) - \overline{T'^j}\\right)$$

    where:

    * $T' = T(t) - \overline{T}$ is the surface temperature anomaly
    * $LW^{\\uparrow}$ is upward longwave radiation at the surface i.e. $\\sigma T^4$ or `lwup_sfc` output by Isca.
        Units: $Wm^{-2}$.
    * $LW^{\\downarrow}$ is downward longwave radiation at the surface i.e. `lwdn_sfc` output by Isca.
        Units: $Wm^{-2}$.
    * $LH^{\\uparrow}$ is the upward latent heat radiation at the surface i.e. `flux_lhe` output by Isca.
        Units: $Wm^{-2}$.
    * $SH^{\\uparrow}$ is the upward sensible heat radiation at the surface i.e. `flux_t` output by Isca.
        Units: $Wm^{-2}$.

    Args:
        time: `float [n_time]`</br>
            Time in days (assumes periodic e.g. annual mean, so `time = np.arange(360)`).
        temp: `float [n_time]`</br>
            Surface temperature, $T$ for each day in `time`. Assumes periodic so `temp[0]` is the temperature
            at time `time[-1]+1`.
        lambda_const: `float [n_lambda+1]`</br>
            The constants $\lambda_i$ used in the approximation.</br>
            `lambda_const[0]` is $\lambda_0$ and `lambda_const[i]` is $\lambda_{i}$ for $i>0$.
        lambda_time_lag: `float [n_lambda]`</br>
            The constants $\Lambda_i$ used in the approximation.</br>
            `lambda_time_lag[0]` is $\Lambda_1$ and `lambda_time_lag[i]` is $\Lambda_{i+1}$ for $i>0$.
        lambda_nl: `float [n_lambda_nl]`
            The constants $\lambda_{nl}$ used in the approximation. `[0]` is squared contribution, `[1]` is cubed, ...
        temp_anom_nl: `float [n_time]`</br>
            The value of $T'(t)$ to use in the calculation of non-linear part:
            $\sum_{j=2}^{j_{max}}\lambda_{nl_j}\\left(T'^{j}(t) - \overline{T'^j}\\right)$.
            If `None`, will compute from `temp`.

    Returns:
        `float [n_time]`</br>
            The approximation $\Gamma^{\\uparrow} \\approx \lambda_0 + \sum_{i=1}^{N_{\lambda}}\lambda_i T(t-\Lambda_i)
            + \lambda_{sq}\\left(T(t) - \overline{T}\\right)^2$ with units of $Wm^{-2}$.
    """
    if lambda_nl is not None:
        #  deal with case when float given as lambda_nl
        if not hasattr(lambda_nl, "__len__"):
            lambda_nl = np.asarray([lambda_nl])
    n_lambda = len(lambda_const) - 1
    if lambda_time_lag is None:
        lambda_temp = np.sum(lambda_const[1:]) * temp
    else:
        # To apply phase shift, need spline so can compute temp anomaly at times outside range in `time`.
        temp_spline_fit = CubicSpline(np.append(time, time[-1] + 1), np.append(temp, temp[0]),
                                      bc_type='periodic')
        lambda_temp = np.zeros_like(temp)
        for i in range(n_lambda):
            lambda_temp += lambda_const[1 + i] * temp_spline_fit(time - lambda_time_lag[i])
    gamma_nl = np.zeros_like(temp)
    if lambda_nl is not None:
        if temp_anom_nl is None:
            temp_anom_nl = temp - np.mean(temp)
        n_lambda_nl = len(lambda_nl)
        for j in range(n_lambda_nl):
            gamma_nl += lambda_nl[j] * (temp_anom_nl ** (j + 2) - np.mean(temp_anom_nl ** (j + 2)))
    return lambda_const[0] + lambda_temp + gamma_nl

get_param_dimensionless(var, heat_capacity=None, n_year_days=None, sw_fourier_amp1=None, lambda_const=None, day_seconds=86400, invert=False)

Returns dimensionless versions of empirical fitting parameter, without changing the sign. There are three possibilities for how var is made dimensionless:

  • If provide heat_capacity, \(C\), and n_year_days, \(\mathcal{T}\), then will assume var is \(\lambda\) or \(\lambda_{phase}\). Will return \(\lambda' = \frac{\lambda}{2\pi fC}\) where \(f=1/\mathcal{T}\).
  • If provide sw_fourier_amp1, \(F_1\) and lambda_const, \(\lambda\), then will assume var is \(\lambda_{sq}\). Will return \(\lambda_{sq}' = \frac{\lambda_{sq}F_1}{2\lambda^2}\).
  • If provide sw_fourier_amp1, \(F_2\), only then will assume var is \(\Lambda_{cos}\) or \(\Lambda_{sin}\): \(\Lambda_{cos}' = \Lambda_{cos}/F_1\).

Parameters:

Name Type Description Default
var Union[float, ndarray, DataArray]

Parameter to make dimensionless.

 required
heat_capacity Optional[Union[float, ndarray, DataArray]]

\(C\), the heat capacity of the surface in units of \(JK^{-1}m^{-2}\).
Obtained from mixed layer depth of ocean using get_heat_capacity.

 None
n_year_days Optional[int]

Number of days in the year.

 None
sw_fourier_amp1 Optional[Union[float, ndarray, DataArray]]

The first harmonic amplitude Fourier coefficient for shortwave radiation, \(F_1\).

 None
lambda_const Optional[Union[float, ndarray, DataArray]]

The linear constant \(\lambda\) used in the approximation for \(\Gamma^{\uparrow}. In normal form with units of Wm\)^{-2}\(K\)^{-1}$.

 None
day_seconds float

Duration of a day in seconds.

 86400
invert bool

If True, will return the parameter with dimensions given the dimensionless version.

 False

Returns:

Name Type Description
var_dim Union[float, ndarray, DataArray]

Dimensionless version of var.
Source code in isca_tools/thesis/surface_energy_budget.py 
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
def get_param_dimensionless(var: Union[float, np.ndarray, xr.DataArray],
                            heat_capacity: Optional[Union[float, np.ndarray, xr.DataArray]] = None,
                            n_year_days: Optional[int] = None,
                            sw_fourier_amp1: Optional[Union[float, np.ndarray, xr.DataArray]] = None,
                            lambda_const: Optional[Union[float, np.ndarray, xr.DataArray]] = None,
                            day_seconds: float = 86400, invert: bool = False
                            ) -> Union[float, np.ndarray, xr.DataArray]:
    """
    Returns dimensionless versions of empirical fitting parameter, without changing the sign.
    There are three possibilities for how `var` is made dimensionless:

    * If provide `heat_capacity`, $C$, and `n_year_days`, $\mathcal{T}$, then will assume `var` is $\lambda$ or
    $\lambda_{phase}$. Will return $\lambda' = \\frac{\lambda}{2\pi fC}$ where $f=1/\mathcal{T}$.
    * If provide `sw_fourier_amp1`, $F_1$ and `lambda_const`, $\lambda$, then will
    assume `var` is $\lambda_{sq}$. Will return $\lambda_{sq}' = \\frac{\lambda_{sq}F_1}{2\lambda^2}$.
    * If provide `sw_fourier_amp1`, $F_2$, only then will assume `var` is $\Lambda_{cos}$ or $\Lambda_{sin}$:
    $\Lambda_{cos}' = \Lambda_{cos}/F_1$.

    Args:
        var: Parameter to make dimensionless.
        heat_capacity: $C$, the heat capacity of the surface in units of $JK^{-1}m^{-2}$.</br>
            Obtained from mixed layer depth of ocean using
            [`get_heat_capacity`](../utils/radiation.md#isca_tools.utils.radiation.get_heat_capacity).
        n_year_days: Number of days in the year.
        sw_fourier_amp1: The first harmonic amplitude Fourier coefficient for shortwave radiation, $F_1$.
        lambda_const: The linear constant $\lambda$ used in the approximation for
            $\Gamma^{\\uparrow}. In normal form with units of Wm$^{-2}$K$^{-1}$.
        day_seconds: Duration of a day in seconds.
        invert: If `True`, will return the parameter with dimensions given the dimensionless version.

    Returns:
        var_dim: Dimensionless version of `var`.
    """
    # May be issue here with var changing sign, but don't think so as only negative is sw_fourier_amp1 and square that
    # But in southern hemisphere, may be different with sw_fourier_amp2
    if (heat_capacity is not None) and (n_year_days is not None):
        # lambda_phase or lambda_const
        f = 1 / (n_year_days * day_seconds)
        if not invert:
            var = var / (2 * np.pi * f * heat_capacity)
        else:
            var = var * (2 * np.pi * f * heat_capacity)
    elif (sw_fourier_amp1 is not None) and (lambda_const is not None):
        # lambda_sq
        if not invert:
            var = var / (2 * lambda_const ** 2) * sw_fourier_amp1
        else:
            var = var * (2 * lambda_const ** 2) / sw_fourier_amp1
    elif lambda_const is None:
        # lambda_cos and lambda_sin
        if not invert:
            var = var / sw_fourier_amp1
        else:
            var = var * sw_fourier_amp1
    else:
        raise ValueError("Combination of parameters provided not correct for any parameter")
    return var

get_temp_extrema_analytic(sw_fourier_amp1, heat_capacity, lambda_const, lambda_phase=0, lambda_sq=0, sw_fourier_amp2=0, n_year_days=360, day_seconds=86400)

This will return the analytic solution for the times and amplitudes of the extrema of the fourier solution of \(T'\) in the following form of the surface energy budget:

\[ C\frac{\partial T'}{\partial t} = F(t) - \lambda_0 - \lambda T'(t) - \lambda_{phase}T'(t-\mathcal{T}/4) - \lambda_{sq} T'^2(t) \]

Parameters:

Name Type Description Default
sw_fourier_amp1 Union[float, ndarray]

float [n_regions]
The first harmonic amplitude fourier coefficients for shortwave radiation, \(SW^{\downarrow}\): \(F_1\).

 required
heat_capacity float

\(C\), the heat capacity of the surface in units of \(JK^{-1}m^{-2}\).
Obtained from mixed layer depth of ocean using get_heat_capacity.

 required
lambda_const Union[float, ndarray]

The constant \(\lambda\) used in the approximation for \(\Gamma^{\uparrow} = LW^{\uparrow} - LW^{\downarrow} + LH^{\uparrow} + SH^{\uparrow}\).

 required
lambda_phase Union[float, ndarray]

The constants \(\lambda_{phase}\) used in the approximation for \(\Gamma^{\uparrow}\).

 0
lambda_sq Union[float, ndarray]

The constant \(\lambda_{sq}\) used in the approximation for \(\Gamma^{\uparrow}\).

 0
sw_fourier_amp2 Union[float, ndarray]

float [n_regions]
The second harmonic amplitude fourier coefficients for shortwave radiation, \(SW^{\downarrow}\): \(F_1\).

 0
n_year_days int

Number of days in a year.

 360
day_seconds float

Duration of a day in seconds.

 86400

Returns:

Name Type Description
time_extrema1 Union[float, ndarray]

Time in days of extrema to occur first
time_extrema2 Union[float, ndarray]

Time in days of extrema to occur last
amp_extrema1 Union[float, ndarray]

Absolute amplitude of extrema to occur first
amp_extrema2 Union[float, ndarray]

Absolute amplitude of extrema to occur last
Source code in isca_tools/thesis/surface_energy_budget.py 
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
def get_temp_extrema_analytic(sw_fourier_amp1: Union[float, np.ndarray], heat_capacity: float,
                              lambda_const: Union[float, np.ndarray], lambda_phase: Union[float, np.ndarray] = 0,
                              lambda_sq: Union[float, np.ndarray] = 0, sw_fourier_amp2: Union[float, np.ndarray] = 0,
                              n_year_days: int = 360, day_seconds: float = 86400
                              ) -> Tuple[Union[float, np.ndarray], Union[float, np.ndarray], Union[float, np.ndarray],
Union[float, np.ndarray]]:
    """
    This will return the analytic solution for the times and amplitudes of the extrema of the fourier
    solution of $T'$ in the following form of the surface energy budget:

    $$
    C\\frac{\partial T'}{\partial t} = F(t) - \lambda_0 - \lambda T'(t) - \lambda_{phase}T'(t-\mathcal{T}/4) -
    \lambda_{sq} T'^2(t)
    $$

    Args:
        sw_fourier_amp1: `float [n_regions]`</br>
            The first harmonic amplitude fourier coefficients for shortwave radiation, $SW^{\\downarrow}$: $F_1$.
        heat_capacity: $C$, the heat capacity of the surface in units of $JK^{-1}m^{-2}$.</br>
            Obtained from mixed layer depth of ocean using
            [`get_heat_capacity`](../utils/radiation.md#isca_tools.utils.radiation.get_heat_capacity).
        lambda_const: The constant $\lambda$ used in the approximation for
            $\Gamma^{\\uparrow} = LW^{\\uparrow} - LW^{\\downarrow} + LH^{\\uparrow} + SH^{\\uparrow}$.</br>
        lambda_phase: The constants $\lambda_{phase}$ used in the approximation for $\Gamma^{\\uparrow}$.
        lambda_sq: The constant $\lambda_{sq}$ used in the approximation for $\Gamma^{\\uparrow}$.
        sw_fourier_amp2: `float [n_regions]`</br>
            The second harmonic amplitude fourier coefficients for shortwave radiation, $SW^{\\downarrow}$: $F_1$.
        n_year_days: Number of days in a year.
        day_seconds: Duration of a day in seconds.

    Returns:
        time_extrema1: Time in days of extrema to occur first
        time_extrema2: Time in days of extrema to occur last
        amp_extrema1: Absolute amplitude of extrema to occur first
        amp_extrema2: Absolute amplitude of extrema to occur last
    """
    f = 1 / (n_year_days * day_seconds)
    if sw_fourier_amp2 == 0:
        if lambda_sq != 0:
            raise ValueError('Cannot solve for non-zero lambda_sq with single harmonic - '
                             'get extrema numerically instead')
        tan_phase = (2 * np.pi * heat_capacity * f - lambda_phase) / lambda_const
        time_extrema1 = np.arctan(tan_phase) / (2 * np.pi) * n_year_days
        time_extrema2 = time_extrema1 + n_year_days / 2
        amp_extrema1 = np.abs(sw_fourier_amp1 / lambda_const) / (np.sqrt(1 + tan_phase ** 2))
        amp_extrema2 = amp_extrema1
    return time_extrema1, time_extrema2, amp_extrema1, amp_extrema2

get_temp_extrema_numerical(time, temp, smooth_window=1, smooth_method='convolve', order='time')

Given the temperature temp, this will return the times and amplitudes of the maxima and minima. The extrema will be returned in time order i.e. if minima occurs first, it will be returned first.

Parameters:

Name Type Description Default
time ndarray

float [n_time]
Time in days (assumes periodic e.g. annual mean, so time = np.arange(360))

 required
temp ndarray

float [n_time]
Value of temperature at each time. Again, assume periodic.

 required
smooth_window int

Number of time steps to use to smooth temp before finding extrema. Smaller equals more accurate fit. 1 is perfect fit.

 1
smooth_method str

convolve or spline
If convolve, will smooth via convolution with window of length smooth_window. If spline, will fit a spline using every smooth_window values of time.

 'convolve'
order Literal['time', 'min', 'max']

If time, will return extrema in time order. If min, will return min first. If max will return max first.

 'time'

Returns:

Name Type Description
time_extrema1 float

Time in days of extrema to occur first (if order='time')
time_extrema2 float

Time in days of extrema to occur last (if order='time')
amp_extrema1 float

Absolute amplitude of extrema to occur first (if order='time')
amp_extrema2 float

Absolute amplitude of extrema to occur last (if order='time')
Source code in isca_tools/thesis/surface_energy_budget.py 
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
def get_temp_extrema_numerical(time: np.ndarray, temp: np.ndarray, smooth_window: int = 1,
                               smooth_method: str = 'convolve',
                               order: Literal['time', 'min', 'max'] = 'time') -> Tuple[float, float, float, float]:
    """
    Given the temperature `temp`, this will return the times and amplitudes of the maxima and minima. The extrema
    will be returned in time order i.e. if minima occurs first, it will be returned first.

    Args:
        time: `float [n_time]`</br>
            Time in days (assumes periodic e.g. annual mean, so `time = np.arange(360)`)
        temp: `float [n_time]`</br>
            Value of temperature at each time. Again, assume periodic.
        smooth_window: Number of time steps to use to smooth `temp` before finding extrema.
            Smaller equals more accurate fit. `1` is perfect fit.
        smooth_method: `convolve` or `spline`</br>
            If `convolve`, will smooth via convolution with window of length `smooth_window`.
            If `spline`, will fit a spline using every `smooth_window` values of `time`.
        order: If `time`, will return extrema in time order. If `min`, will return `min` first.
            If `max` will return `max` first.

    Returns:
        time_extrema1: Time in days of extrema to occur first (if `order='time'`)
        time_extrema2: Time in days of extrema to occur last (if `order='time'`)
        amp_extrema1: Absolute amplitude of extrema to occur first (if `order='time'`)
        amp_extrema2: Absolute amplitude of extrema to occur last (if `order='time'`)
    """
    time_extrema = {}
    amp_extrema = {}
    for key in ['min', 'max']:
        var_use, spline_use = numerical.get_var_extrema_date(time, temp - np.mean(temp),
                                                             smooth_window=smooth_window, type=key, max_extrema=1,
                                                             smooth_method=smooth_method)
        time_extrema[key] = var_use[0]
        amp_extrema[key] = np.abs(spline_use(time_extrema[key]))
    # Put output in time order
    if order == 'min':
        return time_extrema['min'], time_extrema['max'], amp_extrema['min'], amp_extrema['max']
    elif order == 'max':
        return time_extrema['max'], time_extrema['min'], amp_extrema['max'], amp_extrema['min']
    elif order == 'time':
        if time_extrema['min'] <= time_extrema['max']:
            return time_extrema['min'], time_extrema['max'], amp_extrema['min'], amp_extrema['max']
        else:
            return time_extrema['max'], time_extrema['min'], amp_extrema['max'], amp_extrema['min']

get_temp_extrema_theory(heat_capacity, sw_amp1, sw_amp2, lambda_const, lambda_phase=0, lambda_sq=0, lambda_cos=0, lambda_sin=0, extrema_ind=1, n_year_days=360, day_seconds=86400, numerical=False, approx_level=None, coef_nl_phase=False)

Find the phase shift and amplitude of the extrema of the analytical temperature solution relative to the first-harmonic extremum, i.e. solve for points where \(\partial T/\partial\Delta = 0\).

This function uses the shifted representation \(T_s(\Delta) - \overline{T}_s = \pm a_1 \sqrt{1-y^2} + a_2 (1-2y^2) + a_3 y \sqrt{1-y^2}\) with \(y = \sin(2\pi f \Delta)\), and decomposes the contribution of each empirical parameter to the timing and amplitude of the extrema. Depending on the options, it either uses analytic coefficients or computes them numerically by solving \(\partial T/\partial\Delta = 0\).

Parameters:

Name Type Description Default
heat_capacity float

\(C\), the surface heat capacity in units of \(JK^{-1}m^{-2}\).

 required
sw_amp1 float

\(F_1\), the amplitude of the first harmonic of downward shortwave radiation at the surface (the \(n=1\) Fourier coefficient of \(SW^{\downarrow}\)), in units of \(Wm^{-2}\).

 required
sw_amp2 float

\(F_2\), the amplitude of the second harmonic of downward shortwave radiation at the surface (the \(n=2\) Fourier coefficient of \(SW^{\downarrow}\)), in units of \(Wm^{-2}\).

 required
lambda_const float

The linear constant \(\lambda\) used in the approximation for \(\Gamma^{\uparrow} = LW^{\uparrow} - LW^{\downarrow} + LH^{\uparrow} + SH^{\uparrow}\).

 required
lambda_phase float

The constant \(\lambda_{phase}\) multiplying the phase-lag term \(\tfrac{1}{2}\lambda_{phase}(T(t-\mathcal{T}/4) - T(t+\mathcal{T}/4))\) in \(\Gamma^{\uparrow}\).

 0
lambda_sq float

The constant \(\lambda_{sq}\) multiplying the quadratic term \(-\lambda_{sq}T'^2\) in \(\Gamma^{\uparrow}\).

 0
lambda_cos float

The constant \(\Lambda_{cos}\) multiplying the \(\cos(4\pi t/\mathcal{T})\) term in \(\Gamma^{\uparrow}\).

 0
lambda_sin float

The constant \(\Lambda_{sin}\) multiplying the \(\sin(4\pi t/\mathcal{T})\) term in \(\Gamma^{\uparrow}\).

 0
extrema_ind Literal[1, 2]

Index of the extremum to diagnose. Use 1 for the first extremum (e.g. maximum) and 2 for the second extremum (e.g. minimum) relative to the first-harmonic solution.

 1
n_year_days int

Number of days in one period \(\mathcal{T}\) (e.g. 360), used to define the annual frequency \(f=1/\mathcal{T}\).

 360
day_seconds float

Duration of one day in seconds, used to convert from days to seconds when defining \(f\).

 86400
numerical bool

If False, use analytic expressions for the coefficients controlling the timing and amplitude of the extrema as functions of the dimensionless parameters. If True, compute these contributions numerically by explicitly solving \(\partial T/\partial\Delta = 0\).

 False
approx_level Optional[Literal['linear', 'linear_phase']]

Approximation level passed through to get_temp_shift_params when computing the extrema numerically.

  • None - use the exact expressions for \(a_1\), \(a_2\), and \(a_3\) (where defined).
  • linear - use the linear approximation in the small dimensionless parameters \(\{S_2/S_1, \lambda_{sq}', \Lambda_{cos}', \Lambda_{sin}'\}\).
  • linear_phase - as linear, but with the \(\lambda_{phase}'\) dependence factored out.
 None
coef_nl_phase bool

If True, compute timing and amplitude coefficients for each mechanism by differentiating \(a_1\), \(a_2\), and \(a_3\) with respect to the dimensional parameters using get_temp_shift_params, and construct both linear and nonlinear (quadratic) phase contributions explicitly. If False, use pre-derived analytic expressions for these coefficients.

 False

Returns:

Name Type Description
final_answer_linear float

The linear-order contribution to the dimensionless extremum shift \(y\) (or equivalently phase shift) obtained by summing all single-parameter contributions.
final_answer_nl float

The total dimensionless extremum shift \(y\) including both linear and nonlinear cross-terms.
info_cont dict

Dictionary containing the individual contributions of each mechanism to the extremum timing. Keys include e.g. 'sw', 'square', 'cos', 'sin' for linear terms and 'nl_sw', 'nl_square', 'nl_cos', 'nl_sin', as well as mixed terms such as name_nl('sw', 'cos') when cross-mechanism nonlinearities are included.
coef dict

Dictionary of analytic coefficients that multiply the dimensionless parameters in info_cont. These describe how each mechanism changes the extremum timing per unit of the corresponding dimensionless parameter.
final_answer_linear_amp float

The linear-order multiplicative factor for the extremum amplitude relative to the reference case, i.e. \(T_{ext}/T_{ext,0}\) including only linear contributions.
final_answer_nl_amp float

The total multiplicative factor for the extremum amplitude relative to the reference case, including both linear and nonlinear contributions.
info_cont_amp dict

Dictionary containing the individual contributions of each mechanism to the extremum amplitude factor, analogous to info_cont but for amplitude rather than timing.
coef_amp dict

Dictionary of analytic coefficients that multiply the dimensionless parameters in info_cont_amp. These describe how each mechanism changes the extremum amplitude per unit of the corresponding dimensionless parameter.
Source code in isca_tools/thesis/surface_energy_budget.py 
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
def get_temp_extrema_theory(heat_capacity: float, sw_amp1: float, sw_amp2: float, lambda_const: float,
                            lambda_phase: float = 0, lambda_sq: float = 0, lambda_cos: float = 0,
                            lambda_sin: float = 0, extrema_ind: Literal[1, 2] = 1,
                            n_year_days: int = 360,
                            day_seconds: float = 86400, numerical: bool = False,
                            approx_level: Optional[Literal['linear', 'linear_phase']] = None,
                            coef_nl_phase: bool = False) -> Tuple[
    float, float, dict, dict, float, float, dict, dict]:
    """
    Find the phase shift and amplitude of the extrema of the analytical temperature solution
    relative to the first-harmonic extremum, i.e. solve for points where $\\partial T/\\partial\\Delta = 0$.

    This function uses the shifted representation
    $T_s(\\Delta) - \\overline{T}_s = \\pm a_1 \\sqrt{1-y^2} + a_2 (1-2y^2) + a_3 y \\sqrt{1-y^2}$
    with $y = \\sin(2\\pi f \\Delta)$, and decomposes the contribution of each empirical
    parameter to the timing and amplitude of the extrema. Depending on the options, it either
    uses analytic coefficients or computes them numerically by solving $\\partial T/\\partial\\Delta = 0$.

    Args:
        heat_capacity:
            $C$, the surface heat capacity in units of $JK^{-1}m^{-2}$.
        sw_amp1:
            $F_1$, the amplitude of the first harmonic of downward shortwave radiation at the surface
            (the $n=1$ Fourier coefficient of $SW^{\\downarrow}$), in units of $Wm^{-2}$.
        sw_amp2:
            $F_2$, the amplitude of the second harmonic of downward shortwave radiation at the surface
            (the $n=2$ Fourier coefficient of $SW^{\\downarrow}$), in units of $Wm^{-2}$.
        lambda_const:
            The linear constant $\\lambda$ used in the approximation for
            $\\Gamma^{\\uparrow} = LW^{\\uparrow} - LW^{\\downarrow} + LH^{\\uparrow} + SH^{\\uparrow}$.
        lambda_phase:
            The constant $\\lambda_{phase}$ multiplying the phase-lag term
            $\\tfrac{1}{2}\\lambda_{phase}(T(t-\\mathcal{T}/4) - T(t+\\mathcal{T}/4))$ in $\\Gamma^{\\uparrow}$.
        lambda_sq:
            The constant $\\lambda_{sq}$ multiplying the quadratic term $-\\lambda_{sq}T'^2$ in $\\Gamma^{\\uparrow}$.
        lambda_cos:
            The constant $\\Lambda_{cos}$ multiplying the $\\cos(4\\pi t/\\mathcal{T})$ term in $\\Gamma^{\\uparrow}$.
        lambda_sin:
            The constant $\\Lambda_{sin}$ multiplying the $\\sin(4\\pi t/\\mathcal{T})$ term in $\\Gamma^{\\uparrow}$.
        extrema_ind:
            Index of the extremum to diagnose. Use `1` for the first extremum (e.g. maximum) and
            `2` for the second extremum (e.g. minimum) relative to the first-harmonic solution.
        n_year_days:
            Number of days in one period $\\mathcal{T}$ (e.g. 360), used to define the annual frequency $f=1/\\mathcal{T}$.
        day_seconds:
            Duration of one day in seconds, used to convert from days to seconds when defining $f$.
        numerical:
            If `False`, use analytic expressions for the coefficients controlling the timing and amplitude
            of the extrema as functions of the dimensionless parameters.
            If `True`, compute these contributions numerically by explicitly solving $\\partial T/\\partial\\Delta = 0$.
        approx_level:
            Approximation level passed through to `get_temp_shift_params` when computing the extrema
            numerically.

            * `None` - use the exact expressions for $a_1$, $a_2$, and $a_3$ (where defined).
            * `linear` - use the linear approximation in the small dimensionless parameters
              $\\{S_2/S_1, \\lambda_{sq}', \\Lambda_{cos}', \\Lambda_{sin}'\\}$.
            * `linear_phase` - as `linear`, but with the $\\lambda_{phase}'$ dependence factored out.
        coef_nl_phase:
            If `True`, compute timing and amplitude coefficients for each mechanism by differentiating
            $a_1$, $a_2$, and $a_3$ with respect to the dimensional parameters using `get_temp_shift_params`,
            and construct both linear and nonlinear (quadratic) phase contributions explicitly.
            If `False`, use pre-derived analytic expressions for these coefficients.

    Returns:
        final_answer_linear:
            The linear-order contribution to the dimensionless extremum shift $y$ (or equivalently phase shift)
            obtained by summing all single-parameter contributions.
        final_answer_nl:
            The total dimensionless extremum shift $y$ including both linear and nonlinear cross-terms.
        info_cont:
            Dictionary containing the individual contributions of each mechanism to the extremum timing.
            Keys include e.g. `'sw'`, `'square'`, `'cos'`, `'sin'` for linear terms and
            `'nl_sw'`, `'nl_square'`, `'nl_cos'`, `'nl_sin'`, as well as mixed terms such as
            `name_nl('sw', 'cos')` when cross-mechanism nonlinearities are included.
        coef:
            Dictionary of analytic coefficients that multiply the dimensionless parameters in `info_cont`.
            These describe how each mechanism changes the extremum timing per unit of the corresponding
            dimensionless parameter.
        final_answer_linear_amp:
            The linear-order multiplicative factor for the extremum amplitude relative to the reference case,
            i.e. $T_{ext}/T_{ext,0}$ including only linear contributions.
        final_answer_nl_amp:
            The total multiplicative factor for the extremum amplitude relative to the reference case,
            including both linear and nonlinear contributions.
        info_cont_amp:
            Dictionary containing the individual contributions of each mechanism to the extremum amplitude factor,
            analogous to `info_cont` but for amplitude rather than timing.
        coef_amp:
            Dictionary of analytic coefficients that multiply the dimensionless parameters in `info_cont_amp`.
            These describe how each mechanism changes the extremum amplitude per unit of the corresponding
            dimensionless parameter.
    """
    f = 1 / (n_year_days * day_seconds)
    x = float(2 * np.pi * f * heat_capacity / lambda_const)
    # Get dimensionless parameters
    param = {'phase': get_param_dimensionless(lambda_phase, heat_capacity=heat_capacity, n_year_days=n_year_days),
             'square': get_param_dimensionless(lambda_sq, sw_fourier_amp1=sw_amp1, lambda_const=lambda_const),
             'cos': get_param_dimensionless(lambda_cos, sw_fourier_amp1=sw_amp1),
             'sin': get_param_dimensionless(lambda_sin, sw_fourier_amp1=sw_amp1),
             'sw': sw_amp2 / sw_amp1}

    if coef_nl_phase:
        def get_coef_from_a(param_name: Literal['sw_amp2', 'lambda_sq', 'lambda_sin']):
            param_use = {'sw_amp2': 0, 'lambda_sq': 0, 'lambda_sin': 0}
            # set dimensionless param to 1 so da/dparam so derivative is equal to a
            if param_name == 'sw_amp2':
                param_use[param_name] = sw_amp1  # so ratio is equal to 1
            elif param_name == 'lambda_sq':
                param_use[param_name] = get_param_dimensionless(1, sw_fourier_amp1=sw_amp1,
                                                                lambda_const=lambda_const, invert=True)
            elif param_name == 'lambda_sin':
                param_use[param_name] = get_param_dimensionless(1, sw_fourier_amp1=sw_amp1, invert=True)
            else:
                raise ValueError(f"Unknown parameter: {param_name}")
            # Use linear here as to compare to other function in extrema_sanity_check notebook
            # makes fairly significant difference if use linear_phase but None and linear are basically the same
            a_1, a_2, a_3 = get_temp_shift_params(heat_capacity, sw_amp1, lambda_const=lambda_const,
                                                  lambda_phase=lambda_phase, lambda_cos=0, n_year_days=n_year_days,
                                                  day_seconds=day_seconds, approx_level='linear',
                                                  **param_use)
            coef_linear = a_3 / a_1
            coef_amp_linear = a_2 / a_1
            if extrema_ind == 2:
                coef_linear *= -1
                coef_amp_linear *= -1
            coef_square = -4 / a_1 ** 2 * a_2 * a_3
            coef_amp_square = 0.5 * coef_linear ** 2
            return coef_linear, coef_square, coef_amp_linear, coef_amp_square

        coef = {}
        coef_amp = {}
        for key, arg in [('sw', 'sw_amp2'), ('square', 'lambda_sq'), ('sin', 'lambda_sin')]:
            var = get_coef_from_a(arg)
            coef[key] = var[0]
            coef[f'nl_{key}'] = var[1]
            coef_amp[key] = var[2]
            coef_amp[f'nl_{key}'] = var[3]

    else:
        # timing
        lambda_ph_mod = param['phase'] / (1 + x ** 2)
        prefactor = 1 / np.sqrt(1 + x ** 2) / (1 + 4 * x ** 2)
        coef = {'sw': prefactor * 4 * x * (x ** 2 - (x ** 4 + 3 * x ** 2 + 1) * lambda_ph_mod),
                'square': prefactor * 4 * x * (1 + x ** 2 * lambda_ph_mod),
                'sin': prefactor * 2 * (3 * x ** 2 + 1 + x ** 2 * (x ** 2 - 1) * lambda_ph_mod)}
        if extrema_ind == 1:
            for key in coef:
                coef[key] *= -1  # first harmonic, all coefficients take negative value
                # Compute coef_amp from coef hence don't change sign

        prefactor_nl = 16 * x / (1 + x ** 2) / (1 + 4 * x ** 2) ** 2
        coef['nl_sw'] = prefactor_nl * (
                x ** 2 * (3 * x ** 2 + 1) - (2 * x ** 6 + 11 * x ** 4 + 6 * x ** 2 + 1) * lambda_ph_mod)
        coef['nl_square'] = -prefactor_nl * (1 + 2 * x ** 2 * lambda_ph_mod)
        coef['nl_sin'] = -coef['nl_sw']
        # Amplitude
        coef_amp = {'sw': -0.5 * coef['sin'], 'square': coef['square'] / (4 * x), 'sin': 0.5 * coef['sw'],
                    'nl_sw': 0.5 * prefactor_nl * x ** 3 * (
                            x ** 2 - 2 * (x ** 4 + 3 * x ** 2 + 1) * lambda_ph_mod),
                    'nl_square': -0.5 * x * coef['nl_square'],
                    'nl_sin': prefactor_nl / 8 / x * (3 * x ** 2 + 1) * (
                            3 * x ** 2 + 1 + 2 * x ** 2 * (x ** 2 - 1) * lambda_ph_mod)}

    coef['cos'] = -coef['sw']
    coef['nl_cos'] = coef['nl_sw']
    coef[name_nl('sw', 'cos')] = -2 * coef['nl_sw']
    coef_amp['cos'] = -coef_amp['sw']
    coef_amp['nl_cos'] = coef_amp['nl_sw']
    coef_amp[name_nl('sw', 'cos')] = -2 * coef_amp['nl_sw']

    info_cont = {}
    info_cont_amp = {}
    if not numerical:
        for key in coef:
            if key == name_nl('sw', 'cos'):
                info_cont[key] = coef[key] * param['sw'] * param['cos']
                info_cont_amp[key] = coef_amp[key] * param['sw'] * param['cos']
            elif 'nl' in key:
                info_cont[key] = coef[key] * param[key.replace('nl_', '')] ** 2
                info_cont_amp[key] = coef_amp[key] * param[key.replace('nl_', '')] ** 2
            else:
                info_cont[key] = coef[key] * param[key]
                info_cont_amp[key] = coef_amp[key] * param[key]
    else:
        def get_y_extrema(param_sw=0, param_square=0, param_cos=0, param_sin=0):
            a_1, a_2, a_3 = get_temp_shift_params(heat_capacity, sw_amp1, param_sw * sw_amp1, lambda_const,
                                                  lambda_phase,
                                                  param_square, param_cos, param_sin, n_year_days=n_year_days,
                                                  day_seconds=day_seconds, approx_level=approx_level)
            # Solve dT/dt=0 explicitly
            if extrema_ind == 2:
                a_1 *= -1  # set negative for first extrema
            dT_dt_func = lambda y: (-a_1 * y - 4 * a_2 * y * np.sqrt(
                1 - y ** 2) + a_3 * (1 - 2 * y ** 2))
            T_func = lambda y: a_1 * np.sqrt(1 - y ** 2) + a_2 * (1 - 2 * y ** 2) + a_3 * y * np.sqrt(1 - y ** 2)
            y_extrema = optimize.least_squares(dT_dt_func, 0, bounds=(-1, 1))['x'][0]
            return y_extrema, T_func(y_extrema) / a_1

        param_ref = {'param_sw': 0, 'param_square': 0, 'param_cos': 0, 'param_sin': 0}
        if approx_level is None:
            param_ref['param_sw'] = 1e-7  # cannot be exactly zero in exact case
        param_with_dim = {'sw': sw_amp2 / sw_amp1, 'square': lambda_sq, 'cos': lambda_cos,
                          'sin': lambda_sin}  # parameters with dimensions
        for key in param_ref:
            key_short = key.replace('param_', '')
            # Initialize all variables as zero
            param_use = copy.deepcopy(param_ref)
            param_use[key] = param_with_dim[key_short]
            info_cont[key_short], info_cont_amp[key_short] = get_y_extrema(**param_use)
            info_cont_amp[key_short] -= 1
        # Adds a nonlinear combination of mechanisms
        # Get non-linear contributions where only two mechanisms are active - include all permutations
        for key1, key2 in itertools.combinations(param_ref, 2):
            key1_short = key1.replace('param_', '')
            key2_short = key2.replace('param_', '')
            param_use = copy.deepcopy(param_ref)
            param_use[key1] = param_with_dim[key1_short]
            param_use[key2] = param_with_dim[key2_short]
            info_cont[name_nl(key1_short, key2_short)], info_cont_amp[name_nl(key1_short, key2_short)] = get_y_extrema(
                **param_use)
            # Subtract the contribution from the linear mechanisms, so only non-linear contribution remains
            info_cont[name_nl(key1_short, key2_short)] -= info_cont[key1_short] + info_cont[key2_short]
            info_cont_amp[name_nl(key1_short, key2_short)] -= info_cont_amp[key1_short] + info_cont_amp[key2_short] + 1

    final_answer_linear = np.asarray(sum([info_cont[key] for key in info_cont if 'nl' not in key]))
    final_answer_nl = np.asarray(sum([info_cont[key] for key in info_cont]))
    final_answer_linear_amp = 1 + np.asarray(sum([info_cont_amp[key] for key in info_cont_amp if 'nl' not in key]))
    final_answer_nl_amp = 1 + np.asarray(sum([info_cont_amp[key] for key in info_cont_amp]))
    if numerical:
        var = get_y_extrema(param_with_dim['sw'], param_with_dim['square'], param_with_dim['cos'],
                            param_with_dim['sin'])
        info_cont['nl_residual'] = var[0] - final_answer_nl
        info_cont_amp['nl_residual'] = var[1] - final_answer_nl_amp
    return final_answer_linear, final_answer_nl, info_cont, coef, final_answer_linear_amp, final_answer_nl_amp, \
        info_cont_amp, coef_amp

get_temp_fourier(time, swdn, heat_capacity, lambda_const, lambda_time_lag=None, lambda_nl=None, n_harmonics=2, include_sw_phase=False, numerical=False, day_seconds=86400, single_harmonic_nl=False, return_anomaly=True)

OUTDATED FUNCTION - NOW USE get_temp_fourier_analytical or get_temp_fourier_numerical

Seeks a fourier solution of the form \(T(t) = \frac{T_0}{2} + \sum_{n=1}^{N} T_n\cos(2n\pi ft - \phi_n)\) to the linearized surface energy budget of the general form:

\[ C\frac{\partial T}{\partial t} = F(t) - \lambda_0 - \sum_{i=1}^{N_{\lambda}}\lambda_i T(t-\Lambda_i) - \sum_{j=2}^{j_{max}}\lambda_{nl_j}\left(T'^{j}(t) - \overline{T'^j}\right) \]

where:

  • \(T' = T(t) - \overline{T}\) is the surface temperature anomaly
  • \(C\) is the heat capacity of the surface
  • \(\overline{T} = T_0/2\) is the mean temperature.
  • \(F(t) = \frac{F_0}{2} + \sum_{n=1}^{N} F_n\cos(2n\pi ft - \varphi_n)\) is the Fourier representation of the downward shortwave radiation at the surface, \(SW^{\downarrow}\).
  • \(\lambda_0 + \sum_{i=1}^{N_{\lambda}}\lambda_i T(t-\Lambda_i) + \sum_{j=2}^{j_{max}}\lambda_{nl_j}\left(T'^{j}(t) - \overline{T'^j}\right)\) is the approximation for \(\Gamma^{\uparrow} = LW^{\uparrow} - LW^{\downarrow} + LH^{\uparrow} + SH^{\uparrow}\)

The solution is exact if \(\lambda_{nl_j}=0 \forall j\) and has the form:

  • \(T_0 = (F_0-2\lambda_0)/\sum_{i=1}^{N_{\lambda}}\lambda_i\)
  • \(T_n = \frac{F_n \cos(\varphi_n)\sqrt{1+\tan^2\phi_n}}{ (2\pi nfC - \sum_i \lambda_i \sin \Phi_{ni})\tan\phi_n + \sum_i \lambda_i \cos \Phi_{ni}}\)
  • \(\tan \phi_n = \frac{2\pi nfC + \tan \varphi_n \sum_i \lambda_i \cos \Phi_{ni} - \sum_i \lambda_i \sin \Phi_{ni}}{-2\pi nfC \tan \varphi_n + \sum_i \lambda_i \cos \Phi_{ni} - \tan \varphi_n \sum_i \lambda_i \sin \Phi_{ni}}\)
  • \(\Phi_{ni} = 2\pi nf \Lambda_i\) (units of radians, whereas \(\Lambda_i\) is in units of time - days).

If \(\lambda_{nl_j}\neq 0\), an approximate numerical solution will be obtained, still of the Fourier form.

Parameters:

Name Type Description Default
time ndarray

float [n_time]
Time in days (assumes periodic e.g. annual mean, so time = np.arange(360)).

 required
swdn ndarray

float [n_time]
Downward shortwave radiation at the surface, \(SW^{\downarrow}\). I.e. swdn_sfc output by Isca. Units: \(Wm^{-2}\).

 required
heat_capacity float

\(C\), the heat capacity of the surface in units of \(JK^{-1}m^{-2}\).
Obtained from mixed layer depth of ocean using get_heat_capacity.

 required
lambda_const ndarray

float [n_lambda+1]
The constants \(\lambda_i\) used in the approximation for \(\Gamma^{\uparrow} = LW^{\uparrow} - LW^{\downarrow} + LH^{\uparrow} + SH^{\uparrow}\).
lambda_const[0] is \(\lambda_0\) and lambda_const[i] is \(\lambda_{i}\) for \(i>0\).

 required
lambda_time_lag Optional[ndarray]

float [n_lambda]
The constants \(\Lambda_i\) used in the approximation for \(\Gamma^{\uparrow}\).
lambda_time_lag[0] is \(\Lambda_1\) and lambda_time_lag[i] is \(\Lambda_{i+1}\) for \(i>0\).

 None
lambda_nl Optional[Union[float, ndarray]]

float [n_lambda_nl] The constants \(\lambda_{nl}\) used in the approximation for \(\Gamma^{\uparrow}\). [0] is squared contribution, [1] is cubed, ...

 None
n_harmonics int

Number of harmonics to use to fit fourier series for both \(T(t)\) and \(F(t)\), \(N\).

 2
include_sw_phase bool

If False, will set all phase factors, \(\varphi_n=0\), in Fourier expansion of \(F(t)\).
These phase factors are usually very small, and it makes the solution for \(T(t)\) more simple if they are set to 0, hence the option.

 False
numerical bool

If True, will compute solution for \(T(t)\) numerically using scipy.optimize.curve_fit rather than using the analytic solution. Will always return numerical solution if lambda_sq \(\neq 0\).

 False
day_seconds float

Duration of a day in seconds.

 86400
single_harmonic_nl bool

If True, the \(\lambda_{nl_j}T'^j\) terms in \(\Gamma^{\uparrow}\) will only use the first harmonic, not all harmonics.

 False
return_anomaly bool

If True, the first return variable, temp_fourier will be the temperature anomaly, i.e. it will not include \(T_0\).

 True

Returns:

Name Type Description
temp_fourier ndarray

float [n_time]
The Fourier series solution that was found for surface temperature, \(T\).
temp_fourier_amp ndarray

float [n_harmonics+1]
The amplitude Fourier coefficients for surface temperature: \(T_n\).
temp_fourier_phase ndarray

float [n_harmonics]
The phase Fourier coefficients for surface temperature: \(\phi_n\).
sw_fourier_amp ndarray

float [n_harmonics+1]
The amplitude Fourier coefficients for shortwave radiation, \(SW^{\downarrow}\): \(F_n\).
sw_fourier_phase ndarray

float [n_harmonics]
The phase Fourier coefficients for shortwave radiation, \(SW^{\downarrow}\): \(\varphi_n\).
Source code in isca_tools/thesis/surface_energy_budget.py 
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
def get_temp_fourier(time: np.ndarray, swdn: np.ndarray, heat_capacity: float,
                     lambda_const: np.ndarray, lambda_time_lag: Optional[np.ndarray] = None,
                     lambda_nl: Optional[Union[float, np.ndarray]] = None, n_harmonics: int = 2,
                     include_sw_phase: bool = False, numerical: bool = False,
                     day_seconds: float = 86400,
                     single_harmonic_nl: bool = False, return_anomaly: bool = True) -> Tuple[
    np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """
    OUTDATED FUNCTION - NOW USE `get_temp_fourier_analytical` or `get_temp_fourier_numerical`

    Seeks a fourier solution of the form $T(t) = \\frac{T_0}{2} + \\sum_{n=1}^{N} T_n\\cos(2n\\pi ft - \\phi_n)$
    to the linearized surface energy budget of the general form:

    $$
    C\\frac{\partial T}{\partial t} = F(t) - \lambda_0 - \sum_{i=1}^{N_{\lambda}}\lambda_i T(t-\Lambda_i) -
    \sum_{j=2}^{j_{max}}\lambda_{nl_j}\\left(T'^{j}(t) - \overline{T'^j}\\right)
    $$

    where:

    * $T' = T(t) - \overline{T}$ is the surface temperature anomaly
    * $C$ is the heat capacity of the surface
    * $\overline{T} = T_0/2$ is the mean temperature.
    * $F(t) = \\frac{F_0}{2} + \\sum_{n=1}^{N} F_n\\cos(2n\\pi ft - \\varphi_n)$ is the Fourier representation
    of the downward shortwave radiation at the surface, $SW^{\downarrow}$.
    * $\lambda_0 + \sum_{i=1}^{N_{\lambda}}\lambda_i T(t-\Lambda_i) +
    \sum_{j=2}^{j_{max}}\lambda_{nl_j}\\left(T'^{j}(t) - \overline{T'^j}\\right)$ is the approximation for
    $\Gamma^{\\uparrow} = LW^{\\uparrow} - LW^{\\downarrow} + LH^{\\uparrow} + SH^{\\uparrow}$

    The solution is exact if $\lambda_{nl_j}=0 \\forall j$ and has the form:

    * $T_0 = (F_0-2\lambda_0)/\sum_{i=1}^{N_{\lambda}}\lambda_i$
    * $T_n = \\frac{F_n \cos(\\varphi_n)\sqrt{1+\\tan^2\phi_n}}{
    (2\pi nfC - \sum_i \lambda_i \sin \Phi_{ni})\\tan\phi_n + \sum_i \lambda_i \cos \Phi_{ni}}$
    * $\\tan \phi_n = \\frac{2\pi nfC + \\tan \\varphi_n \sum_i \lambda_i \cos \Phi_{ni} -
    \sum_i \lambda_i \sin \Phi_{ni}}{-2\pi nfC \\tan \\varphi_n + \sum_i \lambda_i \cos \Phi_{ni} -
    \\tan \\varphi_n \sum_i \lambda_i \sin \Phi_{ni}}$
    * $\Phi_{ni} = 2\pi nf \Lambda_i$ (units of radians, whereas $\Lambda_i$ is in units of time - days).

    If $\lambda_{nl_j}\\neq 0$, an approximate numerical solution will be obtained, still of the Fourier form.

    Args:
        time: `float [n_time]`</br>
            Time in days (assumes periodic e.g. annual mean, so `time = np.arange(360)`).
        swdn: `float [n_time]`</br>
            Downward shortwave radiation at the surface, $SW^{\\downarrow}$. I.e. `swdn_sfc` output by Isca.
            Units: $Wm^{-2}$.
        heat_capacity: $C$, the heat capacity of the surface in units of $JK^{-1}m^{-2}$.</br>
            Obtained from mixed layer depth of ocean using
            [`get_heat_capacity`](../utils/radiation.md#isca_tools.utils.radiation.get_heat_capacity).
        lambda_const: `float [n_lambda+1]`</br>
            The constants $\lambda_i$ used in the approximation for
            $\Gamma^{\\uparrow} = LW^{\\uparrow} - LW^{\\downarrow} + LH^{\\uparrow} + SH^{\\uparrow}$.</br>
            `lambda_const[0]` is $\lambda_0$ and `lambda_const[i]` is $\lambda_{i}$ for $i>0$.
        lambda_time_lag: `float [n_lambda]`</br>
            The constants $\Lambda_i$ used in the approximation for $\Gamma^{\\uparrow}$.</br>
            `lambda_time_lag[0]` is $\Lambda_1$ and `lambda_time_lag[i]` is $\Lambda_{i+1}$ for $i>0$.
        lambda_nl: `float [n_lambda_nl]`
            The constants $\lambda_{nl}$ used in the approximation for $\Gamma^{\\uparrow}$.
            `[0]` is squared contribution, `[1]` is cubed, ...
        n_harmonics: Number of harmonics to use to fit fourier series for both $T(t)$ and $F(t)$, $N$.
        include_sw_phase: If `False`, will set all phase factors, $\\varphi_n=0$, in Fourier expansion of $F(t)$.</br>
            These phase factors are usually very small, and it makes the solution for $T(t)$ more simple if they
            are set to 0, hence the option.
        numerical: If `True`, will compute solution for $T(t)$ numerically using [`scipy.optimize.curve_fit`](
            https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html) rather than using the
            analytic solution. Will always return numerical solution if `lambda_sq` $\\neq 0$.
        day_seconds: Duration of a day in seconds.
        single_harmonic_nl: If `True`, the $\lambda_{nl_j}T'^j$ terms in $\Gamma^{\\uparrow}$ will only
            use the first harmonic, not all harmonics.
        return_anomaly: If `True`, the first return variable, `temp_fourier` will be
            the temperature anomaly, i.e. it will not include $T_0$.

    Returns:
        temp_fourier: `float [n_time]`</br>
            The Fourier series solution that was found for surface temperature, $T$.
        temp_fourier_amp: `float [n_harmonics+1]`</br>
            The amplitude Fourier coefficients for surface temperature: $T_n$.
        temp_fourier_phase: `float [n_harmonics]`</br>
            The phase Fourier coefficients for surface temperature: $\phi_n$.
        sw_fourier_amp: `float [n_harmonics+1]`</br>
            The amplitude Fourier coefficients for shortwave radiation, $SW^{\\downarrow}$: $F_n$.
        sw_fourier_phase: `float [n_harmonics]`</br>
            The phase Fourier coefficients for shortwave radiation, $SW^{\\downarrow}$: $\\varphi_n$.
    """
    if lambda_nl is not None:
        #  deal with case when float given as lambda_nl
        if not hasattr(lambda_nl, "__len__"):
            lambda_nl = np.asarray([lambda_nl])
    n_year_days = len(time)
    n_lambda = len(lambda_const) - 1
    if lambda_time_lag is None:
        lambda_time_lag = np.zeros(n_lambda)
    else:
        if len(lambda_time_lag) != n_lambda:
            raise ValueError(f'Size of lambda_time_lag should be {n_lambda} not {len(lambda_time_lag)}.')

    # Get fourier representation of SW radiation
    sw_fourier_amp, sw_fourier_phase = fourier.get_fourier_fit(time, swdn, n_harmonics)[1:]
    if not include_sw_phase:
        sw_fourier_phase = np.zeros(n_harmonics)
    sw_fourier = fourier.fourier_series(time, sw_fourier_amp, sw_fourier_phase)
    sw_tan = np.tan(sw_fourier_phase)
    sw_cos = np.cos(sw_fourier_phase)
    sw_sin = np.sin(sw_fourier_phase)
    f = 1 / (
            n_year_days * day_seconds)  # must have frequency in units of s^{-1} to deal with phase stuff in radians

    if numerical or (lambda_nl is not None and not single_harmonic_nl):
        if not numerical:
            warnings.warn('Analytic solution not possible with lambda_nl non zero and single_hamrmonic_nl=False')

        def fit_func(time_array, *args):
            fourier_amp_coef = np.asarray([args[i] for i in range(n_harmonics + 1)])
            fourier_phase_coef = np.asarray([args[i] for i in range(n_harmonics + 1, len(args))])
            return swdn_from_temp_fourier(time_array, fourier_amp_coef, fourier_phase_coef, heat_capacity, lambda_const,
                                          lambda_time_lag, lambda_nl, day_seconds, single_harmonic_nl)

        # force positive phase coefficient to match analytic solution
        bounds_lower = [-np.inf] * (n_harmonics + 1) + [0] * n_harmonics
        bounds_upper = [np.inf] * (n_harmonics + 1) + [2 * np.pi] * n_harmonics

        # Starting solution is 1 harmonic analytical solution
        p0 = np.zeros(2 * n_harmonics + 1)
        lambda_phase_const = lambda_time_lag * 2 * np.pi / n_year_days
        lambda_cos = np.sum(lambda_const[1:] * np.cos(lambda_phase_const))
        lambda_sin = np.sum(lambda_const[1:] * np.sin(lambda_phase_const))
        p0[0] = (sw_fourier_amp[0] - 2 * lambda_const[0]) / np.sum(lambda_const[1:])
        p0[n_harmonics + 1] = np.arctan(
            (2 * np.pi * f * heat_capacity + sw_tan[0] * lambda_cos - lambda_sin) / (
                    -2 * np.pi * f * heat_capacity * sw_tan[0] + lambda_cos - sw_tan[0] * lambda_sin))
        p0[1] = sw_fourier_amp[1] * sw_cos[0] / np.cos(p0[n_harmonics + 1]) / (
                (2 * np.pi * f * heat_capacity - lambda_sin) * np.tan(
            p0[n_harmonics + 1]) + lambda_cos)

        try:
            args_found = optimize.curve_fit(fit_func, time, sw_fourier, p0,
                                            bounds=(bounds_lower, bounds_upper))[0]
        except RuntimeError:
            warnings.warn('Hit Runtime Error, trying without bounds')
            args_found = optimize.curve_fit(fit_func, time, sw_fourier, p0)[0]
        temp_fourier_amp = args_found[:n_harmonics + 1]
        temp_fourier_phase = args_found[n_harmonics + 1:]
    else:
        temp_fourier_amp = np.zeros(n_harmonics + 1)
        temp_fourier_phase = np.zeros(n_harmonics)
        temp_fourier_amp[0] = (sw_fourier_amp[0] - 2 * lambda_const[0]) / np.sum(lambda_const[1:])

        for n in range(1, n_harmonics + 1):
            # convert lambda time lag from days to units of radians
            lambda_phase_const = lambda_time_lag * 2 * n * np.pi / n_year_days
            lambda_cos = np.sum(lambda_const[1:] * np.cos(lambda_phase_const))
            lambda_sin = np.sum(lambda_const[1:] * np.sin(lambda_phase_const))

            temp_fourier_phase[n - 1] = np.arctan(
                (2 * np.pi * n * f * heat_capacity + sw_tan[n - 1] * lambda_cos - lambda_sin) / (
                        -2 * np.pi * n * f * heat_capacity * sw_tan[n - 1] + lambda_cos - sw_tan[n - 1] * lambda_sin))
            temp_fourier_amp[n] = sw_fourier_amp[n] * sw_cos[n - 1] / np.cos(
                temp_fourier_phase[n - 1]) / (
                                          (2 * np.pi * n * f * heat_capacity - lambda_sin) * np.tan(
                                      temp_fourier_phase[n - 1]) + lambda_cos)
            if n == 1 and lambda_nl is not None and single_harmonic_nl:
                if len(lambda_nl) > 1:
                    raise ValueError(f'Analytic solution only possible with 1 not {len(lambda_nl)} non-linear terms')
                # Get analytic solution when include squared term in energy budget
                # Assuming squared term dominated by 1st harmonic
                sw_tan[1] = (sw_fourier_amp[2] * sw_sin[1] -
                             0.5 * lambda_nl[0] * temp_fourier_amp[1] ** 2 * np.sin(2 * temp_fourier_phase[0])) / (
                                    sw_fourier_amp[2] * sw_cos[1] -
                                    0.5 * lambda_nl[0] * temp_fourier_amp[1] ** 2 * np.cos(2 * temp_fourier_phase[0]))
                sw_cos[1] -= 0.5 * lambda_nl[0] * temp_fourier_amp[1] ** 2 * np.cos(2 * temp_fourier_phase[0]
                                                                                    ) / sw_fourier_amp[2]
    if return_anomaly:
        temp_fourier = fourier.fourier_series(time, np.append([0], temp_fourier_amp[1:]),
                                              temp_fourier_phase)
    else:
        temp_fourier = fourier.fourier_series(time, temp_fourier_amp, temp_fourier_phase)
    return temp_fourier, temp_fourier_amp, temp_fourier_phase, sw_fourier_amp, sw_fourier_phase

get_temp_fourier_analytic(time, swdn_sfc, heat_capacity, lambda_const, lambda_phase=0, lambda_sq=0, lambda_cos=0, lambda_sin=0, n_harmonics_sw=2, n_harmonics_temp=None, include_sw_phase=False, day_seconds=86400, pad_coefs_phase=False)

Seeks a fourier solution of the form \(T'(t) = \sum_{n=1}^{N} T_n\cos(2n\pi t/\mathcal{T} - \phi_n)\) to the surface energy budget of the general form:

\[ \begin{align} \begin{split} C\frac{\partial T'}{\partial t} = &F(t) - \lambda_0 - \frac{1}{2}\lambda_{phase}(T'(t-\mathcal{T}/4) - T'(t+\mathcal{T}/4)) - \\ &\lambda T'^(t) - \lambda_{sq} T'^{2}(t) - \Lambda_{cos}\cos(4\pi t/\mathcal{T}) - \Lambda_{sin}\sin(4\pi t/\mathcal{T}) \end{split} \end{align} \]

where:

  • \(T' = T(t) - \overline{T}\) is the surface temperature anomaly
  • \(C\) is the heat capacity of the surface
  • \(\overline{T} = T_0/2\) is the mean temperature.
  • \(F(t) = \frac{F_0}{2} + \sum_{n=1}^{N} F_n\cos(2n\pi t/\mathcal{T} - \varphi_n)\) is the Fourier representation of the downward shortwave radiation at the surface, \(SW^{\downarrow}\).
  • \(\lambda_0 + \lambda T' - \lambda_{sq} T'^{2} + \frac{1}{2}\lambda_{phase}(T'(t-\mathcal{T}/4) - T'(t+\mathcal{T}/4)) +\)
    \(\Lambda_{cos}\cos(4\pi t/\mathcal{T}) + \Lambda_{sin}\sin(4\pi t/\mathcal{T})\) is the approximation for \(\Gamma^{\uparrow} = LW^{\uparrow} - LW^{\downarrow} + LH^{\uparrow} + SH^{\uparrow}\).
  • \(\mathcal{T}\) is the period i.e. one year.

The solution is exact if \(\lambda_{sq}=0\) and has the form:

  • \(T_0 = (F_0-2\lambda_0)/\sum_{i=1}^{N_{\lambda}}\lambda_i\)
  • \(T_n = \frac{F_n \cos(\varphi_n)\sqrt{1+\tan^2\phi_n}}{ (2\pi nfC - \sum_i \lambda_i \sin \Phi_{ni})\tan\phi_n + \sum_i \lambda_i \cos \Phi_{ni}}\)
  • \(\tan \phi_n = \frac{2\pi nfC + \tan \varphi_n \sum_i \lambda_i \cos \Phi_{ni} - \sum_i \lambda_i \sin \Phi_{ni}}{-2\pi nfC \tan \varphi_n + \sum_i \lambda_i \cos \Phi_{ni} - \tan \varphi_n \sum_i \lambda_i \sin \Phi_{ni}}\)
  • \(i=1, 2\) with \(\lambda_1=\lambda\) and \(\lambda_2=\lambda_{phase}\).
  • \(\Phi_{n1}=0\) and \(\Phi_{n2} = n\pi/2\) (from \(2n\pi / \mathcal{T} \times \mathcal{T}/4\)).

If \(\lambda_{sq}\neq 0\), an approximate analytical solution will be obtained, assuming \(T'^2(t) \approx (T_1\cos(2\pi ft - \phi_1))^2\).

Parameters:

Name Type Description Default
time ndarray

float [n_time]
Time in days (assumes periodic e.g. annual mean, so time = np.arange(360)).

 required
swdn_sfc ndarray

float [n_time]
Downward shortwave radiation at the surface, \(SW^{\downarrow}\). I.e. swdn_sfc output by Isca. Units: \(Wm^{-2}\).

 required
heat_capacity float

\(C\), the heat capacity of the surface in units of \(JK^{-1}m^{-2}\).
Obtained from mixed layer depth of ocean using get_heat_capacity.

 required
lambda_const float

The constant \(\lambda\) used in the approximation for \(\Gamma^{\uparrow} = LW^{\uparrow} - LW^{\downarrow} + LH^{\uparrow} + SH^{\uparrow}\).

 required
lambda_phase float

The constants \(\lambda_{phase}\) used in the approximation for \(\Gamma^{\uparrow}\).

 0
lambda_sq float

The constant \(\lambda_{sq}\) used in the approximation for \(\Gamma^{\uparrow}\).

 0
lambda_cos float

The constant \(\Lambda_{cos}\) used in the approximation for \(\Gamma^{\uparrow}\).

 0
lambda_sin float

The constant \(\Lambda_{sin}\) used in the approximation for \(\Gamma^{\uparrow}\).

 0
n_harmonics_sw int

Number of harmonics to use to fit fourier series for \(SW^{\downarrow}\). Cannot exceed n_harmonics_temp as extra harmonics would not be used.

 2
n_harmonics_temp Optional[int]

Number, \(N\), of harmonics in fourier solution of temperature anomaly. If not given, will set to n_harmonics_sw.

 None
include_sw_phase bool

If False, will set all phase factors, \(\varphi_n=0\), in Fourier expansion of \(SW^{\downarrow}\).
These phase factors are usually very small, and it makes the solution for \(T'(t)\) more simple if they are set to 0, hence the option.

 False
day_seconds float

Duration of a day in seconds.

 86400
pad_coefs_phase bool

If True, will set temp_fourier_phase and sw_fourier_phase to length n_harmonics+1 to match _fourier_amp, with the first value as zero. Otherwise, it will be size n_harmonics.

 False

Returns:

Name Type Description
temp_fourier ndarray

float [n_time]
The Fourier series solution that was found for surface temperature, \(T\).
temp_fourier_amp ndarray

float [n_harmonics+1]
The amplitude Fourier coefficients for surface temperature: \(T_n\).
temp_fourier_phase ndarray

float [n_harmonics]
The phase Fourier coefficients for surface temperature: \(\phi_n\). If pad_coefs_phase, it will be of the length n_harmonics+1, with the first value as zero.
sw_fourier_amp ndarray

float [n_harmonics+1]
The amplitude Fourier coefficients for shortwave radiation, \(SW^{\downarrow}\): \(F_n\).
sw_fourier_phase ndarray

float [n_harmonics]
The phase Fourier coefficients for shortwave radiation, \(SW^{\downarrow}\): \(\varphi_n\). If pad_coefs_phase, it will be of the length n_harmonics+1, with the first value as zero.
Source code in isca_tools/thesis/surface_energy_budget.py 
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
def get_temp_fourier_analytic(time: np.ndarray, swdn_sfc: np.ndarray, heat_capacity: float,
                              lambda_const: float, lambda_phase: float = 0,
                              lambda_sq: float = 0, lambda_cos: float = 0, lambda_sin: float = 0,
                              n_harmonics_sw: int = 2,
                              n_harmonics_temp: Optional[int] = None,
                              include_sw_phase: bool = False,
                              day_seconds: float = 86400, pad_coefs_phase: bool = False) -> Tuple[
    np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """
    Seeks a fourier solution of the form $T'(t) = \\sum_{n=1}^{N} T_n\\cos(2n\\pi t/\mathcal{T} - \\phi_n)$
    to the surface energy budget of the general form:

    $$
    \\begin{align}
    \\begin{split}
    C\\frac{\partial T'}{\partial t} = &F(t) - \lambda_0 -
    \\frac{1}{2}\lambda_{phase}(T'(t-\mathcal{T}/4) - T'(t+\mathcal{T}/4)) -  \\\\
    &\lambda T'^(t) - \lambda_{sq} T'^{2}(t) - \Lambda_{cos}\\cos(4\\pi t/\mathcal{T}) -
    \Lambda_{sin}\\sin(4\\pi t/\mathcal{T})
    \\end{split}
    \\end{align}
    $$

    where:

    * $T' = T(t) - \overline{T}$ is the surface temperature anomaly
    * $C$ is the heat capacity of the surface
    * $\overline{T} = T_0/2$ is the mean temperature.
    * $F(t) = \\frac{F_0}{2} + \\sum_{n=1}^{N} F_n\\cos(2n\\pi t/\mathcal{T} - \\varphi_n)$ is the Fourier representation
    of the downward shortwave radiation at the surface, $SW^{\downarrow}$.
    * $\lambda_0 + \lambda T' - \lambda_{sq} T'^{2} +
    \\frac{1}{2}\lambda_{phase}(T'(t-\mathcal{T}/4) - T'(t+\mathcal{T}/4)) +$</br>
    $\Lambda_{cos}\\cos(4\\pi t/\mathcal{T}) + \Lambda_{sin}\\sin(4\\pi t/\mathcal{T})$
    is the approximation for $\Gamma^{\\uparrow} = LW^{\\uparrow} - LW^{\\downarrow} + LH^{\\uparrow} + SH^{\\uparrow}$.
    * $\mathcal{T}$ is the period i.e. one year.

    The solution is exact if $\lambda_{sq}=0$ and has the form:

    * $T_0 = (F_0-2\lambda_0)/\sum_{i=1}^{N_{\lambda}}\lambda_i$
    * $T_n = \\frac{F_n \cos(\\varphi_n)\sqrt{1+\\tan^2\phi_n}}{
    (2\pi nfC - \sum_i \lambda_i \sin \Phi_{ni})\\tan\phi_n + \sum_i \lambda_i \cos \Phi_{ni}}$
    * $\\tan \phi_n = \\frac{2\pi nfC + \\tan \\varphi_n \sum_i \lambda_i \cos \Phi_{ni} -
    \sum_i \lambda_i \sin \Phi_{ni}}{-2\pi nfC \\tan \\varphi_n + \sum_i \lambda_i \cos \Phi_{ni} -
    \\tan \\varphi_n \sum_i \lambda_i \sin \Phi_{ni}}$
    * $i=1, 2$ with $\lambda_1=\lambda$ and $\lambda_2=\lambda_{phase}$.
    * $\Phi_{n1}=0$ and $\Phi_{n2} = n\pi/2$ (from $2n\pi / \mathcal{T} \\times \mathcal{T}/4$).

    If $\lambda_{sq}\\neq 0$, an approximate analytical solution will be obtained, assuming
    $T'^2(t) \\approx (T_1\\cos(2\\pi ft - \\phi_1))^2$.

    Args:
        time: `float [n_time]`</br>
            Time in days (assumes periodic e.g. annual mean, so `time = np.arange(360)`).
        swdn_sfc: `float [n_time]`</br>
            Downward shortwave radiation at the surface, $SW^{\\downarrow}$. I.e. `swdn_sfc` output by Isca.
            Units: $Wm^{-2}$.
        heat_capacity: $C$, the heat capacity of the surface in units of $JK^{-1}m^{-2}$.</br>
            Obtained from mixed layer depth of ocean using
            [`get_heat_capacity`](../utils/radiation.md#isca_tools.utils.radiation.get_heat_capacity).
        lambda_const: The constant $\lambda$ used in the approximation for
            $\Gamma^{\\uparrow} = LW^{\\uparrow} - LW^{\\downarrow} + LH^{\\uparrow} + SH^{\\uparrow}$.</br>
        lambda_phase: The constants $\lambda_{phase}$ used in the approximation for $\Gamma^{\\uparrow}$.
        lambda_sq: The constant $\lambda_{sq}$ used in the approximation for $\Gamma^{\\uparrow}$.
        lambda_cos: The constant $\Lambda_{cos}$ used in the approximation for $\Gamma^{\\uparrow}$.
        lambda_sin: The constant $\Lambda_{sin}$ used in the approximation for $\Gamma^{\\uparrow}$.
        n_harmonics_sw: Number of harmonics to use to fit fourier series for $SW^{\\downarrow}$.
            Cannot exceed `n_harmonics_temp` as extra harmonics would not be used.
        n_harmonics_temp: Number, $N$, of harmonics in fourier solution of temperature anomaly. If not given, will
            set to `n_harmonics_sw`.
        include_sw_phase: If `False`, will set all phase factors, $\\varphi_n=0$, in Fourier expansion of
            $SW^{\\downarrow}$.</br>
            These phase factors are usually very small, and it makes the solution for $T'(t)$ more simple if they
            are set to 0, hence the option.
        day_seconds: Duration of a day in seconds.
        pad_coefs_phase: If `True`, will set `temp_fourier_phase` and `sw_fourier_phase`
            to length `n_harmonics+1` to match `_fourier_amp`, with the first value as zero.
            Otherwise, it will be size `n_harmonics`.

    Returns:
        temp_fourier: `float [n_time]`</br>
            The Fourier series solution that was found for surface temperature, $T$.
        temp_fourier_amp: `float [n_harmonics+1]`</br>
            The amplitude Fourier coefficients for surface temperature: $T_n$.
        temp_fourier_phase: `float [n_harmonics]`</br>
            The phase Fourier coefficients for surface temperature: $\phi_n$.
            If `pad_coefs_phase`, it will be of the length `n_harmonics+1`, with the first value as zero.
        sw_fourier_amp: `float [n_harmonics+1]`</br>
            The amplitude Fourier coefficients for shortwave radiation, $SW^{\\downarrow}$: $F_n$.
        sw_fourier_phase: `float [n_harmonics]`</br>
            The phase Fourier coefficients for shortwave radiation, $SW^{\\downarrow}$: $\\varphi_n$.
            If `pad_coefs_phase`, it will be of the length `n_harmonics+1`, with the first value as zero.
    """
    n_year_days = len(time)
    if n_harmonics_temp is None:
        n_harmonics_temp = n_harmonics_sw

    if n_harmonics_temp == 1 and lambda_sq != 0:
        raise ValueError('Cannot solve for non-zero lambda_sq with single harmonic - '
                         'use get_temp_fourier_numerical instead')

    # Get fourier representation of SW radiation
    if n_harmonics_sw > n_harmonics_temp:
        raise ValueError('Cannot have more harmonics for swdn_sfc than have for temperature')
    sw_fourier_amp = np.zeros(n_harmonics_temp + 1)
    sw_fourier_phase = np.zeros(n_harmonics_temp)
    sw_fourier_amp[:n_harmonics_sw + 1], sw_fourier_phase[:n_harmonics_sw] = \
        fourier.get_fourier_fit(time, swdn_sfc, n_harmonics_sw)[1:]
    if not include_sw_phase:
        sw_fourier_phase = np.zeros(n_harmonics_temp)
    sw_tan = np.tan(sw_fourier_phase)
    sw_cos = np.cos(sw_fourier_phase)
    sw_sin = np.sin(sw_fourier_phase)
    if n_harmonics_temp >= 2:
        # For 2nd harmonic, have modification from cos and sin factors
        sw_cos[1] = (sw_fourier_amp[2] * sw_cos[1] - lambda_cos) / sw_fourier_amp[2]
        sw_sin[1] = (sw_fourier_amp[2] * sw_sin[1] - lambda_sin) / sw_fourier_amp[2]
        sw_tan[1] = sw_sin[1] / sw_cos[1]
    f = 1 / (n_year_days * day_seconds)  # must have frequency in units of s^{-1} to deal with phase stuff in radians

    temp_fourier_amp = np.zeros(n_harmonics_temp + 1)  # 1st component will remain 0 so mean of temperature is 0
    temp_fourier_phase = np.zeros(n_harmonics_temp)

    # 0.25 multiplier accounts for the phase delay i.e. shift of 90 days in 360 day year - same as in polyfit_phase
    phase_prefactor = 0.25
    for n in range(1, n_harmonics_temp + 1):
        lambda_term_cos = lambda_const
        lambda_term_sin = 0.5 * lambda_phase * (np.sin(phase_prefactor * 2 * n * np.pi) -
                                                np.sin(-phase_prefactor * 2 * n * np.pi))

        temp_fourier_phase[n - 1] = np.arctan(
            (2 * np.pi * n * f * heat_capacity + sw_tan[n - 1] * lambda_term_cos - lambda_term_sin) / (
                    -2 * np.pi * n * f * heat_capacity * sw_tan[n - 1] + lambda_term_cos - sw_tan[
                n - 1] * lambda_term_sin))
        temp_fourier_amp[n] = sw_fourier_amp[n] * sw_cos[n - 1] / np.cos(
            temp_fourier_phase[n - 1]) / (
                                      (2 * np.pi * n * f * heat_capacity - lambda_term_sin) * np.tan(
                                  temp_fourier_phase[n - 1]) + lambda_term_cos)

        if n == 1 and lambda_sq != 0:
            if include_sw_phase:
                raise ValueError('Solution not possible with sw phase and lambda_sq.')
            # Get analytic solution when include squared term in energy budget
            # Assuming squared term dominated by 1st harmonic
            sw_tan[1] = (sw_fourier_amp[2] * sw_sin[1] -
                         0.5 * lambda_sq * temp_fourier_amp[1] ** 2 * np.sin(2 * temp_fourier_phase[0])) / (
                                sw_fourier_amp[2] * sw_cos[1] -
                                0.5 * lambda_sq * temp_fourier_amp[1] ** 2 * np.cos(2 * temp_fourier_phase[0]))
            if sw_fourier_amp[2] == 0:
                sw_fourier_amp[2] = 1  # this term will cancel out, when do sw_fourier_amp[2]*sw_cos[1]
                # so is just so don't divide by zero below - exact number not important
            sw_cos[1] -= 0.5 * lambda_sq * temp_fourier_amp[1] ** 2 * np.cos(2 * temp_fourier_phase[0]
                                                                             ) / sw_fourier_amp[2]
    temp_fourier = fourier.fourier_series(time, temp_fourier_amp, temp_fourier_phase)
    if pad_coefs_phase:
        temp_fourier_phase = np.hstack((np.zeros(1), temp_fourier_phase))
        sw_fourier_phase = np.hstack((np.zeros(1), sw_fourier_phase))
    return temp_fourier, temp_fourier_amp, temp_fourier_phase, sw_fourier_amp, sw_fourier_phase

get_temp_fourier_analytic2(time, swdn_sfc, heat_capacity, lambda_const, lambda_phase=0, lambda_sq=0, lambda_cos=0, lambda_sin=0, lambda_0=None, n_harmonics=2, day_seconds=86400, pad_coefs_phase=False)

This is the same as get_temp_fourier_analytic but is constrained to only work with n_harmonics=1 and n_harmonics=2 and enforces include_sw_phase=False. Also, it is coded in a way to reflect the algebra more closely.

Seeks a fourier solution of the form \(T(t) = \frac{T_0}{2} + \sum_{n=1}^{2} T_n\cos(2n\pi t/\mathcal{T} - \phi_n)\) to the surface energy budget of the general form:

\[ \begin{align} \begin{split} C\frac{\partial T}{\partial t} = &F(t) - \lambda_0 - \frac{1}{2}\lambda_{phase}(T(t-\mathcal{T}/4) - T(t+\mathcal{T}/4)) - \\ &\lambda T(t) - \lambda_{sq} T'^{2}(t) - \Lambda_{cos}\cos(4\pi t/\mathcal{T}) - \Lambda_{sin}\sin(4\pi t/\mathcal{T}) \end{split} \end{align} \]

where:

  • \(T' = T(t) - \overline{T}\) is the surface temperature anomaly
  • \(C\) is the heat capacity of the surface
  • \(\overline{T} = T_0/2\) is the mean temperature.
  • \(F(t) = \frac{F_0}{2} + \sum_{n=1}^{N} F_n\cos(2n\pi t/\mathcal{T})\) is the Fourier representation of the downward shortwave radiation at the surface, \(SW^{\downarrow}\).
  • \(\lambda_0 + \lambda T - \lambda_{sq} T'^{2} + \frac{1}{2}\lambda_{phase}(T(t-\mathcal{T}/4) - T(t+\mathcal{T}/4)) +\)
    \(\Lambda_{cos}\cos(4\pi t/\mathcal{T}) + \Lambda_{sin}\sin(4\pi t/\mathcal{T})\) is the approximation for \(\Gamma^{\uparrow} = LW^{\uparrow} - LW^{\downarrow} + LH^{\uparrow} + SH^{\uparrow}\).
  • \(\mathcal{T}\) is the period i.e. one year.

The solution is exact if \(\lambda_{sq}=0\), and otherwise approximate assuming first harmonic dominates \(T'^2\) such that \(T'^2(t) \approx (T_1\cos(2\pi ft - \phi_1))^2 + T_2^2/2\):

  • \(T_0 = (F_0-2\lambda_0)/\lambda - \frac{\lambda_{sq}}{\lambda}(T_1^2 + T_2^2)\)
  • \(\tan \phi_1 = x(1 - \frac{\lambda_{phase}}{2\pi fC})\) where \(x=\frac{2\pi fC}{\lambda}\) and \(f=1/\mathcal{T}\).
  • \(T_1 = \frac{F_1}{\lambda}\frac{1}{1+\tan^2\phi_1}\)
  • \(\tan \phi_2 = 2x \frac{1-\frac{1}{2x}\frac{\alpha_2}{1-\alpha_1}}{1+2x\frac{\alpha_2}{1-\alpha_1}}\)
  • \(T_2 = \frac{F_2(1-\Lambda_{cos}')}{\lambda}\frac{1+\tan^2\phi_2}{1+2x\tan\phi_2}(1-\alpha_1)\)

where we use the following dimensionless parameters:

  • \(\Lambda_{cos}' = \Lambda_{cos}/F_1\)
  • \(\Lambda_{sin}' = \Lambda_{sin}/F_1\)
  • \(\lambda_{sq}' = \frac{\lambda_{sq}F_1}{2\lambda^2}\)
  • \(\alpha_1 = \frac{\lambda_{sq}'}{\frac{F_2}{F_1}-\Lambda_{cos}'}\frac{1-\tan^2(\phi_1)}{(1+\tan^2(\phi_1))^2}\)
  • \(\alpha_2 = \frac{\Lambda_{sin}'}{\frac{F_2}{F_1}-\Lambda_{cos}'} + \frac{2\lambda_{sq}'}{\frac{F_2}{F_1}-\Lambda_{cos}'}\frac{\tan(\phi_1)}{(1+\tan^2(\phi_1))^2}\)

Parameters:

Name Type Description Default
time ndarray

float [n_time]
Time in days (assumes periodic e.g. annual mean, so time = np.arange(360)).

 required
swdn_sfc ndarray

float [n_time]
Downward shortwave radiation at the surface, \(SW^{\downarrow}\). I.e. swdn_sfc output by Isca. Units: \(Wm^{-2}\).

 required
heat_capacity float

\(C\), the heat capacity of the surface in units of \(JK^{-1}m^{-2}\).
Obtained from mixed layer depth of ocean using get_heat_capacity.

 required
lambda_const float

The linear constant \(\lambda\) used in the approximation for \(\Gamma^{\uparrow} = LW^{\uparrow} - LW^{\downarrow} + LH^{\uparrow} + SH^{\uparrow}\).

 required
lambda_phase float

The constants \(\lambda_{phase}\) used in the approximation for \(\Gamma^{\uparrow}\).

 0
lambda_sq float

The constant \(\lambda_{sq}\) used in the approximation for \(\Gamma^{\uparrow}\).

 0
lambda_cos float

The constant \(\Lambda_{cos}\) used in the approximation for \(\Gamma^{\uparrow}\).

 0
lambda_sin float

The constant \(\Lambda_{sin}\) used in the approximation for \(\Gamma^{\uparrow}\).

 0
lambda_0 Optional[float]

The constant \(\lambda_0\) used in the approximation for \(\Gamma^{\uparrow}\). Leave as None to set \(T_0=0\) i.e., return the anomaly \(T_s-\overline{T}_s\).

 None
n_harmonics Literal[1, 2]

Number of harmonics to use to fit fourier series for \(SW^{\downarrow}\). Cannot exceed n_harmonics_temp as extra harmonics would not be used.

 2
day_seconds float

Duration of a day in seconds.

 86400
pad_coefs_phase bool

If True, will set temp_fourier_phase and sw_fourier_phase to length n_harmonics+1 to match _fourier_amp, with the first value as zero. Otherwise, it will be size n_harmonics.

 False

Returns:

Name Type Description
temp_fourier ndarray

float [n_time]
The Fourier series solution that was found for surface temperature, \(T\).
temp_fourier_amp ndarray

float [n_harmonics+1]
The amplitude Fourier coefficients for surface temperature: \(T_n\).
temp_fourier_phase ndarray

float [n_harmonics]
The phase Fourier coefficients for surface temperature: \(\phi_n\). If pad_coefs_phase, it will be of the length n_harmonics+1, with the first value as zero.
sw_fourier_amp ndarray

float [n_harmonics+1]
The amplitude Fourier coefficients for shortwave radiation, \(SW^{\downarrow}\): \(F_n\).
sw_fourier_phase ndarray

float [n_harmonics]
The phase Fourier coefficients for shortwave radiation, \(SW^{\downarrow}\): \(\varphi_n\). If pad_coefs_phase, it will be of the length n_harmonics+1, with the first value as zero.
Source code in isca_tools/thesis/surface_energy_budget.py 
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
def get_temp_fourier_analytic2(time: np.ndarray, swdn_sfc: np.ndarray, heat_capacity: float,
                               lambda_const: float, lambda_phase: float = 0,
                               lambda_sq: float = 0, lambda_cos: float = 0, lambda_sin: float = 0,
                               lambda_0: Optional[float] = None,
                               n_harmonics: Literal[1, 2] = 2,
                               day_seconds: float = 86400, pad_coefs_phase: bool = False) -> Tuple[
    np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """
    This is the same as `get_temp_fourier_analytic` but is constrained to only work with `n_harmonics=1` and
    `n_harmonics=2` and enforces `include_sw_phase=False`.
    Also, it is coded in a way to reflect the algebra more closely.

    Seeks a fourier solution of the form $T(t) = \\frac{T_0}{2} + \\sum_{n=1}^{2} T_n\\cos(2n\\pi t/\mathcal{T} - \\phi_n)$
    to the surface energy budget of the general form:

    $$
    \\begin{align}
    \\begin{split}
    C\\frac{\partial T}{\partial t} = &F(t) - \lambda_0 -
    \\frac{1}{2}\lambda_{phase}(T(t-\mathcal{T}/4) - T(t+\mathcal{T}/4)) -  \\\\
    &\lambda T(t) - \lambda_{sq} T'^{2}(t) - \Lambda_{cos}\\cos(4\\pi t/\mathcal{T}) -
    \Lambda_{sin}\\sin(4\\pi t/\mathcal{T})
    \\end{split}
    \\end{align}
    $$

    where:

    * $T' = T(t) - \overline{T}$ is the surface temperature anomaly
    * $C$ is the heat capacity of the surface
    * $\overline{T} = T_0/2$ is the mean temperature.
    * $F(t) = \\frac{F_0}{2} + \\sum_{n=1}^{N} F_n\\cos(2n\\pi t/\mathcal{T})$ is the Fourier representation
    of the downward shortwave radiation at the surface, $SW^{\downarrow}$.
    * $\lambda_0 + \lambda T - \lambda_{sq} T'^{2} +
    \\frac{1}{2}\lambda_{phase}(T(t-\mathcal{T}/4) - T(t+\mathcal{T}/4)) +$</br>
    $\Lambda_{cos}\\cos(4\\pi t/\mathcal{T}) + \Lambda_{sin}\\sin(4\\pi t/\mathcal{T})$
    is the approximation for $\Gamma^{\\uparrow} = LW^{\\uparrow} - LW^{\\downarrow} + LH^{\\uparrow} + SH^{\\uparrow}$.
    * $\mathcal{T}$ is the period i.e. one year.

    The solution is exact if $\lambda_{sq}=0$, and otherwise approximate assuming first harmonic dominates $T'^2$
    such that $T'^2(t) \\approx (T_1\\cos(2\\pi ft - \\phi_1))^2 + T_2^2/2$:

    * $T_0 = (F_0-2\lambda_0)/\lambda - \\frac{\lambda_{sq}}{\lambda}(T_1^2 + T_2^2)$
    * $\\tan \phi_1 = x(1 - \\frac{\lambda_{phase}}{2\pi fC})$ where $x=\\frac{2\pi fC}{\lambda}$ and $f=1/\mathcal{T}$.
    * $T_1 = \\frac{F_1}{\lambda}\\frac{1}{1+\\tan^2\phi_1}$
    * $\\tan \phi_2 = 2x \\frac{1-\\frac{1}{2x}\\frac{\\alpha_2}{1-\\alpha_1}}{1+2x\\frac{\\alpha_2}{1-\\alpha_1}}$
    * $T_2 = \\frac{F_2(1-\Lambda_{cos}')}{\lambda}\\frac{1+\\tan^2\phi_2}{1+2x\\tan\phi_2}(1-\\alpha_1)$

    where we use the following dimensionless parameters:

    * $\Lambda_{cos}' = \Lambda_{cos}/F_1$
    * $\Lambda_{sin}' = \Lambda_{sin}/F_1$
    * $\lambda_{sq}' = \\frac{\lambda_{sq}F_1}{2\lambda^2}$
    * $\\alpha_1 = \\frac{\lambda_{sq}'}{\\frac{F_2}{F_1}-\Lambda_{cos}'}\\frac{1-\\tan^2(\phi_1)}{(1+\\tan^2(\phi_1))^2}$
    * $\\alpha_2 = \\frac{\Lambda_{sin}'}{\\frac{F_2}{F_1}-\Lambda_{cos}'} +
        \\frac{2\lambda_{sq}'}{\\frac{F_2}{F_1}-\Lambda_{cos}'}\\frac{\\tan(\phi_1)}{(1+\\tan^2(\phi_1))^2}$

    Args:
        time: `float [n_time]`</br>
            Time in days (assumes periodic e.g. annual mean, so `time = np.arange(360)`).
        swdn_sfc: `float [n_time]`</br>
            Downward shortwave radiation at the surface, $SW^{\\downarrow}$. I.e. `swdn_sfc` output by Isca.
            Units: $Wm^{-2}$.
        heat_capacity: $C$, the heat capacity of the surface in units of $JK^{-1}m^{-2}$.</br>
            Obtained from mixed layer depth of ocean using
            [`get_heat_capacity`](../utils/radiation.md#isca_tools.utils.radiation.get_heat_capacity).
        lambda_const: The linear constant $\lambda$ used in the approximation for
            $\Gamma^{\\uparrow} = LW^{\\uparrow} - LW^{\\downarrow} + LH^{\\uparrow} + SH^{\\uparrow}$.</br>
        lambda_phase: The constants $\lambda_{phase}$ used in the approximation for $\Gamma^{\\uparrow}$.
        lambda_sq: The constant $\lambda_{sq}$ used in the approximation for $\Gamma^{\\uparrow}$.
        lambda_cos: The constant $\Lambda_{cos}$ used in the approximation for $\Gamma^{\\uparrow}$.
        lambda_sin: The constant $\Lambda_{sin}$ used in the approximation for $\Gamma^{\\uparrow}$.
        lambda_0: The constant $\lambda_0$ used in the approximation for $\Gamma^{\\uparrow}$.
            Leave as `None` to set $T_0=0$ i.e., return the anomaly $T_s-\overline{T}_s$.
        n_harmonics: Number of harmonics to use to fit fourier series for $SW^{\\downarrow}$.
            Cannot exceed `n_harmonics_temp` as extra harmonics would not be used.
        day_seconds: Duration of a day in seconds.
        pad_coefs_phase: If `True`, will set `temp_fourier_phase` and `sw_fourier_phase`
            to length `n_harmonics+1` to match `_fourier_amp`, with the first value as zero.
            Otherwise, it will be size `n_harmonics`.

    Returns:
        temp_fourier: `float [n_time]`</br>
            The Fourier series solution that was found for surface temperature, $T$.
        temp_fourier_amp: `float [n_harmonics+1]`</br>
            The amplitude Fourier coefficients for surface temperature: $T_n$.
        temp_fourier_phase: `float [n_harmonics]`</br>
            The phase Fourier coefficients for surface temperature: $\phi_n$.
            If `pad_coefs_phase`, it will be of the length `n_harmonics+1`, with the first value as zero.
        sw_fourier_amp: `float [n_harmonics+1]`</br>
            The amplitude Fourier coefficients for shortwave radiation, $SW^{\\downarrow}$: $F_n$.
        sw_fourier_phase: `float [n_harmonics]`</br>
            The phase Fourier coefficients for shortwave radiation, $SW^{\\downarrow}$: $\\varphi_n$.
            If `pad_coefs_phase`, it will be of the length `n_harmonics+1`, with the first value as zero.
    """
    n_year_days = len(time)
    f = 1 / (n_year_days * day_seconds)  # must have frequency in units of s^{-1} to deal with phase stuff in radians

    if n_harmonics == 1 and lambda_sq != 0:
        raise ValueError('Cannot solve for non-zero lambda_sq with single harmonic - '
                         'use get_temp_fourier_numerical instead')
    # Get fourier representation of SW radiation. Note I force phase coefs to zero
    sw_fourier_amp = fourier.get_fourier_fit(time, swdn_sfc, n_harmonics)[1]

    temp_fourier_amp = np.zeros(n_harmonics + 1)  # 1st component will remain 0 so mean of temperature is 0
    temp_fourier_phase = np.zeros(n_harmonics + 1)

    # 1st Harmonic
    heat_cap_eff = heat_capacity * (1 - lambda_phase / (2 * np.pi * f * heat_capacity))
    x = 2 * np.pi * f * heat_capacity / lambda_const
    tan_phase1 = 2 * np.pi * f * heat_cap_eff / lambda_const
    temp_fourier_phase[1] = np.arctan(tan_phase1)
    temp_fourier_amp[1] = sw_fourier_amp[1] / lambda_const / np.sqrt(1 + tan_phase1 ** 2)

    # 2nd Harmonic - not exact if lambda_sq!=0, as approx T^2 given by first harmonic squared only
    if n_harmonics == 2:
        sw_amp_ratio = sw_fourier_amp[2] / sw_fourier_amp[1]
        # Put empirical params in dimensionless form
        lambda_cos_dim = get_param_dimensionless(lambda_cos, sw_fourier_amp1=sw_fourier_amp[1])
        lambda_sin_dim = get_param_dimensionless(lambda_sin, sw_fourier_amp1=sw_fourier_amp[1])
        lambda_sq_dim = get_param_dimensionless(lambda_sq, sw_fourier_amp1=sw_fourier_amp[1],
                                                lambda_const=lambda_const)

        # Combine to form other dimensionless factors
        alpha_1 = lambda_sq_dim / (sw_amp_ratio - lambda_cos_dim) * (1 - tan_phase1 ** 2) / (1 + tan_phase1 ** 2) ** 2
        alpha_2 = lambda_sin_dim / (sw_amp_ratio - lambda_cos_dim) + 2 * lambda_sq_dim / (
                sw_amp_ratio - lambda_cos_dim) * tan_phase1 / (
                          1 + tan_phase1 ** 2) ** 2
        phase_mod_factor = (1 - 1 / 2 / x * alpha_2 / (1 - alpha_1)) / (1 + 2 * x * (alpha_2 / (1 - alpha_1)))
        amp_mod_factor = (sw_amp_ratio - lambda_cos_dim) * (1 - alpha_1)

        # Combine usual phase and amp factors with modification factors
        tan_phase2 = 2 * x * phase_mod_factor
        temp_fourier_phase[2] = np.arctan(tan_phase2)
        temp_fourier_amp[2] = sw_fourier_amp[1] / lambda_const * np.sqrt(1 + tan_phase2 ** 2) / (
                1 + 2 * x * tan_phase2) * amp_mod_factor

        # sw_amp2_eff = sw_fourier_amp[2]-lambda_cos
        # lambda_sin_eff = lambda_sin/sw_amp2_eff
        # lambda_sq_eff = lambda_sq/2 * temp_fourier_amp[1]**2 / (2*sw_amp2_eff)
        # alpha_1 = lambda_sq_eff * (1-x**2)/(1+x**2)
        # alpha_2 = lambda_sin_eff + 2 * lambda_sq_eff * x/(1+x**2)
        # phase_mod_factor = (1-lambda_const/(4*np.pi*f*heat_capacity)*(alpha_2/(1-alpha_1)))/(
        #         1+4*np.pi*f*heat_capacity/lambda_const*(alpha_2/(1-alpha_1)))
        # x2 = phase_mod_factor * 4*np.pi*f*heat_capacity/lambda_const        # tan of phase_2
        # temp_fourier_phase[2] = np.arctan(x2)
        # temp_fourier_amp[2] = sw_amp2_eff * np.sqrt(1+x2**2) / (lambda_const + 4*np.pi*f*heat_capacity*x2) * (1-alpha_1)

    # 0th Harmonic
    if lambda_0 is not None:
        temp_fourier_amp[0] = (sw_fourier_amp[0] - 2 * lambda_0) / lambda_const
        if n_harmonics == 2:
            # If include squared term, then get an extra contribution to T_0
            temp_fourier_amp[0] -= lambda_sq / lambda_const * (temp_fourier_amp[1] ** 2 + temp_fourier_amp[2] ** 2)

    temp_fourier = fourier.fourier_series(time, temp_fourier_amp, temp_fourier_phase, pad_coefs_phase=True)
    if not pad_coefs_phase:
        temp_fourier_phase = temp_fourier_phase[1:]
    sw_fourier_phase = np.zeros_like(temp_fourier_phase)
    return temp_fourier, temp_fourier_amp, temp_fourier_phase, sw_fourier_amp, sw_fourier_phase

get_temp_fourier_numerical(time, temp_anom, gamma, swdn_sfc, heat_capacity, n_harmonics_sw=2, n_harmonics_temp=None, deg_gamma_fit=8, phase_gamma_fit=True, resample=False, gamma_fourier_term=False, include_sw_phase=False, day_seconds=86400)

This uses scipy.optimize.curve_fit to numerically seek a fourier solution of the form \(T'(t) = \sum_{n=1}^{N} T_n\cos(2n\pi t/\mathcal{T} - \phi_n)\) to the linearized surface energy budget of the general form:

\[ \begin{align} \begin{split} C\frac{\partial T'}{\partial t} = &SW^{\downarrow}(t) - \lambda_0 - \frac{1}{2}\lambda_{phase}(T'(t-\mathcal{T}/4) - T'(t+\mathcal{T}/4)) - \\ &\sum_{j=1}^{N_{\Gamma}}\lambda_j T'^{j}(t) - \sum_{n=2}^N (\Lambda_{n, cos}\cos(2n\pi t/\mathcal{T}) + \Lambda_{n, sin}\sin(2n\pi t/\mathcal{T})) \end{split} \end{align} \]

where:

  • \(T' = T(t) - \overline{T}\) is the surface temperature anomaly
  • \(C\) is the heat capacity of the surface
  • \(\overline{T} = T_0/2\) is the mean temperature.
  • \(SW^{\downarrow}\) is the downward shortwave radiation at the surface, \(SW^{\downarrow}\).
  • \(\lambda_0 + \frac{1}{2}\lambda_{phase}(T'(t-\mathcal{T}/4) - T'(t+\mathcal{T}/4)) + \sum_{j=1}^{N_{\Gamma}}\lambda_j T'^{j}(t) +\)
    \(\sum_{n=2}^N (\Lambda_{n, cos}\cos(2n\pi t/\mathcal{T}) + \Lambda_{n, sin}\sin(2n\pi t/\mathcal{T}))\) is the approximation for \(\Gamma^{\uparrow} = LW^{\uparrow} - LW^{\downarrow} + LH^{\uparrow} + SH^{\uparrow}\).
  • \(\mathcal{T}\) is the period i.e. one year.

Parameters:

Name Type Description Default
time ndarray

float [n_time]
Time in days (assumes periodic e.g. annual mean, so time = np.arange(360)).

 required
temp_anom ndarray

float [n_time]
Surface temperature anomaly, \(T'(t)\), for each day in time. Used for approximating \(\Gamma^{\uparrow}\).
Assumes periodic so temp[0] is the temperature at time step immediately after time[-1].

 required
gamma ndarray

float [n_time]
Simulated value of \(\Gamma^{\uparrow} = LW^{\uparrow} - LW^{\downarrow} + LH^{\uparrow} + SH^{\uparrow}\) where \(LW\) is longwave, \(LH\) is latent heat and \(SH\) is sensible heat.
Units: \(Wm^{-2}\).

 required
swdn_sfc ndarray

float [n_time]
Downward shortwave radiation at the surface.
Units: \(Wm^{-2}\).

 required
heat_capacity float

\(C\), the heat capacity of the surface in units of \(JK^{-1}m^{-2}\).
Obtained from mixed layer depth of ocean using get_heat_capacity.

 required
n_harmonics_sw int

Number of harmonics to use to fit fourier series for \(SW^{\downarrow}\). Cannot exceed n_harmonics_temp as extra harmonics would not be used. Set to None, to use no approximation for \(SW^{\downarrow}\) - but we weary with comparing to analytic solution in this case.

 2
n_harmonics_temp Optional[int]

Number, \(N\), of harmonics in fourier solution of temperature anomaly. If not given, will set to n_harmonics_sw.

 None
deg_gamma_fit int

Power, \(N_{\Gamma}\), to go up to in polyomial approximation of \(\Gamma^{\uparrow}\) seeked.

 8
phase_gamma_fit bool

If False will set \(\lambda_{phase}=0\). Otherwise, will use polyfit_phase to estimate it.

 True
resample bool

If True, will use resample_data to make data evenly spaced in \(x\) before calling np.polyfit, when obtaining \(\Gamma\) coefficients.

 False
gamma_fourier_term bool

Whether to fit the Fourier contribution \(\sum_{n=2}^N (\Lambda_{n, cos}\cos(2n\pi t/\mathcal{T}) + \Lambda_{n, sin}\sin(2n\pi t/\mathcal{T}))\) to \(\Gamma^{\uparrow}\) with fourier_harmonics=np.arange(2, n_harmonics+1) in polyfit_phase. Idea behind this is to account for contribution of \(\Gamma^{\uparrow}\) that is not temperature dependent.

 False
include_sw_phase bool

If False, will set all phase factors, \(\varphi_n=0\), in Fourier expansion of \(SW^{\downarrow}\).
These phase factors are usually very small, and it makes the analytic solution for \(T'(t)\) more simple if they are set to 0, hence the option. Only use if n_harmonics_sw not None.

 False
day_seconds float

Duration of a day in seconds.

 86400

Returns:

Name Type Description
ndarray

temp_fourier float [n_time]
The Fourier series solution that was found for surface temperature anomaly.
Units: \(K\).
temp_fourier_amp ndarray

float [n_harmonics+1]
The amplitude Fourier coefficients for surface temperature: \(T_n\).
temp_fourier_phase ndarray

float [n_harmonics]
The phase Fourier coefficients for surface temperature: \(\phi_n\).
Source code in isca_tools/thesis/surface_energy_budget.py 
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
def get_temp_fourier_numerical(time: np.ndarray, temp_anom: np.ndarray, gamma: np.ndarray,
                               swdn_sfc: np.ndarray, heat_capacity: float,
                               n_harmonics_sw: int = 2, n_harmonics_temp: Optional[int] = None,
                               deg_gamma_fit: int = 8, phase_gamma_fit: bool = True,
                               resample: bool = False,
                               gamma_fourier_term: bool = False,
                               include_sw_phase: bool = False,
                               day_seconds: float = 86400) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """
    This uses [`scipy.optimize.curve_fit`](
    https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html) to numerically seek
    a fourier solution of the form $T'(t) = \\sum_{n=1}^{N} T_n\\cos(2n\\pi t/\mathcal{T} - \\phi_n)$
    to the linearized surface energy budget of the general form:

    $$
    \\begin{align}
    \\begin{split}
    C\\frac{\partial T'}{\partial t} = &SW^{\\downarrow}(t) - \lambda_0 -
    \\frac{1}{2}\lambda_{phase}(T'(t-\mathcal{T}/4) - T'(t+\mathcal{T}/4)) -  \\\\
    &\\sum_{j=1}^{N_{\Gamma}}\lambda_j T'^{j}(t) - \\sum_{n=2}^N (\Lambda_{n, cos}\\cos(2n\\pi t/\mathcal{T}) +
    \Lambda_{n, sin}\\sin(2n\\pi t/\mathcal{T}))
    \\end{split}
    \\end{align}
    $$

    where:

    * $T' = T(t) - \overline{T}$ is the surface temperature anomaly
    * $C$ is the heat capacity of the surface
    * $\overline{T} = T_0/2$ is the mean temperature.
    * $SW^{\downarrow}$ is the downward shortwave radiation at the surface, $SW^{\downarrow}$.
    * $\lambda_0 + \\frac{1}{2}\lambda_{phase}(T'(t-\mathcal{T}/4) - T'(t+\mathcal{T}/4)) +
    \\sum_{j=1}^{N_{\Gamma}}\lambda_j T'^{j}(t) +$</br>
    $\\sum_{n=2}^N (\Lambda_{n, cos}\\cos(2n\\pi t/\mathcal{T}) + \Lambda_{n, sin}\\sin(2n\\pi t/\mathcal{T}))$
    is the approximation for $\Gamma^{\\uparrow} = LW^{\\uparrow} - LW^{\\downarrow} + LH^{\\uparrow} + SH^{\\uparrow}$.
    * $\mathcal{T}$ is the period i.e. one year.


    Args:
        time: `float [n_time]`</br>
            Time in days (assumes periodic e.g. annual mean, so `time = np.arange(360)`).
        temp_anom: `float [n_time]`</br>
            Surface temperature anomaly, $T'(t)$, for each day in `time`. Used for approximating
            $\Gamma^{\\uparrow}$.</br>
            Assumes periodic so temp[0] is the temperature at time step immediately after `time[-1]`.
        gamma: `float [n_time]`</br>
            Simulated value of $\Gamma^{\\uparrow} = LW^{\\uparrow} - LW^{\\downarrow} + LH^{\\uparrow} + SH^{\\uparrow}$
            where $LW$ is longwave, $LH$ is latent heat and $SH$ is sensible heat.</br>
            Units: $Wm^{-2}$.
        swdn_sfc: `float [n_time]`</br>
            Downward shortwave radiation at the surface.</br>
            Units: $Wm^{-2}$.
        heat_capacity: $C$, the heat capacity of the surface in units of $JK^{-1}m^{-2}$.</br>
            Obtained from mixed layer depth of ocean using
            [`get_heat_capacity`](../utils/radiation.md#isca_tools.utils.radiation.get_heat_capacity).
        n_harmonics_sw: Number of harmonics to use to fit fourier series for $SW^{\\downarrow}$.
            Cannot exceed `n_harmonics_temp` as extra harmonics would not be used.</n>
            Set to None, to use no approximation for $SW^{\\downarrow}$ - but we weary with comparing to analytic
            solution in this case.
        n_harmonics_temp: Number, $N$, of harmonics in fourier solution of temperature anomaly. If not given, will
            set to `n_harmonics_sw`.
        deg_gamma_fit: Power, $N_{\Gamma}$, to go up to in polyomial approximation of $\Gamma^{\\uparrow}$ seeked.
        phase_gamma_fit: If `False` will set $\lambda_{phase}=0$.
            Otherwise, will use [`polyfit_phase`](../utils/numerical.md#isca_tools.utils.numerical.polyfit_phase)
            to estimate it.
        resample: If `True`, will use [`resample_data`](../utils/numerical.md#isca_tools.utils.numerical.resample_data)
            to make data evenly spaced in $x$ before calling `np.polyfit`, when obtaining $\Gamma$ coefficients.
        gamma_fourier_term: Whether to fit the Fourier contribution
            $\\sum_{n=2}^N (\Lambda_{n, cos}\\cos(2n\\pi t/\mathcal{T}) + \Lambda_{n, sin}\\sin(2n\\pi t/\mathcal{T}))$
             to $\Gamma^{\\uparrow}$ with `fourier_harmonics=np.arange(2, n_harmonics+1)` in
            [`polyfit_phase`](../utils/numerical.md#isca_tools.utils.numerical.polyfit_phase). Idea behind this
            is to account for contribution of $\Gamma^{\\uparrow}$ that is not temperature dependent.
        include_sw_phase: If `False`, will set all phase factors, $\\varphi_n=0$, in Fourier expansion of
            $SW^{\\downarrow}$.</br>
            These phase factors are usually very small, and it makes the analytic solution for $T'(t)$ more simple
            if they are set to 0, hence the option. Only use if `n_harmonics_sw` not `None`.
        day_seconds: Duration of a day in seconds.

    Returns:
        temp_fourier `float [n_time]`</br>
            The Fourier series solution that was found for surface temperature anomaly.</br>
            Units: $K$.
        temp_fourier_amp: `float [n_harmonics+1]`</br>
            The amplitude Fourier coefficients for surface temperature: $T_n$.
        temp_fourier_phase: `float [n_harmonics]`</br>
            The phase Fourier coefficients for surface temperature: $\phi_n$.
    """
    if n_harmonics_temp is None:
        n_harmonics_temp = n_harmonics_sw
    if gamma_fourier_term and phase_gamma_fit:
        gamma_approx_coefs, gamma_fourier_term_coefs_amp, gamma_fourier_term_coefs_phase = \
            numerical.polyfit_phase(temp_anom, gamma, deg_gamma_fit, resample=resample,
                                    include_phase=phase_gamma_fit, fourier_harmonics=np.arange(2, n_harmonics_temp + 1))
    else:
        gamma_approx_coefs = numerical.polyfit_phase(temp_anom, gamma, deg_gamma_fit, resample=resample,
                                                     include_phase=phase_gamma_fit)
        gamma_fourier_term_coefs_amp = None
        gamma_fourier_term_coefs_phase = None

    def fit_func(time_array, *args):
        # first n_harmonic values in args are the temperature amplitude coefficients (excluding T_0)
        # last n_harmonic values in args are the temperature phase coefficients
        fourier_amp_coef = np.asarray([args[i] for i in range(n_harmonics_temp)])
        # make first coefficient 0 so gives anomaly
        fourier_amp_coef = np.append([0], fourier_amp_coef)
        fourier_phase_coef = np.asarray([args[i] for i in range(n_harmonics_temp, len(args))])
        temp_anom_fourier = fourier.fourier_series(time_array, fourier_amp_coef, fourier_phase_coef)
        dtemp_dt_fourier = fourier.fourier_series_deriv(time_array, fourier_amp_coef,
                                                        fourier_phase_coef, day_seconds)
        if phase_gamma_fit:
            gamma_approx = numerical.polyval_phase(gamma_approx_coefs, temp_anom_fourier,
                                                   coefs_fourier_amp=gamma_fourier_term_coefs_amp,
                                                   coefs_fourier_phase=gamma_fourier_term_coefs_phase)
        else:
            gamma_approx = np.polyval(gamma_approx_coefs, temp_anom_fourier)
        return heat_capacity * dtemp_dt_fourier + gamma_approx

    # Starting guess is linear 1 harmonic linear analytical solution given gamma=lambda_const_guess*temp
    # so only 1st harmonic coefficients of amplitude and phase are needed, rest are set to zero
    p0 = np.zeros(2 * n_harmonics_temp)
    f = 1 / (time.size * day_seconds)
    p0[n_harmonics_temp] = np.arctan((2 * np.pi * f * heat_capacity) / gamma_approx_coefs[-2])
    if n_harmonics_sw is None:
        sw_fourier_amp = fourier.get_fourier_fit(time, swdn_sfc, 1)[1]
        # find temperature solution which minimises error to full insolation, no fourier approx
        sw_fourier_fit = swdn_sfc
    else:
        if include_sw_phase:
            sw_fourier_fit, sw_fourier_amp = fourier.get_fourier_fit(time, swdn_sfc, n_harmonics_sw)[:2]
        else:
            sw_fourier_amp, sw_fourier_phase = fourier.get_fourier_fit(time, swdn_sfc, n_harmonics_sw)[1:]
            sw_fourier_fit = fourier.fourier_series(time, sw_fourier_amp, sw_fourier_phase * 0)
    p0[0] = sw_fourier_amp[1] / np.cos(p0[n_harmonics_temp]) / (
            (2 * np.pi * f * heat_capacity) * np.tan(p0[n_harmonics_temp]) + gamma_approx_coefs[-2])
    args_found = optimize.curve_fit(fit_func, time, sw_fourier_fit, p0)[0]
    temp_fourier_amp = np.append([0], args_found[:n_harmonics_temp])
    temp_fourier_phase = args_found[n_harmonics_temp:]
    return fourier.fourier_series(time, temp_fourier_amp, temp_fourier_phase), temp_fourier_amp, temp_fourier_phase

get_temp_shift_params(heat_capacity, sw_amp1, sw_amp2, lambda_const, lambda_phase, lambda_sq, lambda_cos, lambda_sin, n_year_days=360, day_seconds=86400, approx_level=None)

If shifting the solution to get_temp_fourier_analytic2 with two harmonics in time from \(T_s(t)\) to \(T_s(t_1+\Delta)\) where \(t_1\) is the extrema of the first harmonic, this returns the values of \(a_1\), \(a_2\), and \(a_3\) such that:

$$ T_s(\Delta) - \overline{T}_s = \pm a_1 \sqrt{1-y^2} + a_2 (1-2y^2) + a_3 y \sqrt{1-y^2} $$ where \(y=\sin(2\pi f\Delta)\).

Parameters:

Name Type Description Default
heat_capacity Union[ndarray, DataArray, float]

\(C\), the surface heat capacity, in units of \(JK^{-1}m^{-2}\). Typically obtained from a mixed-layer depth via get_heat_capacity.

 required
sw_amp1 Union[ndarray, DataArray, float]

\(F_1\), the amplitude of the first harmonic of downward shortwave radiation at the surface, i.e. the \(n=1\) Fourier coefficient of \(SW^{\downarrow}\) in units of \(Wm^{-2}\).

 required
sw_amp2 Union[ndarray, DataArray, float]

\(F_2\), the amplitude of the second harmonic of downward shortwave radiation at the surface, i.e. the \(n=2\) Fourier coefficient of \(SW^{\downarrow}\) in units of \(Wm^{-2}\).

 required
lambda_const Union[ndarray, DataArray, float]

The linear constant \(\lambda\) used in the approximation for \(\Gamma^{\uparrow} = LW^{\uparrow} - LW^{\downarrow} + LH^{\uparrow} + SH^{\uparrow}\).

 required
lambda_phase Union[ndarray, DataArray, float]

The constant \(\lambda_{phase}\) multiplying the phase-lag term \(\tfrac{1}{2}\lambda_{phase}(T(t-\mathcal{T}/4) - T(t+\mathcal{T}/4))\) in \(\Gamma^{\uparrow}\).

 required
lambda_sq Union[ndarray, DataArray, float]

The constant \(\lambda_{sq}\) multiplying the quadratic term \(-\lambda_{sq}T'^2\) in \(\Gamma^{\uparrow}\).

 required
lambda_cos Union[ndarray, DataArray, float]

The constant \(\Lambda_{cos}\) multiplying the \(\cos(4\pi t/\mathcal{T})\) term in \(\Gamma^{\uparrow}\).

 required
lambda_sin Union[ndarray, DataArray, float]

The constant \(\Lambda_{sin}\) multiplying the \(\sin(4\pi t/\mathcal{T})\) term in \(\Gamma^{\uparrow}\).

 required
n_year_days int

Number of days in one period \(\mathcal{T}\) (e.g. 360), used to define the annual frequency \(f=1/\mathcal{T}\).

 360
day_seconds float

Duration of one day in seconds, used to convert from days to seconds when defining \(f\).

 86400
approx_level Optional[Literal['linear', 'linear_phase']]

Three options:

  • None - will return the exact values. Will produce an error if \(S_2 = \Lambda_{\cos}=0\), which does not happen with the other approximations.
  • linear - \(a_2\) and \(a_3\) will be in the form \(\sum_{\chi} \gamma(\lambda, C, S_1, \lambda_{ph}')\chi\) with \(\chi \in \{\frac{S_2}{S_1}, \lambda_{sq}', \Lambda_{\cos}', \Lambda_{\sin}'\}\). Not exact but very good approximation.
  • linear_phase - \(a_2\) and \(a_3\) will be in the form \(\sum_{\chi} \gamma(\lambda, C, S_1) (1+\gamma_2(\lambda, C)\lambda_{ph}')\chi\) i.e., will extract the \(\lambda_{ph}'\) dependence.
 None

Returns:

Name Type Description
a_1 Union[ndarray, DataArray, float]

The parameter \(a_1\) controlling the leading-order dependence on \(\sqrt{1-y^2}\).
a_2 Union[ndarray, DataArray, float]

The parameter \(a_2\) controlling the \((1-2y^2)\) contribution.
a_3 Union[ndarray, DataArray, float]

The parameter \(a_3\) controlling the \(y\sqrt{1-y^2}\) contribution.
Source code in isca_tools/thesis/surface_energy_budget.py 
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
def get_temp_shift_params(heat_capacity: Union[np.ndarray, xr.DataArray, float],
                          sw_amp1: Union[np.ndarray, xr.DataArray, float],
                          sw_amp2: Union[np.ndarray, xr.DataArray, float],
                          lambda_const: Union[np.ndarray, xr.DataArray, float],
                          lambda_phase: Union[np.ndarray, xr.DataArray, float],
                          lambda_sq: Union[np.ndarray, xr.DataArray, float],
                          lambda_cos: Union[np.ndarray, xr.DataArray, float],
                          lambda_sin: Union[np.ndarray, xr.DataArray, float],
                          n_year_days: int = 360,
                          day_seconds: float = 86400,
                          approx_level: Optional[Literal['linear', 'linear_phase']] = None
                          ) -> Tuple[Union[np.ndarray, xr.DataArray, float], Union[np.ndarray, xr.DataArray, float],
Union[np.ndarray, xr.DataArray, float]]:
    """
    If shifting the solution to `get_temp_fourier_analytic2` with two harmonics in time from $T_s(t)$ to
    $T_s(t_1+\Delta)$ where $t_1$ is the extrema of the first harmonic, this returns the values of $a_1$, $a_2$,
    and $a_3$ such that:

    $$
    T_s(\Delta) - \overline{T}_s = \pm a_1 \sqrt{1-y^2} + a_2 (1-2y^2) + a_3 y \sqrt{1-y^2}
    $$
    where $y=\sin(2\pi f\Delta)$.

    Args:
        heat_capacity: $C$, the surface heat capacity, in units of $JK^{-1}m^{-2}$.
            Typically obtained from a mixed-layer depth via `get_heat_capacity`.
        sw_amp1:
            $F_1$, the amplitude of the first harmonic of downward shortwave radiation at the surface,
            i.e. the $n=1$ Fourier coefficient of $SW^{\downarrow}$ in units of $Wm^{-2}$.
        sw_amp2:
            $F_2$, the amplitude of the second harmonic of downward shortwave radiation at the surface,
            i.e. the $n=2$ Fourier coefficient of $SW^{\downarrow}$ in units of $Wm^{-2}$.
        lambda_const:
            The linear constant $\lambda$ used in the approximation for
            $\Gamma^{\\uparrow} = LW^{\\uparrow} - LW^{\\downarrow} + LH^{\\uparrow} + SH^{\\uparrow}$.
        lambda_phase:
            The constant $\lambda_{phase}$ multiplying the phase-lag term
            $\\tfrac{1}{2}\lambda_{phase}(T(t-\mathcal{T}/4) - T(t+\mathcal{T}/4))$ in $\Gamma^{\\uparrow}$.
        lambda_sq:
            The constant $\lambda_{sq}$ multiplying the quadratic term $-\\lambda_{sq}T'^2$ in $\Gamma^{\\uparrow}$.
        lambda_cos:
            The constant $\Lambda_{cos}$ multiplying the $\cos(4\pi t/\mathcal{T})$ term in $\\Gamma^{\\uparrow}$.
        lambda_sin:
            The constant $\Lambda_{sin}$ multiplying the $\sin(4\pi t/\mathcal{T})$ term in $\Gamma^{\\uparrow}$.
        n_year_days:
            Number of days in one period $\mathcal{T}$ (e.g. 360), used to define the annual frequency $f=1/\mathcal{T}$.
        day_seconds:
            Duration of one day in seconds, used to convert from days to seconds when defining $f$.
        approx_level: Three options:

            * `None` - will return the exact values. Will produce an error if $S_2 = \Lambda_{\cos}=0$, which
            does not happen with the other approximations.
            * `linear` - $a_2$ and $a_3$ will be in the form $\sum_{\chi} \gamma(\lambda, C, S_1, \lambda_{ph}')\chi$
            with $\chi \in \{\\frac{S_2}{S_1}, \lambda_{sq}', \Lambda_{\cos}', \Lambda_{\sin}'\}$. Not exact but
            very good approximation.
            * `linear_phase` - $a_2$ and $a_3$ will be in the form $\sum_{\chi} \gamma(\lambda, C, S_1)
            (1+\gamma_2(\lambda, C)\lambda_{ph}')\chi$ i.e., will extract the $\lambda_{ph}'$ dependence.

    Returns:
        a_1: The parameter $a_1$ controlling the leading-order dependence on $\sqrt{1-y^2}$.
        a_2: The parameter $a_2$ controlling the $(1-2y^2)$ contribution.
        a_3: The parameter $a_3$ controlling the $y\sqrt{1-y^2}$ contribution.
    """
    f = 1 / (n_year_days * day_seconds)
    x = 2 * np.pi * f * heat_capacity / lambda_const
    # Make params dimensionless
    lambda_cos = get_param_dimensionless(lambda_cos, sw_fourier_amp1=sw_amp1)
    lambda_sin = get_param_dimensionless(lambda_sin, sw_fourier_amp1=sw_amp1)
    lambda_sq = get_param_dimensionless(lambda_sq, sw_fourier_amp1=sw_amp1, lambda_const=lambda_const)
    lambda_ph = get_param_dimensionless(lambda_phase, heat_capacity=heat_capacity, n_year_days=n_year_days)

    x1 = x * (1 - lambda_ph)
    sw_amp_ratio = sw_amp2 / sw_amp1
    sw_amp_ratio_mod = sw_amp_ratio - lambda_cos
    a_1 = sw_amp1 / lambda_const / np.sqrt(1 + x1 ** 2)

    if approx_level is None:
        if _all_zero(sw_amp_ratio_mod):
            raise ValueError('Exact solution not possible with sw_amp2 and lambda_cos both zero')
        # Get alpha_1 and alpha_2 to compute x2 - copied from `get_temp_fourier_analytic2` function
        alpha_1 = lambda_sq / sw_amp_ratio_mod * (1 - x1 ** 2) / (1 + x1 ** 2) ** 2
        alpha_2 = lambda_sin / sw_amp_ratio_mod + 2 * lambda_sq / sw_amp_ratio_mod * x1 / (1 + x1 ** 2) ** 2
        phase_mod_factor = (1 - 1 / 2 / x * alpha_2 / (1 - alpha_1)) / (1 + 2 * x * (alpha_2 / (1 - alpha_1)))
        # Combine usual phase and amp factors with modification factors
        x2 = 2 * x * phase_mod_factor
        prefactor = sw_amp1 / lambda_const / (1 + x1 ** 2) / (1 + 2 * x * x2) * sw_amp_ratio_mod * (1 - alpha_1)
        a_2 = prefactor * (1 - x1 ** 2 + 2 * x1 * x2)
        a_3 = -2 * prefactor * (2 * x1 - x2 + x1 ** 2 * x2)
    elif approx_level == 'linear':
        prefactor = sw_amp1 / lambda_const / (1 + x1 ** 2) / (1 + 4 * x ** 2)
        a_2 = prefactor * (
                sw_amp_ratio_mod * (4 * x1 * x - x1 ** 2 + 1) - lambda_sin * 2 * (x * x1 ** 2 + x1 - x) - lambda_sq)
        a_3 = -2 * prefactor * (
                sw_amp_ratio_mod * 2 * (x1 - x + x * x1 ** 2) + lambda_sin * (
                4 * x * x1 - x1 ** 2 + 1) + lambda_sq * 2 * x)
    elif approx_level == 'linear_phase':
        prefactor = sw_amp1 / lambda_const / (1 + x ** 2) / (1 + 4 * x ** 2)
        lambda_ph_mod = lambda_ph / (1 + x ** 2)
        # Approx 1/a_1 and then invert as 1/a_1 used elsewhere
        a_1 = 1 / (lambda_const / sw_amp1 * np.sqrt(1 + x ** 2) * (1 - x ** 2 * lambda_ph_mod))
        a_2 = prefactor * ((3 * x ** 2 + 1 + 4 * x ** 4 * lambda_ph_mod) * sw_amp_ratio_mod -
                           2 * x * (x ** 2 - (3 * x ** 2 + 1) * lambda_ph_mod) * lambda_sin -
                           (1 + 2 * x ** 2 * lambda_ph_mod) * lambda_sq)
        a_3 = -2 * prefactor * ((2 * x ** 3 - 2 * x * (3 * x ** 2 + 1) * lambda_ph_mod) * sw_amp_ratio_mod +
                                (3 * x ** 2 + 1 + 4 * x ** 4 * lambda_ph_mod) * lambda_sin +
                                2 * x * (1 + 2 * x ** 2 * lambda_ph_mod) * lambda_sq)
    else:
        raise ValueError(f"Unknown approx_level: {approx_level}")
    return a_1, a_2, a_3

phase_coef_conversion(coef_linear, coef_phase, to_time=True, pos_amp=False, neg_amp=False, take_linear_sign=True)

Method to convert between interpretation of phase empirical coefficients. If to_time=True, expect input to be \(\lambda\) and \(\lambda_{ph}\), and will return \(\lambda_{mod}\), and \(2\pi ft_{ph}\). If to_time=False, expect the opposite.

Conversion performed according to:

\[\begin{align} \lambda &= \lambda_{\text{mod}}\cos(2 \pi ft_{\text{ph}}) \\ \lambda_{\text{ph}} &= \lambda_{\text{mod}}\sin(2 \pi ft_{\text{ph}}) \end{align}\]

Parameters:

Name Type Description Default
coef_linear Union[float, ndarray, DataArray]

Cosine/linear-amplitude coefficient associated with the in-phase component. If to_time=True this is interpreted as \(\lambda\); if to_time=False this is interpreted as \(\lambda_{\mathrm{mod}}\). Can be a scalar, numpy.ndarray, or xarray.DataArray.

 required
coef_phase Union[float, ndarray, DataArray]

Phase-related coefficient. If to_time=True this is interpreted as \(\lambda_{\mathrm{ph}}\) (sine coefficient); if to_time=False this is interpreted as the phase angle \(2\pi f t_{\mathrm{ph}}\) (radians) used inside \(\cos(\cdot)\) / \(\sin(\cdot)\). Can be a scalar, numpy.ndarray, or xarray.DataArray.

 required
to_time bool

Direction of conversion. If True, convert \((\lambda,\lambda_{\mathrm{ph}})\rightarrow (\lambda_{\mathrm{mod}},2\pi f t_{\mathrm{ph}})\). If False, convert \((\lambda_{\mathrm{mod}},2\pi f t_{\mathrm{ph}}) \rightarrow (\lambda,\lambda_{\mathrm{ph}})\).

 True
pos_amp bool

Only makes a difference if to_time=True.
If pos_amp=True: amp_coef >= 0, phase in (\([-pi, pi]\) (via atan2). If pos_amp=False: allow signed amp_coef and shift phase by +/-pi so that phase is constrained to \([-pi/2, pi/2]\) (minimize |phase|).

 False
neg_amp bool

Force amplitude to be negative.

 False
take_linear_sign bool

Only makes a difference if to_time=True.
If True, will make coef_linear_out the same sign as coef_linear.

 True

Returns:

Name Type Description
coef_linear_out Union[float, ndarray, DataArray]

Output amplitude coefficient, same container type as inputs. Interpreted as \(\lambda_{\mathrm{mod}}\) if to_time=True; otherwise \(\lambda\).
coef_phase_out Union[float, ndarray, DataArray]

Output phase quantity, same container type as inputs. Interpreted as \(2\pi f t_{\mathrm{ph}}\) (radians) if to_time=True; otherwise \(\lambda_{\mathrm{ph}}\).
Source code in isca_tools/thesis/surface_energy_budget.py 
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
def phase_coef_conversion(coef_linear: Union[float, np.ndarray, xr.DataArray],
                          coef_phase: Union[float, np.ndarray, xr.DataArray],
                          to_time: bool = True,
                          pos_amp: bool = False,
                          neg_amp: bool = False,
                          take_linear_sign: bool = True
                          ) -> Tuple[Union[float, np.ndarray, xr.DataArray], Union[float, np.ndarray, xr.DataArray]]:
    """
    Method to convert between interpretation of phase empirical coefficients. If `to_time=True`, expect
    input to be $\lambda$ and $\lambda_{ph}$, and will return
    $\lambda_{mod}$, and $2\pi ft_{ph}$. If `to_time=False`, expect the opposite.

    Conversion performed according to:

    \\begin{align}
        \lambda &= \lambda_{\\text{mod}}\cos(2 \pi ft_{\\text{ph}})
        \\\\
        \lambda_{\\text{ph}} &= \lambda_{\\text{mod}}\sin(2 \pi ft_{\\text{ph}})
    \end{align}

    Args:
        coef_linear:
            Cosine/linear-amplitude coefficient associated with the in-phase component.
            If ``to_time=True`` this is interpreted as \(\lambda\); if ``to_time=False`` this is
            interpreted as \(\lambda_{\mathrm{mod}}\).
            Can be a scalar, ``numpy.ndarray``, or ``xarray.DataArray``.
        coef_phase:
            Phase-related coefficient.
            If ``to_time=True`` this is interpreted as \(\lambda_{\mathrm{ph}}\) (sine coefficient);
            if ``to_time=False`` this is interpreted as the phase angle \(2\pi f t_{\mathrm{ph}}\)
            (radians) used inside \(\cos(\cdot)\) / \(\sin(\cdot)\).
            Can be a scalar, ``numpy.ndarray``, or ``xarray.DataArray``.
        to_time:
            Direction of conversion.
            If ``True``, convert \((\lambda,\lambda_{\mathrm{ph}})\\rightarrow
            (\lambda_{\mathrm{mod}},2\pi f t_{\mathrm{ph}})\).
            If ``False``, convert \((\lambda_{\mathrm{mod}},2\pi f t_{\mathrm{ph}})
            \\rightarrow (\lambda,\lambda_{\mathrm{ph}})\).
        pos_amp:
            Only makes a difference if `to_time=True`.</br>
            If pos_amp=True: amp_coef >= 0, phase in ($[-pi, pi]$ (via atan2).
            If pos_amp=False: allow signed amp_coef and shift phase by +/-pi so that phase
            is constrained to $[-pi/2, pi/2]$ (minimize |phase|).
        neg_amp: Force amplitude to be negative.
        take_linear_sign:
            Only makes a difference if `to_time=True`.</br>
            If `True`, will make coef_linear_out the same sign as coef_linear.

    Returns:
        coef_linear_out:
            Output amplitude coefficient, same container type as inputs.
            Interpreted as \(\lambda_{\mathrm{mod}}\) if ``to_time=True``; otherwise \(\lambda\).
        coef_phase_out:
            Output phase quantity, same container type as inputs.
            Interpreted as \(2\pi f t_{\mathrm{ph}}\) (radians) if ``to_time=True``; otherwise
            \(\lambda_{\mathrm{ph}}\).
    """
    if to_time:
        coef_linear_out, coef_phase_out = fourier.coef_conversion(sin_coef=coef_phase, cos_coef=coef_linear,
                                                                  pos_amp=pos_amp, take_cos_sign=take_linear_sign,
                                                                  neg_amp=neg_amp)
        # if freq is not None:
        #     coef_phase_out = coef_phase_out / (2 * np.pi * freq)  # convert to seconds
    else:
        coef_linear_out = coef_linear * np.cos(coef_phase)
        coef_phase_out = coef_linear * np.sign(coef_phase)
    return coef_linear_out, coef_phase_out

swdn_from_temp_fourier(time, temp_fourier_amp, temp_fourier_phase, heat_capacity, lambda_const, lambda_time_lag=None, lambda_nl=None, day_seconds=86400, single_harmonic_nl=False)

OUTDATED FUNCTION - USED FOR get_temp_fourier but now use get_temp_fourier_numerical

This inverts the linearized surface energy budget to return an approximation for downward shortwave radiation at the surface, \(F(t)\), given a Fourier approximation for surface temperature, \(T(t) = \frac{T_0}{2} + \sum_{n=1}^{N} T_n\cos(2n\pi ft - \phi_n)\):

\[ F(t) \approx C\frac{\partial T}{\partial t} + \lambda_0 + \sum_{i=1}^{N_{\lambda}}\lambda_i T(t-\Lambda_i) + \sum_{j=2}^{j_{max}}\lambda_{nl_j}\left(T'^{j}(t) - \overline{T'^j}\right) \]

Parameters:

Name Type Description Default
time ndarray

float [n_time]
Time in days (assumes periodic e.g. annual mean, so time = np.arange(360)).

 required
temp_fourier_amp ndarray

float [n_harmonics+1]
The amplitude Fourier coefficients for surface temperature: \(T_n\).

 required
temp_fourier_phase ndarray

float [n_harmonics]
The phase Fourier coefficients for surface temperature: \(\phi_n\).

 required
heat_capacity float

\(C\), the heat capacity of the surface in units of \(JK^{-1}m^{-2}\).
Obtained from mixed layer depth of ocean using get_heat_capacity.

 required
lambda_const ndarray

float [n_lambda+1]
The constants \(\lambda_i\) used in the approximation for \(\Gamma^{\uparrow} = LW^{\uparrow} - LW^{\downarrow} + LH^{\uparrow} + SH^{\uparrow}\).
lambda_const[0] is \(\lambda_0\) and lambda_const[i] is \(\lambda_{i}\) for \(i>0\).

 required
lambda_time_lag Optional[ndarray]

float [n_lambda]
The constants \(\Lambda_i\) used in the approximation for \(\Gamma^{\uparrow}\).
lambda_time_lag[0] is \(\Lambda_1\) and lambda_time_lag[i] is \(\Lambda_{i+1}\) for \(i>0\).

 None
lambda_nl Optional[Union[float, ndarray]]

float [n_lambda_nl] The constants \(\lambda_{nl}\) used in the approximation for \(\Gamma^{\uparrow}\). [0] is squared contribution, [1] is cubed, ...

 None
day_seconds float

Duration of a day in seconds.

 86400
single_harmonic_nl bool

If True, the \(\lambda_{nl_j}T'^j\) terms in \(\Gamma^{\uparrow}\) will only use the first harmonic, not all harmonics.

 False

Returns:

Type Description
ndarray

float [n_time]
Approximation for downward shortwave radiation at the surface. Units: \(Wm^{-2}\).
Source code in isca_tools/thesis/surface_energy_budget.py 
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
def swdn_from_temp_fourier(time: np.ndarray, temp_fourier_amp: np.ndarray, temp_fourier_phase: np.ndarray,
                           heat_capacity: float, lambda_const: np.ndarray,
                           lambda_time_lag: Optional[np.ndarray] = None,
                           lambda_nl: Optional[Union[float, np.ndarray]] = None,
                           day_seconds: float = 86400, single_harmonic_nl: bool = False) -> np.ndarray:
    """
    OUTDATED FUNCTION - USED FOR `get_temp_fourier` but now use `get_temp_fourier_numerical`

    This inverts the linearized surface energy budget to return an approximation for downward shortwave radiation
    at the surface, $F(t)$, given a Fourier approximation for surface temperature,
    $T(t) = \\frac{T_0}{2} + \\sum_{n=1}^{N} T_n\\cos(2n\\pi ft - \\phi_n)$:

    $$
    F(t) \\approx C\\frac{\partial T}{\partial t} + \lambda_0 + \sum_{i=1}^{N_{\lambda}}\lambda_i T(t-\Lambda_i) +
    \sum_{j=2}^{j_{max}}\lambda_{nl_j}\\left(T'^{j}(t) - \overline{T'^j}\\right)
    $$


    Args:
        time: `float [n_time]`</br>
            Time in days (assumes periodic e.g. annual mean, so `time = np.arange(360)`).
        temp_fourier_amp: `float [n_harmonics+1]`</br>
            The amplitude Fourier coefficients for surface temperature: $T_n$.
        temp_fourier_phase: `float [n_harmonics]`</br>
            The phase Fourier coefficients for surface temperature: $\phi_n$.
        heat_capacity: $C$, the heat capacity of the surface in units of $JK^{-1}m^{-2}$.</br>
            Obtained from mixed layer depth of ocean using
            [`get_heat_capacity`](../utils/radiation.md#isca_tools.utils.radiation.get_heat_capacity).
        lambda_const: `float [n_lambda+1]`</br>
            The constants $\lambda_i$ used in the approximation for
            $\Gamma^{\\uparrow} = LW^{\\uparrow} - LW^{\\downarrow} + LH^{\\uparrow} + SH^{\\uparrow}$.</br>
            `lambda_const[0]` is $\lambda_0$ and `lambda_const[i]` is $\lambda_{i}$ for $i>0$.
        lambda_time_lag: `float [n_lambda]`</br>
            The constants $\Lambda_i$ used in the approximation for $\Gamma^{\\uparrow}$.</br>
            `lambda_time_lag[0]` is $\Lambda_1$ and `lambda_time_lag[i]` is $\Lambda_{i+1}$ for $i>0$.
        lambda_nl: `float [n_lambda_nl]`
            The constants $\lambda_{nl}$ used in the approximation for $\Gamma^{\\uparrow}$.
            `[0]` is squared contribution, `[1]` is cubed, ...
        day_seconds: Duration of a day in seconds.
        single_harmonic_nl: If `True`, the $\lambda_{nl_j}T'^j$ terms in $\Gamma^{\\uparrow}$ will only
            use the first harmonic, not all harmonics.

    Returns:
        `float [n_time]`</br>
            Approximation for downward shortwave radiation at the surface.
            Units: $Wm^{-2}$.
    """
    if lambda_nl is not None:
        #  deal with case when float given as lambda_nl
        if not hasattr(lambda_nl, "__len__"):
            lambda_nl = np.asarray([lambda_nl])
    n_year_days = len(time)
    temp_fourier = fourier.fourier_series(time, temp_fourier_amp, temp_fourier_phase)
    if single_harmonic_nl:
        temp_anom_nl = fourier.fourier_series(time, [0, temp_fourier_amp[1]],
                                              [temp_fourier_phase[0]])
        gamma = gamma_linear_approx(time, temp_fourier, lambda_const,
                                    lambda_time_lag, lambda_nl, temp_anom_nl)
    else:
        gamma = gamma_linear_approx(time, temp_fourier, lambda_const,
                                    lambda_time_lag, lambda_nl)
    dtemp_dt = fourier.fourier_series_deriv(time, temp_fourier_amp, temp_fourier_phase, day_seconds)
    return heat_capacity * dtemp_dt + gamma