def SPEEDY_band_fraction(T):
    '''Python / numpy implementation of the formula used by SPEEDY and MITgcm
    to partition longwave emissions into 4 spectral bands.

    Input: temperature in Kelvin

    returns: a four-element array of band fraction

    Reproducing here the FORTRAN code from MITgcm/pkg/aim_v23/phy_radiat.F

    .. code-block:: fortran

    	      EPS3=0.95 _d 0

    	      DO JTEMP=200,320
    	        FBAND(JTEMP,0)= EPSLW
    	        FBAND(JTEMP,2)= 0.148 _d 0 - 3.0 _d -6 *(JTEMP-247)**2
    	        FBAND(JTEMP,3)=(0.375 _d 0 - 5.5 _d -6 *(JTEMP-282)**2)*EPS3
    	        FBAND(JTEMP,4)= 0.314 _d 0 + 1.0 _d -5 *(JTEMP-315)**2
    	        FBAND(JTEMP,1)= 1. _d 0 -(FBAND(JTEMP,0)+FBAND(JTEMP,2)
    	     &                           +FBAND(JTEMP,3)+FBAND(JTEMP,4))
    	      ENDDO

    	      DO JB=0,NBAND
    	        DO JTEMP=lwTemp1,199
    	          FBAND(JTEMP,JB)=FBAND(200,JB)
    	        ENDDO
    	        DO JTEMP=321,lwTemp2
    	          FBAND(JTEMP,JB)=FBAND(320,JB)
    	        ENDDO
    	      ENDDO

    '''
    # EPSLW is the fraction of longwave emission that goes directly to space
    #  It is set to zero by default in MITgcm code. We won't use it here.
    Tarray = np.array(T)
    Tarray = np.minimum(Tarray, 230.)
    Tarray = np.maximum(Tarray, 200.)
    num_band = 4
    dims = [num_band]
    dims.extend(Tarray.shape)
    FBAND = np.zeros(dims)

    EPS2=0.95
    FBAND[1,:] = 0.148 - 3.0E-6 *(T-247.)**2
    FBAND[2,:] = (0.375 - 5.5E-6 *(T-282.)**2)*EPS2
    FBAND[3,:] = 0.314 + 1.0E-5 *(T-315.)**2
    FBAND[0,:] = 1. - np.sum(FBAND, axis=0)
    return FBAND