def CoupledChannel_time(nt, nx, xx, dt, avep, sst, tas, hice, sst_boundary, sst_out, tas_out, hice_out, sflx_f_out, sflx_out, forcing, spatial_pattern, ra, Cp, va, vo, Da, Do, Cs, T0, Tf, emissivity, SW0, SW_anom, H, Hb, Cpo, ro, tau_entrainment, Li, ri, use_ocn_tendencies, use_atm_tendencies, atm_DA_tendencies, ocn_DA_tendencies, ice_model, atm_adv, return_coupled_fluxes):
    '\n    Separate time loop to enable numba\n    '
    c = 0
    c2 = 0
    c3 = 0
    n = 1
    for nn in range(nt):
        sflx = (forcing[nn] * spatial_pattern)
        sflx_f_out[c, :] = (sflx_f_out[c, :] + sflx)
        sflx = (sflx + ((((ra * Cp) * va) * Cs) * (sst[(n - 1), :] - tas[(n - 1), :])))
        LW_cooling = ((emissivity * 5.67e-08) * (tas[(n - 1), :] ** 4))
        SW_warming = (SW0 + max((SW_anom * np.sin(((((2 * float(nn)) * dt) * np.pi) / ((360 * 24) * 3600)))), 0.0))
        net_radiation = (- LW_cooling)
        sst_boundary[n] = (((sst_boundary[(n - 1)] + ((SW_warming[0] * dt) / ((H * Cpo) * ro))) - ((((emissivity * 5.67e-08) * (sst_boundary[(n - 1)] ** 4)) * dt) / ((H * Cpo) * ro))) + (((T0 - sst_boundary[(n - 1)]) * dt) / ((360 * 24) * 3600)))
        if atm_adv:
            a_adv = (np.concatenate(((sst_boundary[(n - 1)] - tas[(n - 1), :1]), (tas[(n - 1), :(- 1)] - tas[(n - 1), 1:])), axis=0) * ((va * dt) / xx))
        else:
            a_adv = np.zeros(nx)
        a_diff = (((tas[(n - 1), 2:] + tas[(n - 1), :(- 2)]) - (2 * tas[(n - 1), 1:(- 1)])) * ((Da * dt) / (xx ** 2)))
        a_diff0 = (((tas[((n - 1), 1)] + sst_boundary[(n - 1)]) - (2 * tas[((n - 1), 0)])) * ((Da * dt) / (xx ** 2)))
        a_diff = np.concatenate((np.array([a_diff0]), a_diff, a_diff[(- 1):]), axis=0)
        a_netsflx = (((sflx * dt) / ((Hb * Cp) * ra)) + ((net_radiation * dt) / ((Hb * Cp) * ra)))
        if return_coupled_fluxes:
            atm_DA_tendencies[nn, :] = np.sum((a_adv, a_diff), axis=0)
        if use_atm_tendencies:
            tas[n, :] = ((tas[(n - 1), :] + a_netsflx) + atm_DA_tendencies[c3, :])
        else:
            tas[n, :] = (((tas[(n - 1), :] + a_netsflx) + a_adv) + a_diff)
        o_adv = (np.concatenate(((sst_boundary[(n - 1)] - sst[(n - 1), :1]), (sst[(n - 1), :(- 1)] - sst[(n - 1), 1:])), axis=0) * ((vo * dt) / xx))
        o_diff = (((sst[(n - 1), 2:] + sst[(n - 1), :(- 2)]) - (2 * sst[(n - 1), 1:(- 1)])) * ((Do * dt) / (xx ** 2)))
        o_diff0 = (((sst[((n - 1), 1)] + sst_boundary[(n - 1)]) - (2 * sst[((n - 1), 0)])) * ((Do * dt) / (xx ** 2)))
        o_diff = np.concatenate((np.array([o_diff0]), o_diff, o_diff[(- 1):]), axis=0)
        o_entrain = (((T0 - sst[(n - 1), :]) * dt) / tau_entrainment)
        o_netsflx = ((((- sflx) * dt) / ((H * Cpo) * ro)) + ((SW_warming * dt) / ((H * Cpo) * ro)))
        if return_coupled_fluxes:
            ocn_DA_tendencies[nn, :] = ((o_adv + o_diff) + o_entrain)
        if use_ocn_tendencies:
            sst[n, :] = ((sst[(n - 1), :] + o_netsflx) + ocn_DA_tendencies[c3, :])
        else:
            sst[n, :] = ((((sst[(n - 1), :] + o_netsflx) + o_adv) + o_diff) + o_entrain)
        if ice_model:
            ice_mask = (hice[(n - 1), :] > 0).astype(np.float)
            freezing_mask = (sst[n, :] < Tf).astype(np.float)
            dEdt = ((((H * ro) * Cpo) * (sst[n, :] - sst[(n - 1), :])) / dt)
            excess_freeze = (freezing_mask * np.max([(- dEdt), np.zeros(nx)], axis=0))
            excess_melt = (ice_mask * np.max([dEdt, np.zeros(nx)], axis=0))
            dhice_freeze = ((dt * excess_freeze) / (Li * ri))
            dhice_melt = ((dt * excess_melt) / (Li * ri))
            hice[n, :] = ((hice[(n - 1), :] + dhice_freeze) - dhice_melt)
            hice_melt = ((dhice_melt > 0).astype(np.float) * np.min([dhice_melt, hice[(n - 1), :]], axis=0))
            hice[n, :] = np.max([hice[n, :], np.zeros(nx)], axis=0)
            sst[n, :] = (sst[n, :] + ((dt * excess_freeze) / ((H * Cpo) * ro)))
            sst[n, :] = (sst[n, :] - ((hice_melt * (Li * ri)) / ((ro * Cpo) * H)))
        tas_out[c, :] = (tas_out[c, :] + tas[n, :])
        sst_out[c, :] = (sst_out[c, :] + sst[n, :])
        hice_out[c, :] = (hice_out[c, :] + hice[n, :])
        sflx_out[c, :] = (sflx_out[c, :] + sflx)
        c2 = (c2 + 1)
        c3 = (c3 + 1)
        if ((((nn + 1) * dt) % ((360 * 24) * 3600)) == 0):
            c3 = 0
        if (((((nn + 1) * dt) % avep) == 0) and (nn > 0)):
            tas_out[c, :] = (tas_out[c, :] / c2)
            sst_out[c, :] = (sst_out[c, :] / c2)
            sflx_out[c, :] = (sflx_out[c, :] / c2)
            sflx_f_out[c, :] = (sflx_f_out[c, :] / c2)
            hice_out[c, :] = (hice_out[c, :] / c2)
            c = (c + 1)
            c2 = 0
        tas[0, :] = tas[1, :].copy()
        sst[0, :] = sst[1, :].copy()
        hice[0, :] = hice[1, :].copy()
        sst_boundary[0] = sst_boundary[1].copy()
    hice_out[np.where((hice_out == 0))] = np.nan
    return (tas_out, sst_out, hice_out, sflx_out, sflx_f_out, atm_DA_tendencies, ocn_DA_tendencies)
