subroutine nubeam2toric_wr(fname)
  use nbi_com
  use xplasma_obj_instance
  use ezspline_obj, only : ezspline1
  use ezspline, only : ezspline_free,ezspline_init,ezspline_error,&
       ezspline_setup,ezspline_interp
  use netcdf
  implicit none
  integer :: izone,izone_toric,i,ierr,icdf,isb,idims(1)
  integer :: id_nzone_toric,id_sqrt_nrm_plf,id_bdenss_sqrnrmplf,id_ebavps_sqrnrmplf
  integer :: id_prf_sqrnrmplf,id_toric_prf_sqrnrmplf,id_ebapp2s_sqrnrmplf,id_ebapl2s_sqrnrmplf
  real(fp),allocatable,dimension(:) :: rho_wrk !sqrt(norm tor flux)
  real(fp),allocatable,dimension(:) :: plf_wrk !poloidal flux
  real(fp),allocatable,dimension(:) :: sqr_nrmplf !sqrt(norm poloidal flux)
  real(fp),allocatable,dimension(:) :: bdenss_sqrnrmplf !bdenss on sqrt(norm poloidal flux)
  real(fp),allocatable,dimension(:) :: prf_sqrnrmplf !absoption power from nubeam on sqrt(norm poloidal flux)
  real(fp),allocatable,dimension(:) :: toric_prf_sqrnrmplf !absoption power from toric on sqrt(norm poloidal flux)
  real(fp),allocatable,dimension(:) :: ebapp2s_sqrnrmplf !ebapp2s on sqrt(norm poloidal flux)
  real(fp),allocatable,dimension(:) :: ebapl2s_sqrnrmplf !ebapl2s on sqrt(norm poloidal flux)
  real(fp),allocatable,dimension(:) :: ebavps_sqrnrmplf !ebavps on sqrt(norm poloidal flux)
  real(fp),allocatable,dimension(:,:) :: wrk
  real(fp) :: psi_lcfs
  real(fp), parameter :: tol_max=1.0e-8 ! [nominal: 1d-8]

  character*(*), intent(in) ::fname !netcdf file name for output data
  type(ezspline1) :: spl_bdenss ! 1d interpolation
  type(ezspline1) :: spl_ebapp2s ! 1d interpolation
  type(ezspline1) :: spl_ebapl2s ! 1d interpolation
  type(ezspline1) :: spl_ebavps ! 1d interpolation
  type(ezspline1) :: spl_prf ! 1d interpolation
  type(ezspline1) :: spl_toric_prf ! 1d interpolation
  integer :: bcs_sig(2)
  logical :: il_prf

  isb=1 !species #, only first one to be printed out
  izone=ledge-lcentr+1
  izone_toric=2*izone+1
  allocate(rho_wrk(izone_toric))
  allocate(plf_wrk(izone_toric))
  allocate(sqr_nrmplf(izone_toric))
  allocate(bdenss_sqrnrmplf(izone_toric))
  allocate(ebavps_sqrnrmplf(izone_toric))
  allocate(ebapp2s_sqrnrmplf(izone_toric))
  allocate(ebapl2s_sqrnrmplf(izone_toric))
  il_prf=.false.
  if(maxval(abs(orbrzv_prf(:,isb))).gt.0.0_fp)then
     il_prf=.true.
  endif
  if(il_prf) then
     allocate(prf_sqrnrmplf(izone_toric))
     allocate(toric_prf_sqrnrmplf(izone_toric))
     prf_sqrnrmplf=0.0_fp
     toric_prf_sqrnrmplf=0.0_fp
  endif
  allocate(wrk(lep1,2))
  
  rho_wrk=0._fp
  plf_wrk=0.0_fp
  sqr_nrmplf=0.0_fp
  bdenss_sqrnrmplf=0.0_fp
  ebavps_sqrnrmplf=0.0_fp
  ebapp2s_sqrnrmplf=0.0_fp
  ebapl2s_sqrnrmplf=0.0_fp
  do i=2,izone_toric
     rho_wrk(i)= rho_wrk(i-1)+1._fp/(izone_toric-1.0_fp)!at zone boundary
  enddo
  sqr_nrmplf=rho_wrk
  call eq_gfnum('PSI',id_psi)
  call eq_rgetf(izone_toric, rho_wrk, id_psi, 0, plf_wrk, ierr)
  
  psi_lcfs=maxval(abs(plf_wrk))
  wrk(:,:)=0.0_fp ! build profile at zone boundary 
  wrk(lcentr:ledge,2)=ebavps(lcentr:ledge,isb)
  wrk(lcentr,1)=wrk(lcentr,2)
  wrk(lcp1:ledge,1)=0.5_fp*(wrk(lcentr:ledge-1,2)+ &
       wrk(lcp1:ledge,2))
  wrk(lep1,1)=max(0.5_fp*wrk(ledge,2), &
       2.0_fp*wrk(ledge,2)-wrk(ledge,1))
  bcs_sig = (/0, 0/)     ! not-a-knot boundary conditions
  call ezspline_free(spl_ebavps,ierr)
  call ezspline_init(spl_ebavps,izone+1,bcs_sig,ierr)      ! initiali
  call ezspline_error(ierr)
  spl_ebavps%x1(1:izone+1)=xi(lcentr:lep1,1)
  call ezspline_setup(spl_ebavps,wrk(lcentr:lep1,1),ierr) ! initialize splin
  call ezspline_error(ierr)

  wrk(:,:)=0.0_fp ! build profile at zone boundary 
  wrk(lcentr:ledge,2)=ebapp2s(lcentr:ledge,isb)
  wrk(lcentr,1)=wrk(lcentr,2)
  wrk(lcp1:ledge,1)=0.5_fp*(wrk(lcentr:ledge-1,2)+ &
       wrk(lcp1:ledge,2))
  wrk(lep1,1)=max(0.5_fp*wrk(ledge,2), &
       2.0_fp*wrk(ledge,2)-wrk(ledge,1))
  bcs_sig = (/0, 0/)     ! not-a-knot boundary conditions
  call ezspline_free(spl_ebapp2s,ierr)
  call ezspline_init(spl_ebapp2s,izone+1,bcs_sig,ierr)      ! initiali
  call ezspline_error(ierr)
  spl_ebapp2s%x1(1:izone+1)=xi(lcentr:lep1,1)
  call ezspline_setup(spl_ebapp2s,wrk(lcentr:lep1,1),ierr) ! initialize splin
  call ezspline_error(ierr)

  wrk(:,:)=0.0_fp ! build profile at zone boundary 
  wrk(lcentr:ledge,2)=ebapl2s(lcentr:ledge,isb)
  wrk(lcentr,1)=wrk(lcentr,2)
  wrk(lcp1:ledge,1)=0.5_fp*(wrk(lcentr:ledge-1,2)+ &
       wrk(lcp1:ledge,2))
  wrk(lep1,1)=max(0.5_fp*wrk(ledge,2), &
       2.0_fp*wrk(ledge,2)-wrk(ledge,1))
  bcs_sig = (/0, 0/)     ! not-a-knot boundary conditions
  call ezspline_free(spl_ebapl2s,ierr)
  call ezspline_init(spl_ebapl2s,izone+1,bcs_sig,ierr)      ! initiali
  call ezspline_error(ierr)
  spl_ebapl2s%x1(1:izone+1)=xi(lcentr:lep1,1)
  call ezspline_setup(spl_ebapl2s,wrk(lcentr:lep1,1),ierr) ! initialize splin
  call ezspline_error(ierr)


  call ezspline_free(spl_bdenss,ierr)
  call ezspline_init(spl_bdenss,izone+1,bcs_sig,ierr)      ! initiali
  call ezspline_error(ierr)
  spl_bdenss%x1(1:izone+1)=xi(lcentr:lep1,1) !at zone boundary
  wrk(:,:)=0.0_fp ! build profile at zone boundary 
  wrk(lcentr:ledge,2)=bdenss(lcentr:ledge,isb)/dvolzc(lcentr:ledge)! beam density #/cm-3
  wrk(lcentr,1)=wrk(lcentr,2)
  wrk(lcp1:ledge,1)=0.5_fp*(wrk(lcentr:ledge-1,2)+ &
       wrk(lcp1:ledge,2))
  wrk(lep1,1)=max(0.5_fp*wrk(ledge,2), &
       2.0_fp*wrk(ledge,2)-wrk(ledge,1))
  call ezspline_setup(spl_bdenss,wrk(lcentr:lep1,1),ierr) ! initialize splin
  call ezspline_error(ierr)
  call ezspline_free(spl_prf,ierr)
  call ezspline_init(spl_prf,izone+1,bcs_sig,ierr)      ! initiali
  call ezspline_error(ierr)
  spl_prf%x1(1:izone+1)=xi(lcentr:lep1,1) !at zone boundary
  wrk(:,:)=0.0_fp ! build profile at zone boundary 
  wrk(lcentr:ledge,2)=orbrzv_prf(lcentr:ledge,isb) ! power w
  wrk(lcentr,1)=wrk(lcentr,2)
  wrk(lcp1:ledge,1)=0.5_fp*(wrk(lcentr:ledge-1,2)+ &
       wrk(lcp1:ledge,2))
  wrk(lep1,1)=max(0.5_fp*wrk(ledge,2), &
       2.0_fp*wrk(ledge,2)-wrk(ledge,1))
  call ezspline_setup(spl_prf,wrk(lcentr:lep1,1),ierr) ! initialize splin
  call ezspline_error(ierr)
  call ezspline_free(spl_toric_prf,ierr)
  call ezspline_init(spl_toric_prf,izone+1,bcs_sig,ierr)      ! initiali
  call ezspline_error(ierr)
  spl_prf%x1(1:izone+1)=xi(lcentr:lep1,1) !at zone boundary
  wrk(:,:)=0.0_fp ! build profile at zone boundary 
  wrk(lcentr:ledge,2)=orbrzv_toric_prf(lcentr:ledge,isb) ! power w
  wrk(lcentr,1)=wrk(lcentr,2)
  wrk(lcp1:ledge,1)=0.5_fp*(wrk(lcentr:ledge-1,2)+ &
       wrk(lcp1:ledge,2))
  wrk(lep1,1)=max(0.5_fp*wrk(ledge,2), &
       2.0_fp*wrk(ledge,2)-wrk(ledge,1))
  call ezspline_setup(spl_toric_prf,wrk(lcentr:lep1,1),ierr) ! initialize splin
  call ezspline_error(ierr)
  rho_wrk=0.0_fp
  call eq_rhopsi(izone_toric,psi_lcfs*sqr_nrmplf(:)*sqr_nrmplf(:),rho_wrk,tol_max,ierr)
  do i=1,izone_toric
     call ezspline_interp(spl_bdenss,rho_wrk(i),bdenss_sqrnrmplf(i),ierr)
     call ezspline_interp(spl_ebavps,rho_wrk(i),ebavps_sqrnrmplf(i),ierr)
     call ezspline_interp(spl_ebapp2s,rho_wrk(i),ebapp2s_sqrnrmplf(i),ierr)
     call ezspline_interp(spl_ebapl2s,rho_wrk(i),ebapl2s_sqrnrmplf(i),ierr)
  enddo
     
  if(il_prf) then
     do i=1,izone_toric
        call ezspline_interp(spl_prf,rho_wrk(i),prf_sqrnrmplf(i),ierr)
        call ezspline_interp(spl_toric_prf,rho_wrk(i),toric_prf_sqrnrmplf(i),ierr)
     enddo
  endif

  ierr = nf90_create(path=trim(fname),cmode= or(nf90_clobber,nf90_64bit_offset),ncid=icdf)
  call check(nf90_def_dim(icdf,'nzone_toric',izone_toric,id_nzone_toric),'nzone_torici')

  idims = (/id_nzone_toric/)

  call check(nf90_def_var(icdf,'sqrt_nrm_plf',nf90_double,idims,id_sqrt_nrm_plf),'sqrt_nrm_plf')
  call check(nf90_def_var(icdf,'bdenss_sqrnrmplf',nf90_double,idims,id_bdenss_sqrnrmplf),'bdenss_sqrnrmplf')
  call check(nf90_def_var(icdf,'ebavps_sqrnrmplf',nf90_double,idims,id_ebavps_sqrnrmplf),'ebavps_sqrnrmplf')
  call check(nf90_def_var(icdf,'ebapp2s_sqrnrmplf',nf90_double,idims,id_ebapp2s_sqrnrmplf),'ebapp2s_sqrnrmplf')
  call check(nf90_def_var(icdf,'ebapl2s_sqrnrmplf',nf90_double,idims,id_ebapl2s_sqrnrmplf),'ebapl2s_sqrnrmplf')
  call check(nf90_put_att(icdf,id_sqrt_nrm_plf,'sqrt_nrm_plf',"sqrt(norm poloidal flux)"),'sqrt_nrm_plf')
  ierr=nf90_put_att(icdf,id_bdenss_sqrnrmplf,'bdenss_sqrnrmplf',"bdenss on sqrt(norm poloidal flux) at zone boundary,1/cm-3")
  ierr=nf90_put_att(icdf,id_ebavps_sqrnrmplf,'ebavps_sqrnrmplf',"ebavps on sqrt(norm poloidal flux) at zone boundary,ev")
  ierr=nf90_put_att(icdf,id_ebapp2s_sqrnrmplf,'ebapp2s_sqrnrmplf',"ebapp2s on sqrt(norm poloidal flux) at zone boundary,ev")
  ierr=nf90_put_att(icdf,id_ebapl2s_sqrnrmplf,'ebapl2s_sqrnrmplf',"ebapl2s on sqrt(norm poloidal flux) at zone boundary,ev")
  if(il_prf) then
     call check(nf90_def_var(icdf,'prf_sqrnrmplf',nf90_double,idims,id_prf_sqrnrmplf),'prf_sqrnrmplf')
     call check(nf90_def_var(icdf,'toric_prf_sqrnrmplf',nf90_double,idims,id_toric_prf_sqrnrmplf),&
          'toric_prf_sqrnrmplf')
     ierr=nf90_put_att(icdf,id_prf_sqrnrmplf,'prf_sqrnrmplf',&
          "rf power from nubeam on sqrt(norm poloidal flux) at zone boundary, w")
     ierr=nf90_put_att(icdf,id_toric_prf_sqrnrmplf,'toric_prf_sqrnrmplf',&
          "rf power from toric on sqrt(norm poloidal flux) at zone boundary, w")
  endif
  ierr = nf90_enddef(icdf) ! out of def mode and go to data mode
  call check( nf90_put_var(icdf,id_sqrt_nrm_plf,sqr_nrmplf),'sqrt_nrm_plf')
  call check( nf90_put_var(icdf,id_bdenss_sqrnrmplf,bdenss_sqrnrmplf),'bdenss_sqrnrmplf')
  call check( nf90_put_var(icdf,id_ebavps_sqrnrmplf,ebavps_sqrnrmplf),'ebavps_sqrnrmplf')
  call check( nf90_put_var(icdf,id_ebapp2s_sqrnrmplf,ebapp2s_sqrnrmplf),'ebapp2s_sqrnrmplf')
  call check( nf90_put_var(icdf,id_ebapl2s_sqrnrmplf,ebapl2s_sqrnrmplf),'ebapl2s_sqrnrmplf')
  if(il_prf) then
     call check( nf90_put_var(icdf,id_prf_sqrnrmplf,prf_sqrnrmplf),'prf_sqrnrmplf')
     call check( nf90_put_var(icdf,id_toric_prf_sqrnrmplf,toric_prf_sqrnrmplf),'toric_prf_sqrnrmplf')
  endif
  ierr = nf90_close( icdf)
  deallocate(rho_wrk,plf_wrk,bdenss_sqrnrmplf,ebavps_sqrnrmplf,sqr_nrmplf,wrk)
  deallocate(ebapp2s_sqrnrmplf,ebapl2s_sqrnrmplf)
  if(il_prf) then
     deallocate(prf_sqrnrmplf,toric_prf_sqrnrmplf)
  endif

contains 

  subroutine check(status,varnam)
    integer, intent ( in) :: status 
    character*(*),intent(in) :: varnam ! dimension name input 
 
    if(status /= nf90_noerr) then 
       write(nonlin,*) ' ?nubeam2toric_wr:'& 
      &// trim(nf90_strerror(status))//' for var.= '//trim(varnam) 
       call exit(status) 
    end if 
  end subroutine check 

end subroutine nubeam2toric_wr
