subroutine nbi_aep_get(mnkaep,nsbeam,nonlin,ierr)

  !  if there is a aep (mario podesta model) data, this routine get data from xplasma black box
  !  and put it to nbi_com
  use iso_c_binding, only: fp => c_double
  use nbi_aep_com
  use xplasma_obj_instance
  implicit none
  !--------------------------------
  !  passed:

  character*32 zobj_name  ! obj name in black box
  character*32 zauth
  character*60 zlabel
  integer,intent(in) :: mnkaep,nonlin
  integer,intent(out) :: ierr
  integer :: ido,iertmp
  integer :: ntot, npack
  !----------------------------------------
  !  pointers provided by xplasma:

  integer, dimension(:), pointer :: idata
  real(fp), dimension(:), pointer :: r8data
  integer :: idp,ie,ip,imu,ik,il,ide
  integer :: nsbeam !# of fast ions beam+fusion
  integer,dimension(:), allocatable ::  i_aepmap   !(npack) indecies mapping 
                             ! for packed probability matrix
                             ! mhd transport probability profile data
  real(fp),dimension(:), allocatable ::  aep_d   !(npack) packed probability matrix
                             ! mhd transport probability profile data
    !--------------------------------

  !
  !  start of executable code...
  !

  ierr = 0


  zauth = 'aep_transp'
  zobj_name = 'aep'
  zlabel = 'data for aep operator'

  call xplasma_blkbxid(s,zobj_name,ido)
  
  if(ido.eq.0) then
     ierr = 1
     write(nonlin,*) ' ?nbi_aep_get: "',trim(zobj_name), &
          '" object has been requested but not found.'
     return
  endif

  !-------------------
  !  data exists
  
  call xplasma_blackbox_retrieve(s,ido,ierr, ia_ptr=idata, r8a_ptr=r8data)
  if(ierr.ne.0) then
     write(nonlin,*) ' ?nbi_aep_get: unexpected error, blackbox "', &
          trim(zobj_name),'":'
     call xplasma_error(s,ierr,nonlin)
     return
  endif
 
  nkaep=mnkaep !set nbi_aep_com

  if(idata(2).ne.nkaep) then !nkaep = 0
     write(nonlin,*) &
          ' ?nbi_aep_get: data does not match, blkbx nkaep=',idata(2)
     write(nonlin,*) '            nkaep=',nkaep
     ierr = 2
     return
  endif
  ntaep=idata(3)
  ndeaep=idata(4)
  ndpaep=idata(5)
  neaep=idata(6)
  npaep=idata(7)
  nmuaep=idata(8)
  npack=idata(9)
  write(nonlin,*) 'MARINA npack=',npack
  ntot=ndeaep*ndpaep*neaep*npaep*nmuaep*nkaep
  if(allocated(aepsp)) deallocate(aepsp)
  allocate(aepsp(nkaep))
  if(allocated(aept)) deallocate(aept)
  allocate(aept(ntaep))
  if(allocated(aepa)) deallocate(aepa)
  allocate(aepa(ntaep,nkaep))
  if(allocated(aeptstep)) deallocate(aeptstep)
  allocate(aeptstep(nkaep))
  if(allocated(aepde)) deallocate(aepde)
  allocate(aepde(nkaep,3))
  if(allocated(aepdp)) deallocate(aepdp)
  allocate(aepdp(nkaep,3))
  if(allocated(aepe)) deallocate(aepe)
  allocate(aepe(nkaep,3))
  if(allocated(aepp)) deallocate(aepp)
  allocate(aepp(nkaep,3))
  if(allocated(aepmu)) deallocate(aepmu)
  allocate(aepmu(nkaep,3))
  if(npack.gt.0) then
     if(allocated(i_aepmap)) deallocate(i_aepmap)
     allocate(i_aepmap(npack))
     if(allocated(aep_d)) deallocate(aep_d)
     allocate(aep_d(ntot))
  endif
  il=1
  aepsp(1:nkaep)=r8data(il:il+nkaep-1) 
  il=il+nkaep
  aept(1:ntaep)=r8data(il:il+ntaep-1)
  il=il+ntaep
  do ik=1,nkaep
     aepa(1:ntaep,ik)=r8data(il:il+ntaep-1)
     il=il+ntaep
  enddo
  aeptstep(1:nkaep)=r8data(il:il+nkaep-1)
  il=il+nkaep
  do ik=1,3
     aepde(1:nkaep,ik)=r8data(il:il+nkaep-1)
     il=il+nkaep
  enddo
  do ik=1,3
     aepdp(1:nkaep,ik)=r8data(il:il+nkaep-1)
     il=il+nkaep
  enddo
  do ik=1,3
     aepe(1:nkaep,ik)=r8data(il:il+nkaep-1)
     il=il+nkaep
  enddo
  do ik=1,3
     aepp(1:nkaep,ik)=r8data(il:il+nkaep-1)
     il=il+nkaep
  enddo
  do ik=1,3
     aepmu(1:nkaep,ik)=r8data(il:il+nkaep-1)
     il=il+nkaep
  enddo
  if(npack.gt.0) then
     i_aepmap(1:npack)=int(r8data(il:il+npack-1))
     il=il+npack
     aep_d=0.0_fp
     aep_d(1:npack)=r8data(il:il+npack-1)
     il=il+npack
     do ik=npack,1,-1
        if(ik.ne.i_aepmap(ik)) then
           aep_d(i_aepmap(ik))=aep_d(ik) 
           aep_d(ik)=0.0_fp
        endif
     enddo
     deallocate(i_aepmap)
     if(allocated(aepd)) deallocate(aepd)
     allocate(aepd(ndeaep,ndpaep,neaep,npaep,nmuaep,nkaep))
     aepd=0.0_fp
     aepd=reshape((aep_d),(/ndeaep,ndpaep,neaep,npaep,nmuaep,nkaep/))
     if(allocated(aep_d)) deallocate(aep_d)
  else
! no packing has been used in old version of nubeam
     if(allocated(aepd)) deallocate(aepd)
     allocate(aepd(ndeaep,ndpaep,neaep,npaep,nmuaep,nkaep))
     aepd=0.0_fp
     do ik=1,nkaep
        do imu=1,nmuaep
           do ip=1,npaep
              do ie=1,neaep
                 do idp=1,ndpaep
                    aepd(1:ndeaep,idp,ie,ip,imu,ik)=r8data(il:il+ndeaep-1)
                    il=il+ndeaep
                 enddo
              enddo
           enddo
        enddo
     enddo
  endif
!  all done.

  nullify(idata)
  nullify(r8data)
 
end subroutine nbi_aep_get