Separate time loop to enable numba
def CoupledChannel_time(nt, nx, xx, dt, avep, sst, tas, hice, sst_boundary, sst_out, tas_out, hice_out, sflx_f_out, sflx_out, forcing, spatial_pattern, ra, Cp, va, vo, Da, Do, Cs, T0, Tf, emissivity, SW0, SW_anom, H, Hb, Cpo, ro, tau_entrainment, Li, ri, use_ocn_tendencies, use_atm_tendencies, atm_DA_tendencies, ocn_DA_tendencies, ice_model, atm_adv, return_coupled_fluxes):
    '\n    \n    '
    c = 0
    c2 = 0
    c3 = 0
    n = 1
    for nn in range(nt):
        sflx = (forcing[nn] * spatial_pattern)
        sflx_f_out[c, :] = (sflx_f_out[c, :] + sflx)
        sflx = (sflx + ((((ra * Cp) * va) * Cs) * (sst[(n - 1), :] - tas[(n - 1), :])))
        LW_cooling = ((emissivity * 5.67e-08) * (tas[(n - 1), :] ** 4))
        SW_warming = (SW0 + max((SW_anom * np.sin(((((2 * float(nn)) * dt) * np.pi) / ((360 * 24) * 3600)))), 0.0))
        net_radiation = (- LW_cooling)
        sst_boundary[n] = (((sst_boundary[(n - 1)] + ((SW_warming[0] * dt) / ((H * Cpo) * ro))) - ((((emissivity * 5.67e-08) * (sst_boundary[(n - 1)] ** 4)) * dt) / ((H * Cpo) * ro))) + (((T0 - sst_boundary[(n - 1)]) * dt) / ((360 * 24) * 3600)))
        if atm_adv:
            a_adv = (np.concatenate(((sst_boundary[(n - 1)] - tas[(n - 1), :1]), (tas[(n - 1), :(- 1)] - tas[(n - 1), 1:])), axis=0) * ((va * dt) / xx))
        else:
            a_adv = np.zeros(nx)
        a_diff = (((tas[(n - 1), 2:] + tas[(n - 1), :(- 2)]) - (2 * tas[(n - 1), 1:(- 1)])) * ((Da * dt) / (xx ** 2)))
        a_diff0 = (((tas[((n - 1), 1)] + sst_boundary[(n - 1)]) - (2 * tas[((n - 1), 0)])) * ((Da * dt) / (xx ** 2)))
        a_diff = np.concatenate((np.array([a_diff0]), a_diff, a_diff[(- 1):]), axis=0)
        a_netsflx = (((sflx * dt) / ((Hb * Cp) * ra)) + ((net_radiation * dt) / ((Hb * Cp) * ra)))
        if return_coupled_fluxes:
            atm_DA_tendencies[nn, :] = np.sum((a_adv, a_diff), axis=0)
        if use_atm_tendencies:
            tas[n, :] = ((tas[(n - 1), :] + a_netsflx) + atm_DA_tendencies[c3, :])
        else:
            tas[n, :] = (((tas[(n - 1), :] + a_netsflx) + a_adv) + a_diff)
        o_adv = (np.concatenate(((sst_boundary[(n - 1)] - sst[(n - 1), :1]), (sst[(n - 1), :(- 1)] - sst[(n - 1), 1:])), axis=0) * ((vo * dt) / xx))
        o_diff = (((sst[(n - 1), 2:] + sst[(n - 1), :(- 2)]) - (2 * sst[(n - 1), 1:(- 1)])) * ((Do * dt) / (xx ** 2)))
        o_diff0 = (((sst[((n - 1), 1)] + sst_boundary[(n - 1)]) - (2 * sst[((n - 1), 0)])) * ((Do * dt) / (xx ** 2)))
        o_diff = np.concatenate((np.array([o_diff0]), o_diff, o_diff[(- 1):]), axis=0)
        o_entrain = (((T0 - sst[(n - 1), :]) * dt) / tau_entrainment)
        o_netsflx = ((((- sflx) * dt) / ((H * Cpo) * ro)) + ((SW_warming * dt) / ((H * Cpo) * ro)))
        if return_coupled_fluxes:
            ocn_DA_tendencies[nn, :] = ((o_adv + o_diff) + o_entrain)
        if use_ocn_tendencies:
            sst[n, :] = ((sst[(n - 1), :] + o_netsflx) + ocn_DA_tendencies[c3, :])
        else:
            sst[n, :] = ((((sst[(n - 1), :] + o_netsflx) + o_adv) + o_diff) + o_entrain)
        if ice_model:
            ice_mask = (hice[(n - 1), :] > 0).astype(np.float)
            freezing_mask = (sst[n, :] < Tf).astype(np.float)
            dEdt = ((((H * ro) * Cpo) * (sst[n, :] - sst[(n - 1), :])) / dt)
            excess_freeze = (freezing_mask * np.max([(- dEdt), np.zeros(nx)], axis=0))
            excess_melt = (ice_mask * np.max([dEdt, np.zeros(nx)], axis=0))
            dhice_freeze = ((dt * excess_freeze) / (Li * ri))
            dhice_melt = ((dt * excess_melt) / (Li * ri))
            hice[n, :] = ((hice[(n - 1), :] + dhice_freeze) - dhice_melt)
            hice_melt = ((dhice_melt > 0).astype(np.float) * np.min([dhice_melt, hice[(n - 1), :]], axis=0))
            hice[n, :] = np.max([hice[n, :], np.zeros(nx)], axis=0)
            sst[n, :] = (sst[n, :] + ((dt * excess_freeze) / ((H * Cpo) * ro)))
            sst[n, :] = (sst[n, :] - ((hice_melt * (Li * ri)) / ((ro * Cpo) * H)))
        tas_out[c, :] = (tas_out[c, :] + tas[n, :])
        sst_out[c, :] = (sst_out[c, :] + sst[n, :])
        hice_out[c, :] = (hice_out[c, :] + hice[n, :])
        sflx_out[c, :] = (sflx_out[c, :] + sflx)
        c2 = (c2 + 1)
        c3 = (c3 + 1)
        if ((((nn + 1) * dt) % ((360 * 24) * 3600)) == 0):
            c3 = 0
        if (((((nn + 1) * dt) % avep) == 0) and (nn > 0)):
            tas_out[c, :] = (tas_out[c, :] / c2)
            sst_out[c, :] = (sst_out[c, :] / c2)
            sflx_out[c, :] = (sflx_out[c, :] / c2)
            sflx_f_out[c, :] = (sflx_f_out[c, :] / c2)
            hice_out[c, :] = (hice_out[c, :] / c2)
            c = (c + 1)
            c2 = 0
        tas[0, :] = tas[1, :].copy()
        sst[0, :] = sst[1, :].copy()
        hice[0, :] = hice[1, :].copy()
        sst_boundary[0] = sst_boundary[1].copy()
    hice_out[np.where((hice_out == 0))] = np.nan
    return (tas_out, sst_out, hice_out, sflx_out, sflx_f_out, atm_DA_tendencies, ocn_DA_tendencies)