Writes npa to a HDF5 file
subroutine write_npa
!+ Writes [[libfida:npa]] to a HDF5 file
integer(HID_T) :: fid, gid
integer(HSIZE_T), dimension(3) :: dim3
integer(HSIZE_T), dimension(2) :: dim2
integer(HSIZE_T), dimension(1) :: d
integer :: error
integer, dimension(:), allocatable :: dcount
real(Float64), dimension(:,:), allocatable :: ri, rf
integer :: i, n
character(charlim) :: filename = ''
allocate(dcount(npa_chords%nchan))
do i=1,npa_chords%nchan
dcount(i) = count(npa%part%detector.eq.i)
enddo
filename=trim(adjustl(inputs%result_dir))//"/"//trim(adjustl(inputs%runid))//"_npa.h5"
!Open HDF5 interface
call h5open_f(error)
!Create file overwriting any existing file
call h5fcreate_f(filename, H5F_ACC_TRUNC_F, fid, error)
!Write Flux
d(1) = 1
dim2 = [npa%nenergy, npa%nchan]
dim3 = [npa%nenergy, npa%nchan, particles%nclass]
if(particles%nclass.gt.1) then
call h5ltmake_dataset_int_f(fid,"/nclass", 0, d, [particles%nclass], error)
call h5ltmake_compressed_dataset_double_f(fid,"/flux",3,dim3,npa%flux, error)
call h5ltset_attribute_string_f(fid,"/flux", "description", &
"Neutral flux: flux(energy,chan,class)", error)
else
call h5ltmake_compressed_dataset_double_f(fid,"/flux",2,dim3(1:2),npa%flux(:,:,1), error)
call h5ltset_attribute_string_f(fid,"/flux", "description", &
"Neutral flux: flux(energy,chan)", error)
endif
call h5ltset_attribute_string_f(fid,"/flux", "units","neutrals/(s*dE)", error)
call h5ltmake_dataset_int_f(fid,"/nenergy", 0, d, [npa%nenergy], error)
call h5ltmake_dataset_int_f(fid,"/nchan", 0, d, [npa%nchan], error)
call h5ltmake_compressed_dataset_double_f(fid,"/energy",1,dim2(1:1),&
npa%energy, error)
call h5ltmake_compressed_dataset_double_f(fid,"/radius",1,dim2(2:2),&
npa_chords%radius, error)
call h5ltmake_compressed_dataset_int_f(fid,"/count",1,dim2(2:2), dcount, error)
!Add attributes
call h5ltset_attribute_string_f(fid, "/", "version", version, error)
call h5ltset_attribute_string_f(fid,"/","description", &
"NPA flux calculated by FIDASIM",error)
call h5ltset_attribute_string_f(fid,"/nenergy","description",&
"Number of energy values",error)
call h5ltset_attribute_string_f(fid,"/nchan","description",&
"Number of channels",error)
call h5ltset_attribute_string_f(fid,"/energy","description", &
"Energy array", error)
call h5ltset_attribute_string_f(fid,"/energy","units","keV", error)
call h5ltset_attribute_string_f(fid,"/radius","description", &
"Detector line of sight radius at midplane or tangency point", error)
call h5ltset_attribute_string_f(fid,"/radius","units","cm",error)
call h5ltset_attribute_string_f(fid,"/count","description", &
"Number of particles that hit the detector: count(chan)", error)
deallocate(dcount)
if((npa%npart.ne.0).and.(inputs%calc_npa.ge.2)) then
n = npa%npart
allocate(ri(3,n),rf(3,n))
ri(1,:) = npa%part(1:n)%xi
ri(2,:) = npa%part(1:n)%yi
ri(3,:) = npa%part(1:n)%zi
rf(1,:) = npa%part(1:n)%xf
rf(2,:) = npa%part(1:n)%yf
rf(3,:) = npa%part(1:n)%zf
!Create Group
call h5gcreate_f(fid,"/particles",gid, error)
call h5ltmake_dataset_int_f(gid, "nparticle", 0, d, [npa%npart], error)
d(1) = npa%npart
dim2 = [3, n]
call h5ltmake_compressed_dataset_double_f(gid,"ri",2,dim2, ri, error)
call h5ltmake_compressed_dataset_double_f(gid,"rf",2,dim2, rf, error)
call h5ltmake_compressed_dataset_double_f(gid,"pitch",1,d, &
npa%part(1:n)%pitch, error)
call h5ltmake_compressed_dataset_double_f(gid,"energy",1,d,&
npa%part(1:n)%energy, error)
call h5ltmake_compressed_dataset_double_f(gid,"weight",1,d,&
npa%part(1:n)%weight, error)
call h5ltmake_compressed_dataset_int_f(gid,"detector",1,d,&
npa%part(1:n)%detector, error)
!Add attributes
call h5ltset_attribute_string_f(gid,"nparticle","description", &
"Number of particles that hit a detector", error)
call h5ltset_attribute_string_f(gid,"ri","description", &
"Neutral particle's birth position in machine coordinates: ri([x,y,z],particle)", error)
call h5ltset_attribute_string_f(gid,"ri","units", "cm", error)
call h5ltset_attribute_string_f(gid,"rf","description", &
"Neutral particle's hit position in machine coordinates: rf([x,y,z],particle)", error)
call h5ltset_attribute_string_f(gid,"rf","units", "cm", error)
call h5ltset_attribute_string_f(gid,"pitch","description", &
"Pitch value of the neutral particle: p = v_parallel/v w.r.t. the magnetic field", error)
call h5ltset_attribute_string_f(gid,"energy","description", &
"Energy value of the neutral particle", error)
call h5ltset_attribute_string_f(gid,"energy","units","keV",error)
call h5ltset_attribute_string_f(gid,"weight","description", &
"Neutral particle's contribution to the flux", error)
call h5ltset_attribute_string_f(gid,"weight","units","neutrals/s",error)
call h5ltset_attribute_string_f(gid,"detector","description", &
"Detector that the neutral particle hit", error)
call h5ltset_attribute_string_f(fid,"/particles","coordinate_system", &
"Right-handed cartesian",error)
call h5ltset_attribute_string_f(fid,"/particles","description", &
"Monte Carlo particles",error)
!Close group
call h5gclose_f(gid, error)
deallocate(ri,rf)
endif
!Close file
call h5fclose_f(fid, error)
!Close HDF5 interface
call h5close_f(error)
if(inputs%verbose.ge.1) then
write(*,'(T4,a,a)') 'NPA data written to: ',trim(filename)
endif
end subroutine write_npa