subroutine writefs(ef0)!Fermi surface mode (eigenvalues in full BZ). No output variables.
  use m_lmfinit, only: nsp,nspc
  use m_lattic,only: qlat=>lat_qlat,plat=>lat_plat
  use m_mkqp,only: bz_nabc
  use m_suham,only: ndhamx=>ham_ndhamx,nspx=>ham_nspx
  use m_qplist,only:nkp,qplist
  use m_bandcal,only: evlall
  use m_shortn3_qlat,only: shortn3_qlat,nout,nlatout
  implicit none
  logical:: cmdopt0,allband
  real(8):: ppin(3),ef0
  integer:: ip,i,isp,ififm,nbxx,iq,ib,nkk1,nkk2,nkk3,ifi
  real(8):: rlatp(3,3),xmx2(3),vadd,qshort(3)
  integer:: iout 
  nkk1=bz_nabc(1)
  nkk2=bz_nabc(2)
  nkk3=bz_nabc(3)
  allband  = cmdopt0('--allband')
  i=nsp
  if(nspc==2) i=1
  do isp=1,i
     if(isp==1) open(newunit=ifi, file='fermiup.bxsf')
     if(isp==2) open(newunit=ifi, file='fermidn.bxsf')
     if(isp==1) open(newunit=ififm, file='fermiup.data')
     if(isp==2) open(newunit=ififm, file='fermidn.data')
     write(ififm,*) " # Generated by bndfp.F with --fermisurface"
     write(ififm,*) " # q-point is reduced in 1st BZ"
     write(ifi,*) "BEGIN_INFO"
     write(ifi,*) " # usage xcrysden --bxsf fermiup.bxsf"
     write(ifi,*) " # http://www.xcrysden.org/doc/XSF.html#2l.16"
     write(ifi,'(a,f9.5)') "  Fermi Energy:",ef0
     write(ififm,'(a,f9.5)') "  # Fermi Energy [eV]:",ef0
     write(ififm,*) " # qshort(3) band_energy[eV]"
     write(ifi,*) "END_INFO"
     write(ifi,*)"BEGIN_BLOCK_BANDGRID_3D"
     write(ifi,*)"  this_is_for_xcrysden"
     write(ifi,*)"  BEGIN_BANDGRID_3D_simple_test"
     nbxx=0
     do ib=1,ndhamx
        if(allband .OR. (minval(evlall(ib,isp,:))<ef0+0.5 .AND. &
             maxval(evlall(ib,isp,:))>ef0-0.5)) then
           nbxx=nbxx+1
        endif
     enddo
     write(ifi,"(4x,i8)") nbxx
     write(ifi,"(4x,3i8)") nkk1,nkk2,nkk3
     write(ifi,"(4x,3(f9.5,1x))") 0d0,0d0,0d0
     write(ifi,"(4x,3(f9.5,1x))") qlat(:,1)
     write(ifi,"(4x,3(f9.5,1x))") qlat(:,2)
     write(ifi,"(4x,3(f9.5,1x))") qlat(:,3)
     !call shortn3_initialize(qlat)
     do ib=1,ndhamx
        if(allband .OR. (minval(evlall(ib,isp,:))<ef0+0.5 .AND. &
             maxval(evlall(ib,isp,:))>ef0-0.5)) then
           write(ifi,"(a,i8)")"  BAND: ",ib
           write(ififm,"(a,i8)")"  # BAND: ",ib
           write(ifi,"(3x,10(x,f9.5))") &
                (evlall(ib,isp,iq), iq=1,nkk1*nkk2*nkk3)
           ! to reduce q-point to qshort
           do iq=1,nkk1*nkk2*nkk3
              ppin=matmul(transpose(plat),qplist(:,iq))
              call shortn3_qlat(ppin) !bug (qlist(:,iq) ==>ppin 2022-6-8 tkotani)
              iout=1
              qshort(1:3)=  matmul(qlat(:,:), ppin+nlatout(:,iout))
              write(ififm,"(4f13.5)") qshort(:),evlall(ib,isp,iq) !q-point reduced to 1stBZ
           enddo
        endif
     enddo
     write(ifi,*)
     write(ifi,*)"  END_BANDGRID_3D"
     write(ifi,*)"END_BLOCK_BANDGRID_3D"
  enddo
  close(ifi)
  close(ififm)
end subroutine writefs
