read_mc Subroutine

subroutine read_mc(fid, error)
    !+ Reads in a MC particle fast-ion distribution and puts them in [[libfida:particles]]
    integer(HID_T), intent(inout) :: fid
        !+ HDF5 file ID
    integer, intent(out)          :: error
        !+ Error code

    integer(HSIZE_T), dimension(1) :: dims
    integer(Int32) :: i,j,ii,ir,iz
    real(Float64) :: phi,phi_enter,phi_exit,delta_phi
    real(Float64), dimension(3) :: uvw,ri,vi,e1_xyz,e2_xyz,C_xyz,dum
    integer(Int32), dimension(1) :: minpos
    real(Float64), dimension(:), allocatable :: weight
    type(LocalEMFields) :: fields
    integer :: cnt,num
    logical :: inp
    character(len=32) :: dist_type_name = ''

    if(inputs%verbose.ge.1) then
        write(*,'(a)') '---- Fast-ion distribution settings ----'
    endif

    call h5ltread_dataset_int_scalar_f(fid, "/nparticle", particles%nparticle, error)
    call h5ltread_dataset_int_scalar_f(fid, "/nclass", particles%nclass, error)

    !!ALLOCATE SPACE
    allocate(particles%fast_ion(particles%nparticle))
    allocate(weight(particles%nparticle))

    dims(1) = particles%nparticle
    call h5ltread_dataset_double_f(fid, "/r", particles%fast_ion%r, dims, error)
    call h5ltread_dataset_double_f(fid, "/z", particles%fast_ion%z, dims, error)
    call h5ltread_dataset_int_f(fid, "/class", particles%fast_ion%class, dims, error)

    if(any(particles%fast_ion%class.gt.particles%nclass)) then
        if(inputs%verbose.ge.0) then
            write(*,'(a)') 'READ_MC: Orbit class ID greater then the number of classes'
        endif
        stop
    endif

    if(inputs%dist_type.eq.2) then
        dist_type_name = "Guiding Center Monte Carlo"
        call h5ltread_dataset_double_f(fid, "/energy", particles%fast_ion%energy, dims, error)
        call h5ltread_dataset_double_f(fid, "/pitch", particles%fast_ion%pitch, dims, error)
        particles%fast_ion%vabs  = sqrt(particles%fast_ion%energy/(v2_to_E_per_amu*inputs%ab))
    else
        dist_type_name = "Full Orbit Monte Carlo"
        call h5ltread_dataset_double_f(fid, "/vr", particles%fast_ion%vr, dims, error)
        call h5ltread_dataset_double_f(fid, "/vt", particles%fast_ion%vt, dims, error)
        call h5ltread_dataset_double_f(fid, "/vz", particles%fast_ion%vz, dims, error)
        particles%fast_ion%vabs = sqrt(particles%fast_ion%vr**2 + &
                                       particles%fast_ion%vt**2 + &
                                       particles%fast_ion%vz**2)
        particles%fast_ion%energy = v2_to_E_per_amu*inputs%ab*particles%fast_ion%vabs**2
    endif

    call h5ltread_dataset_double_f(fid, "/weight", weight, dims, error)

    cnt=0
    e1_xyz = matmul(beam_grid%inv_basis,[1.0,0.0,0.0])
    e2_xyz = matmul(beam_grid%inv_basis,[0.0,1.0,0.0])
    !$OMP PARALLEL DO schedule(guided) private(i,ii,j,ir,iz,minpos,fields,uvw,phi,ri,vi, &
    !$OMP& delta_phi,phi_enter,phi_exit,C_xyz)
    particle_loop: do i=1,particles%nparticle
        if(inputs%verbose.ge.2) then
            WRITE(*,'(f7.2,"% completed",a,$)') cnt/real(particles%nparticle)*100,char(13)
        endif
        uvw = [particles%fast_ion(i)%r, 0.d0, particles%fast_ion(i)%z]
        call in_plasma(uvw,inp,machine_coords=.True.)
        if(.not.inp) cycle particle_loop

        phi_enter = 0.0
        phi_exit = 0.0
        dum = [0.d0, 0.d0, particles%fast_ion(i)%z]
        call uvw_to_xyz(dum, C_xyz)
        call circle_grid_intersect(C_xyz,e1_xyz,e2_xyz,particles%fast_ion(i)%r,phi_enter,phi_exit)
        delta_phi = phi_exit-phi_enter
        if(delta_phi.gt.0) then
            particles%fast_ion(i)%cross_grid = .True.
        else
            particles%fast_ion(i)%cross_grid = .False.
            delta_phi = 2*pi
        endif
        particles%fast_ion(i)%phi_enter = phi_enter
        particles%fast_ion(i)%delta_phi = delta_phi
        particles%fast_ion(i)%weight = weight(i)*(delta_phi/(2*pi))/beam_grid%dv

        minpos = minloc(abs(inter_grid%r - particles%fast_ion(i)%r))
        ir = minpos(1)
        minpos = minloc(abs(inter_grid%z - particles%fast_ion(i)%z))
        iz = minpos(1)
        !$OMP CRITICAL(mc_denf)
        equil%plasma(ir,iz)%denf = equil%plasma(ir,iz)%denf + weight(i) / &
                                   (2*pi*particles%fast_ion(i)%r*inter_grid%da)
        !$OMP END CRITICAL(mc_denf)
        cnt=cnt+1
    enddo particle_loop
    !$OMP END PARALLEL DO

    num = count(particles%fast_ion%cross_grid)
    if(num.le.0) then
        if(inputs%verbose.ge.0) then
            write(*,'(a)') 'READ_MC: No mc particles in beam grid'
        endif
        stop
    endif

    if(inputs%verbose.ge.1) then
        write(*,'(T2,"Distribution type: ",a)') dist_type_name
        write(*,'(T2,"Number of mc particles: ",i9)') particles%nparticle
        write(*,'(T2,"Number of orbit classes: ",i6)') particles%nclass
        write(*,*) ''
    endif

end subroutine read_mc