def net_radiation(K, cloudcover, T, albedo=0.23):
    """Computes the net radiation based on incoming solar radiation, cloud cover
    and the temperature.

    Parameters
    ----------
    K : float
        Incoming solar radiation in W/m^2
    cloudcover : float
        In percent (0-1)
    T : float
        Temperature in Kelvin
    albedo : float (optional)
        Albedo of the surface

    Other Parameters
    ----------------
    L_plus : float (internal)
        Incoming longwave
    L_minus : float (internal)
        Outgoing longwave

    Returns
    -------
    Q : float
        Net radiation in W/m^2
    """

    if cloudcover > 1:
        cloudcover = 1
    if cloudcover < 0:
        cloudcover = 0

    c1 = 5.31 * 10**(-13)  # W m**-2 K**-6
    c2 = 60  # W m**-2
    c3 = 0.12
    sigma = 5.67 * 10**(-8)  # W m**-2 K**-4 Stefan Boltzmann constant

    L_plus = c1 * T**6 + c2 * cloudcover
    #L_minus = sigma * T**4 + c3*Q

    Q = ((1 - albedo) * K + L_plus - sigma * T**4) / (1 + c3)

    return Q