def CoupledChannel(C, forcing, T_boundary=None, dt_f=((30 * 24) * 3600), restoring=False, ice_model=True, atm_adv=True, spatial_pattern=None, atm_DA_tendencies=None, ocn_DA_tendencies=None, return_coupled_fluxes=False, random_amp=0.1):
    '\n    This is the main function for the coupled ocean--atm channel model.\n    \n    ## INPUT VARIABLES ##\n    \n    tmax: running time in seconds\n    avep: averaging period for the ouput\n    T0: initial temperature\n    forcing: dimensionless scaling for the heat flux forcing - default strength is 5 W/m2\n    dt_f: timestep of the forcing\n    atm_adv: boolean, advective atmosphere\n    atm_ocn: boolean, advective ocean\n    '
    nt = int((C['tmax'] / C['dt']))
    nt1 = int((C['tmax'] / C['avep']))
    sst = (C['T0'] * np.ones((2, C['nx'])))
    tas = (C['T0'] * np.ones((2, C['nx'])))
    hice = np.zeros((2, C['nx']))
    SW0 = np.tile(C['SW0'][:, np.newaxis], (1, nt))
    naxis = np.tile(np.arange(nt)[(np.newaxis,)], (C['nx'], 1))
    SW_warming = np.max(np.concatenate([(SW0 - (C['SW_anom'] * np.cos((((2 * np.pi) * (naxis * C['dt'])) / ((360 * 24) * 3600)))))[(np.newaxis,)], np.zeros((C['nx'], nt))[(np.newaxis,)]], axis=0), 0)
    if np.all((T_boundary == None)):
        T_boundary = (C['T0'] * np.ones(nt))
    sst_boundary = (T_boundary[0] * np.ones(2))
    if np.all((forcing != None)):
        forcing = np.interp(np.arange(0, (len(forcing) * dt_f), C['dt']), np.arange(0, (len(forcing) * dt_f), dt_f), forcing)
    else:
        forcing = np.zeros((nt + 1))
    sst_out = np.zeros((nt1, C['nx']))
    tas_out = np.zeros((nt1, C['nx']))
    hice_out = np.zeros((nt1, C['nx']))
    sflx_f_out = np.zeros((nt1, C['nx']))
    sflx_out = np.zeros((nt1, C['nx']))
    if np.all((spatial_pattern == None)):
        spatial_pattern = np.ones(C['nx'])
    if np.all((atm_DA_tendencies != None)):
        use_atm_tendencies = True
    else:
        use_atm_tendencies = False
    if np.all((ocn_DA_tendencies != None)):
        use_ocn_tendencies = True
    else:
        use_ocn_tendencies = False
    if return_coupled_fluxes:
        atm_DA_tendencies = np.zeros((nt, C['nx']))
        ocn_DA_tendencies = np.zeros((nt, C['nx']))
    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 + ((((C['ra'] * C['Cp']) * C['va']) * C['Cs']) * (sst[(n - 1), :] - tas[(n - 1), :])))
        LW_cooling = ((C['emissivity'] * 5.67e-08) * (tas[(n - 1), :] ** 4))
        sst_boundary_tendency = ((((SW_warming[(0, nn)] * C['dt']) / ((C['H'] * C['Cpo']) * C['ro'])) - ((((C['emissivity'] * 5.67e-08) * (sst_boundary[(n - 1)] ** 4)) * C['dt']) / ((C['H'] * C['Cpo']) * C['ro']))) + (((T_boundary[nn] - sst_boundary[(n - 1)]) * C['dt']) / C['period']))
        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) * ((C['va'] * C['dt']) / C['xx']))
        else:
            a_adv = 0
        a_diff = (((tas[(n - 1), 2:] + tas[(n - 1), :(- 2)]) - (2 * tas[(n - 1), 1:(- 1)])) * ((C['Da'] * C['dt']) / (C['xx'] ** 2)))
        a_diff0 = (((tas[((n - 1), 1)] + sst_boundary[(n - 1)]) - (2 * tas[((n - 1), 0)])) * ((C['Da'] * C['dt']) / (C['xx'] ** 2)))
        a_diff = np.concatenate([np.array([a_diff0]), a_diff, a_diff[(- 1):]], axis=0)
        a_netsflx = (((sflx * C['dt']) / ((C['Hb'] * C['Cp']) * C['ra'])) - ((LW_cooling * C['dt']) / ((C['Hb'] * C['Cp']) * C['ra'])))
        if return_coupled_fluxes:
            atm_DA_tendencies[nn, :] = (a_adv + a_diff)
        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) * ((C['vo'] * C['dt']) / C['xx']))
        o_diff = (((sst[(n - 1), 2:] + sst[(n - 1), :(- 2)]) - (2 * sst[(n - 1), 1:(- 1)])) * ((C['Do'] * C['dt']) / (C['xx'] ** 2)))
        o_diff0 = (((sst[((n - 1), 1)] + sst_boundary[(n - 1)]) - (2 * sst[((n - 1), 0)])) * ((C['Do'] * C['dt']) / (C['xx'] ** 2)))
        o_diff = np.concatenate([np.array([o_diff0]), o_diff, o_diff[(- 1):]], axis=0)
        if (C['tau_entrainment'] > 0):
            o_entrain = (((C['T0'] - sst[(n - 1), :]) * C['dt']) / C['tau_entrainment'])
        else:
            o_entrain = 0
        o_netsflx = ((((- sflx) * C['dt']) / ((C['H'] * C['Cpo']) * C['ro'])) + ((SW_warming[:, nn] * C['dt']) / ((C['H'] * C['Cpo']) * C['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, :] < C['Tf']).astype(np.float)
            dEdt = ((((C['H'] * C['ro']) * C['Cpo']) * (sst[n, :] - sst[(n - 1), :])) / C['dt'])
            excess_freeze = (freezing_mask * np.max([(- dEdt), np.zeros(C['nx'])], axis=0))
            excess_melt = (ice_mask * np.max([dEdt, np.zeros(C['nx'])], axis=0))
            dhice_freeze = ((C['dt'] * excess_freeze) / (C['Li'] * C['ri']))
            dhice_melt = ((C['dt'] * excess_melt) / (C['Li'] * C['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(C['nx'])], axis=0)
            sst[n, :] = (sst[n, :] + ((C['dt'] * excess_freeze) / ((C['H'] * C['Cpo']) * C['ro'])))
            sst[n, :] = (sst[n, :] - ((hice_melt * (C['Li'] * C['ri'])) / ((C['ro'] * C['Cpo']) * C['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) * C['dt']) % ((360 * 24) * 3600)) == 0):
            c3 = 0
        if (((((nn + 1) * C['dt']) % C['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
            if ((((nn + 1) * C['dt']) % ((360 * 24) * 3600)) == 0):
                print('Year ', (((nn + 1) * C['dt']) / ((360 * 24) * 3600)), sst[(1, int((C['nx'] / 4)))], sst[(1, int(((3 * C['nx']) / 4)))])
        tas[0, :] = tas[1, :].copy()
        sst[0, :] = sst[1, :].copy()
        hice[0, :] = hice[1, :].copy()
        sst_boundary[(n - 1)] = (sst_boundary[(n - 1)] + sst_boundary_tendency)
    hice_out[np.where((hice_out == 0))] = np.nan
    if return_coupled_fluxes:
        return (tas_out, sst_out, hice_out, sflx_out, sflx_f_out, nt1, nt, atm_DA_tendencies, ocn_DA_tendencies)
    else:
        return (tas_out, sst_out, hice_out, sflx_out, sflx_f_out, nt1, nt)
This is the main function for the coupled ocean--atm channel model.

## INPUT VARIABLES ##

tmax: running time in seconds
avep: averaging period for the ouput
T0: initial temperature
forcing: dimensionless scaling for the heat flux forcing - default strength is 5 W/m2
dt_f: timestep of the forcing
atm_adv: boolean, advective atmosphere
atm_ocn: boolean, advective ocean
def CoupledChannel(C, forcing, T_boundary=None, dt_f=((30 * 24) * 3600), restoring=False, ice_model=True, atm_adv=True, spatial_pattern=None, atm_DA_tendencies=None, ocn_DA_tendencies=None, return_coupled_fluxes=False, random_amp=0.1):
    '\n    This is the main function for the coupled ocean--atm channel model.\n    \n    ## INPUT VARIABLES ##\n    \n    tmax: running time in seconds\n    avep: averaging period for the ouput\n    T0: initial temperature\n    forcing: dimensionless scaling for the heat flux forcing - default strength is 5 W/m2\n    dt_f: timestep of the forcing\n    atm_adv: boolean, advective atmosphere\n    atm_ocn: boolean, advective ocean\n    '
    nt = int((C['tmax'] / C['dt']))
    nt1 = int((C['tmax'] / C['avep']))
    sst = (C['T0'] * np.ones((2, C['nx'])))
    tas = (C['T0'] * np.ones((2, C['nx'])))
    hice = np.zeros((2, C['nx']))
    SW0 = np.tile(C['SW0'][:, np.newaxis], (1, nt))
    naxis = np.tile(np.arange(nt)[(np.newaxis,)], (C['nx'], 1))
    SW_warming = np.max(np.concatenate([(SW0 - (C['SW_anom'] * np.cos((((2 * np.pi) * (naxis * C['dt'])) / ((360 * 24) * 3600)))))[(np.newaxis,)], np.zeros((C['nx'], nt))[(np.newaxis,)]], axis=0), 0)
    if np.all((T_boundary == None)):
        T_boundary = (C['T0'] * np.ones(nt))
    sst_boundary = (T_boundary[0] * np.ones(2))
    if np.all((forcing != None)):
        forcing = np.interp(np.arange(0, (len(forcing) * dt_f), C['dt']), np.arange(0, (len(forcing) * dt_f), dt_f), forcing)
    else:
        forcing = np.zeros((nt + 1))
    sst_out = np.zeros((nt1, C['nx']))
    tas_out = np.zeros((nt1, C['nx']))
    hice_out = np.zeros((nt1, C['nx']))
    sflx_f_out = np.zeros((nt1, C['nx']))
    sflx_out = np.zeros((nt1, C['nx']))
    if np.all((spatial_pattern == None)):
        spatial_pattern = np.ones(C['nx'])
    if np.all((atm_DA_tendencies != None)):
        use_atm_tendencies = True
    else:
        use_atm_tendencies = False
    if np.all((ocn_DA_tendencies != None)):
        use_ocn_tendencies = True
    else:
        use_ocn_tendencies = False
    if return_coupled_fluxes:
        atm_DA_tendencies = np.zeros((nt, C['nx']))
        ocn_DA_tendencies = np.zeros((nt, C['nx']))
    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 + ((((C['ra'] * C['Cp']) * C['va']) * C['Cs']) * (sst[(n - 1), :] - tas[(n - 1), :])))
        LW_cooling = ((C['emissivity'] * 5.67e-08) * (tas[(n - 1), :] ** 4))
        sst_boundary_tendency = ((((SW_warming[(0, nn)] * C['dt']) / ((C['H'] * C['Cpo']) * C['ro'])) - ((((C['emissivity'] * 5.67e-08) * (sst_boundary[(n - 1)] ** 4)) * C['dt']) / ((C['H'] * C['Cpo']) * C['ro']))) + (((T_boundary[nn] - sst_boundary[(n - 1)]) * C['dt']) / C['period']))
        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) * ((C['va'] * C['dt']) / C['xx']))
        else:
            a_adv = 0
        a_diff = (((tas[(n - 1), 2:] + tas[(n - 1), :(- 2)]) - (2 * tas[(n - 1), 1:(- 1)])) * ((C['Da'] * C['dt']) / (C['xx'] ** 2)))
        a_diff0 = (((tas[((n - 1), 1)] + sst_boundary[(n - 1)]) - (2 * tas[((n - 1), 0)])) * ((C['Da'] * C['dt']) / (C['xx'] ** 2)))
        a_diff = np.concatenate([np.array([a_diff0]), a_diff, a_diff[(- 1):]], axis=0)
        a_netsflx = (((sflx * C['dt']) / ((C['Hb'] * C['Cp']) * C['ra'])) - ((LW_cooling * C['dt']) / ((C['Hb'] * C['Cp']) * C['ra'])))
        if return_coupled_fluxes:
            atm_DA_tendencies[nn, :] = (a_adv + a_diff)
        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) * ((C['vo'] * C['dt']) / C['xx']))
        o_diff = (((sst[(n - 1), 2:] + sst[(n - 1), :(- 2)]) - (2 * sst[(n - 1), 1:(- 1)])) * ((C['Do'] * C['dt']) / (C['xx'] ** 2)))
        o_diff0 = (((sst[((n - 1), 1)] + sst_boundary[(n - 1)]) - (2 * sst[((n - 1), 0)])) * ((C['Do'] * C['dt']) / (C['xx'] ** 2)))
        o_diff = np.concatenate([np.array([o_diff0]), o_diff, o_diff[(- 1):]], axis=0)
        if (C['tau_entrainment'] > 0):
            o_entrain = (((C['T0'] - sst[(n - 1), :]) * C['dt']) / C['tau_entrainment'])
        else:
            o_entrain = 0
        o_netsflx = ((((- sflx) * C['dt']) / ((C['H'] * C['Cpo']) * C['ro'])) + ((SW_warming[:, nn] * C['dt']) / ((C['H'] * C['Cpo']) * C['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, :] < C['Tf']).astype(np.float)
            dEdt = ((((C['H'] * C['ro']) * C['Cpo']) * (sst[n, :] - sst[(n - 1), :])) / C['dt'])
            excess_freeze = (freezing_mask * np.max([(- dEdt), np.zeros(C['nx'])], axis=0))
            excess_melt = (ice_mask * np.max([dEdt, np.zeros(C['nx'])], axis=0))
            dhice_freeze = ((C['dt'] * excess_freeze) / (C['Li'] * C['ri']))
            dhice_melt = ((C['dt'] * excess_melt) / (C['Li'] * C['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(C['nx'])], axis=0)
            sst[n, :] = (sst[n, :] + ((C['dt'] * excess_freeze) / ((C['H'] * C['Cpo']) * C['ro'])))
            sst[n, :] = (sst[n, :] - ((hice_melt * (C['Li'] * C['ri'])) / ((C['ro'] * C['Cpo']) * C['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) * C['dt']) % ((360 * 24) * 3600)) == 0):
            c3 = 0
        if (((((nn + 1) * C['dt']) % C['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
            if ((((nn + 1) * C['dt']) % ((360 * 24) * 3600)) == 0):
                print('Year ', (((nn + 1) * C['dt']) / ((360 * 24) * 3600)), sst[(1, int((C['nx'] / 4)))], sst[(1, int(((3 * C['nx']) / 4)))])
        tas[0, :] = tas[1, :].copy()
        sst[0, :] = sst[1, :].copy()
        hice[0, :] = hice[1, :].copy()
        sst_boundary[(n - 1)] = (sst_boundary[(n - 1)] + sst_boundary_tendency)
    hice_out[np.where((hice_out == 0))] = np.nan
    if return_coupled_fluxes:
        return (tas_out, sst_out, hice_out, sflx_out, sflx_f_out, nt1, nt, atm_DA_tendencies, ocn_DA_tendencies)
    else:
        return (tas_out, sst_out, hice_out, sflx_out, sflx_f_out, nt1, nt)