subroutine param_fbm
  ! calculate bdens2, eba2pl, eba2pp, vbtor2 at gc
  use nbi_com
  use physconst_mod, only : evptocm_sec
  implicit none
  real(fp) :: zbmnrm,zbmn,zdelebd,zvi2,zvplov,zvplov2,zvppov2,zdnw,z_v2toev
  integer   ine,ina,inz,isb
  
  bdens2 = 0.0_fp
  eba2pl = 0.0_fp
  eba2pp = 0.0_fp
  vbtor2 = 0.0_fp
  z_v2toev = evptocm_sec*evptocm_sec
  zbmnrm=real(nznbma)/(delt)
  do inz = 1, nfbznsi
     zbmn=zbmnrm/bmvol(inz)
     do isb = 1, nsbeam
        do  ine = 1, nznbmea(isb)
           zdelebd=(efbmb(ine+1,isb)-efbmb(ine,isb))
           zvi2=(efbmb(ine+1,isb)+efbmb(ine,isb))/2.
           do ina = 1, nznbma
              zdnw=fbm(ine,ina,inz,isb)
              if(zdnw.gt.0.0_fp) then
                 zvplov=-1.0_fp+ 2.0_fp*(ina-0.5_fp)/nznbma
                 zvplov2=zvplov*zvplov
                 zvppov2=1.0_fp-zvplov2
                 bdens2(inz,isb)=bdens2(inz,isb)+zdnw*zdelebd/zbmn
                 eba2pl(inz,isb)=eba2pl(inz,isb)+zdnw*zvplov2*zvi2*zdelebd/zbmn
                 eba2pp(inz,isb)=eba2pp(inz,isb)+zdnw*zvppov2*zvi2*zdelebd/zbmn
                 vbtor2(inz,isb)=vbtor2(inz,isb)+zdnw*zvplov*&
                      &sqrt(zvi2*z_v2toev/abeams(isb))*zdelebd/zbmn
                 !mg                      &sqrt(zvi2/v2toev)*zdelebd/zbmn
                 !mg, i change v2toev=abeam/(1.384e6_fp)**2 -- conversion factor
                 ! abeam is not known here
                 
              endif
           enddo
        enddo
     enddo
  enddo
  return
end subroutine  param_fbm
