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
real(Float64), dimension(:), allocatable :: weight, energy, pitch
integer, dimension(:), allocatable :: det, orbit_type
integer :: i, npart, c, start_index, end_index
character(charlim) :: filename = ''
logical :: do_write = .True.
#ifdef _MPI
integer :: rank
integer, dimension(:), allocatable :: npart_image
rank = my_rank()
allocate(npart_image(0:num_ranks()-1))
#endif
#ifdef _MPI
if(my_rank().ne.0) do_write = .False.
#endif
filename=trim(adjustl(inputs%result_dir))//"/"//trim(adjustl(inputs%runid))//"_npa.h5"
if(do_write) then
!Open HDF5 interface
call h5open_f(error)
!Create file overwriting any existing file
call h5fcreate_f(filename, H5F_ACC_TRUNC_F, fid, error)
endif
!! Active
allocate(dcount(npa_chords%nchan))
npart = npa%npart
if(npart.gt.0) then
do i=1,npa_chords%nchan
dcount(i) = count(npa%part%detector.eq.i)
enddo
else
dcount = 0
endif
#ifdef _MPI
call parallel_sum(dcount)
call parallel_sum(npart)
#endif
c = 0
#ifdef _MPI
npart_image(:) = 0
npart_image(rank) = npa%npart
call parallel_sum(npart_image)
do i=0,rank-1
c = c + npart_image(i)
enddo
#endif
start_index = 1 + c
end_index = npa%npart + c
if(npart.gt.0) then
allocate(ri(3,npart),rf(3,npart))
allocate(weight(npart),energy(npart),pitch(npart))
allocate(det(npart),orbit_type(npart))
ri = 0.d0 ; rf = 0.d0
weight = 0.d0 ; energy = 0.d0
pitch = 0.d0 ; det = 0
orbit_type = 0
c = 1
do i=start_index, end_index
ri(1,i) = npa%part(c)%xi
ri(2,i) = npa%part(c)%yi
ri(3,i) = npa%part(c)%zi
rf(1,i) = npa%part(c)%xf
rf(2,i) = npa%part(c)%yf
rf(3,i) = npa%part(c)%zf
weight(i) = npa%part(c)%weight
energy(i) = npa%part(c)%energy
pitch(i) = npa%part(c)%pitch
det(i) = npa%part(c)%detector
orbit_type(i) = npa%part(c)%class
c = c + 1
enddo
#ifdef _MPI
call parallel_sum(ri)
call parallel_sum(rf)
call parallel_sum(weight)
call parallel_sum(energy)
call parallel_sum(pitch)
call parallel_sum(det)
call parallel_sum(orbit_type)
#endif
endif
if(do_write.and.(inputs%calc_npa.ge.1)) then
!Write Active 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", &
"Active 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", &
"Active 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)
if((npart.gt.0).and.(inputs%calc_npa.ge.2)) then
!Create Group
call h5gcreate_f(fid,"/particles",gid, error)
call h5ltmake_dataset_int_f(gid, "nparticle", 0, d, [npart], error)
d(1) = npart
dim2 = [3, npart]
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, pitch, error)
call h5ltmake_compressed_dataset_double_f(gid,"energy",1,d, energy, error)
call h5ltmake_compressed_dataset_double_f(gid,"weight",1,d, weight, error)
call h5ltmake_compressed_dataset_int_f(gid,"detector",1,d, det, error)
call h5ltmake_compressed_dataset_int_f(gid,"class",1,d, orbit_type, 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(gid,"class","description", &
"Class of the neutral particle", error)
call h5ltset_attribute_string_f(fid,"/particles","coordinate_system", &
"Right-handed cartesian",error)
call h5ltset_attribute_string_f(fid,"/particles","description", &
"Active NPA Monte Carlo particles",error)
!Close group
call h5gclose_f(gid, error)
endif
endif
deallocate(dcount)
if(npart.gt.0) then
deallocate(ri,rf)
deallocate(energy,pitch,weight,det,orbit_type)
endif
!! Passive
allocate(dcount(npa_chords%nchan))
npart = pnpa%npart
if(npart.gt.0) then
do i=1,npa_chords%nchan
dcount(i) = count(pnpa%part%detector.eq.i)
enddo
else
dcount = 0
endif
#ifdef _MPI
call parallel_sum(dcount)
call parallel_sum(npart)
#endif
c = 0
#ifdef _MPI
npart_image(:) = 0
npart_image(rank) = pnpa%npart
call parallel_sum(npart_image)
do i=0,rank-1
c = c + npart_image(i)
enddo
#endif
start_index = 1 + c
end_index = pnpa%npart + c
if(npart.gt.0) then
allocate(ri(3,npart),rf(3,npart))
allocate(weight(npart),energy(npart),pitch(npart))
allocate(det(npart),orbit_type(npart))
ri = 0.d0 ; rf = 0.d0
weight = 0.d0 ; energy = 0.d0
pitch = 0.d0 ; det = 0
orbit_type = 0
c = 1
do i=start_index, end_index
ri(1,i) = pnpa%part(c)%xi
ri(2,i) = pnpa%part(c)%yi
ri(3,i) = pnpa%part(c)%zi
rf(1,i) = pnpa%part(c)%xf
rf(2,i) = pnpa%part(c)%yf
rf(3,i) = pnpa%part(c)%zf
weight(i) = pnpa%part(c)%weight
energy(i) = pnpa%part(c)%energy
pitch(i) = pnpa%part(c)%pitch
det(i) = pnpa%part(c)%detector
orbit_type(i) = pnpa%part(c)%class
c = c + 1
enddo
#ifdef _MPI
call parallel_sum(ri)
call parallel_sum(rf)
call parallel_sum(weight)
call parallel_sum(energy)
call parallel_sum(pitch)
call parallel_sum(det)
call parallel_sum(orbit_type)
#endif
endif
if(do_write.and.(inputs%calc_pnpa.ge.1)) then
!Write Passive Flux
d(1) = 1
dim2 = [pnpa%nenergy, pnpa%nchan]
dim3 = [pnpa%nenergy, pnpa%nchan, particles%nclass]
if(particles%nclass.gt.1) then
call h5ltmake_compressed_dataset_double_f(fid,"/pflux",3,dim3,pnpa%flux, error)
call h5ltset_attribute_string_f(fid,"/pflux", "description", &
"Passive Neutral flux: pflux(energy,chan,class)", error)
else
call h5ltmake_compressed_dataset_double_f(fid,"/pflux",2,dim3(1:2),pnpa%flux(:,:,1), error)
call h5ltset_attribute_string_f(fid,"/pflux", "description", &
"Passive Neutral flux: pflux(energy,chan)", error)
endif
call h5ltset_attribute_string_f(fid,"/pflux", "units","neutrals/(s*dE)", error)
call h5ltmake_compressed_dataset_int_f(fid,"/pcount",1,dim2(2:2), dcount, error)
call h5ltset_attribute_string_f(fid,"/pcount","description", &
"Number of passive particles that hit the detector: pcount(chan)", error)
if(inputs%calc_npa.le.0) then
call h5ltmake_dataset_int_f(fid,"/nenergy", 0, d, [pnpa%nenergy], error)
call h5ltmake_dataset_int_f(fid,"/nchan", 0, d, [pnpa%nchan], error)
call h5ltmake_compressed_dataset_double_f(fid,"/energy",1,dim2(1:1),&
pnpa%energy, error)
call h5ltmake_compressed_dataset_double_f(fid,"/radius",1,dim2(2:2),&
npa_chords%radius, 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)
endif
if((npart.gt.0).and.(inputs%calc_pnpa.ge.2)) then
!Create Group
call h5gcreate_f(fid,"/passive_particles",gid, error)
call h5ltmake_dataset_int_f(gid, "nparticle", 0, d, [npart], error)
d(1) = npart
dim2 = [3, npart]
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, pitch, error)
call h5ltmake_compressed_dataset_double_f(gid,"energy",1,d, energy, error)
call h5ltmake_compressed_dataset_double_f(gid,"weight",1,d, weight, error)
call h5ltmake_compressed_dataset_int_f(gid,"detector",1,d, det, error)
call h5ltmake_compressed_dataset_int_f(gid,"class",1,d, orbit_type, 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(gid,"class","description", &
"Class of the neutral particle", error)
call h5ltset_attribute_string_f(fid,"/passive_particles","coordinate_system", &
"Right-handed cartesian",error)
call h5ltset_attribute_string_f(fid,"/passive_particles","description", &
"Passive NPA Monte Carlo particles",error)
!Close group
call h5gclose_f(gid, error)
endif
endif
if(do_write) then
!Close file
call h5fclose_f(fid, error)
!Close HDF5 interface
call h5close_f(error)
endif
if(inputs%verbose.ge.1) then
write(*,'(T4,a,a)') 'NPA data written to: ',trim(filename)
endif
#ifdef _MPI
deallocate(npart_image)
#endif
end subroutine write_npa