docs: 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