read_f Subroutine

subroutine read_f(fid, error)
    !+ Reads in the fast-ion distribution function and stores the quantities in [[libfida:fbm]]
    integer(HID_T), intent(inout) :: fid
        !+ HDF5 file ID
    integer, intent(out)          :: error
        !+ Error code

    integer(HSIZE_T), dimension(4) :: dims
    real(Float64) :: dummy(1), denp_tot
    integer :: ir

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

    call h5ltread_dataset_int_scalar_f(fid,"/nenergy", fbm%nenergy, error)
    call h5ltread_dataset_int_scalar_f(fid,"/npitch", fbm%npitch, error)
    call h5ltread_dataset_int_scalar_f(fid,"/nr", fbm%nr, error)
    call h5ltread_dataset_int_scalar_f(fid,"/nz", fbm%nz, error)

    if((fbm%nr.ne.inter_grid%nr).or.(fbm%nz.ne.inter_grid%nz)) then
        if(inputs%verbose.ge.0) then
            write(*,'(a)') "READ_F: Distribution file has incompatable grid dimensions"
        endif
        stop
    endif

    allocate(fbm%energy(fbm%nenergy), fbm%pitch(fbm%npitch), fbm%r(fbm%nr), fbm%z(fbm%nz))
    allocate(fbm%denf(fbm%nr, fbm%nz))
    allocate(fbm%f(fbm%nenergy, fbm%npitch, fbm%nr, fbm%nz))

    dims = [fbm%nenergy, fbm%npitch, fbm%nr, fbm%nz]
    call h5ltread_dataset_double_f(fid, "/energy", fbm%energy, dims(1:1), error)
    call h5ltread_dataset_double_f(fid, "/pitch", fbm%pitch, dims(2:2), error)
    call h5ltread_dataset_double_f(fid, "/r", fbm%r, dims(3:3), error)
    call h5ltread_dataset_double_f(fid, "/z", fbm%z, dims(4:4), error)
    call h5ltread_dataset_double_f(fid, "/denf",fbm%denf, dims(3:4), error)
    call h5ltread_dataset_double_f(fid, "/f", fbm%f, dims, error)
    equil%plasma%denf = fbm%denf

    fbm%dE = abs(fbm%energy(2) - fbm%energy(1))
    fbm%dp = abs(fbm%pitch(2) - fbm%pitch(1))
    fbm%dr = abs(fbm%r(2) - fbm%r(1))
    fbm%dz = abs(fbm%z(2) - fbm%z(1))

    dummy = minval(fbm%energy)
    fbm%emin = dummy(1)
    dummy = maxval(fbm%energy)
    fbm%emax = dummy(1)
    fbm%e_range = fbm%emax - fbm%emin
    dummy = minval(fbm%pitch)
    fbm%pmin = dummy(1)
    dummy = maxval(fbm%pitch)
    fbm%pmax = dummy(1)
    fbm%p_range = fbm%pmax - fbm%pmin

    denp_tot = 0.0
    do ir=1,fbm%nr
        fbm%n_tot = fbm%n_tot + 2*pi*fbm%dr*fbm%dz*sum(fbm%denf(ir,:))*fbm%r(ir)
        denp_tot = denp_tot + 2*pi*fbm%dr*fbm%dz*sum(equil%plasma(ir,:)%denp)*fbm%r(ir)
    enddo

    if(fbm%n_tot.ge.denp_tot) then
        if(inputs%verbose.ge.0) then
            write(*,'(a," (",ES10.3," >=",ES10.3,")")') &
                "READ_F: The total of number of fast ions exceeded the total number of thermal ions.", &
                 fbm%n_tot, denp_tot
            write(*,'(a)') "This is usually caused by zeff being incorrect."
        endif
        stop
    endif

    if(inputs%verbose.ge.1) then
        write(*,'(T2,"Distribution type: ",a)') "Fast-ion Density Function F(energy,pitch,R,Z)"
        write(*,'(T2,"Nenergy = ",i3)') fbm%nenergy
        write(*,'(T2,"Npitch  = ",i3)') fbm%npitch
        write(*,'(T2,"Energy range = [",f5.2,",",f6.2,"]")') fbm%emin,fbm%emax
        write(*,'(T2,"Pitch  range = [",f5.2,",",f5.2,"]")') fbm%pmin,fbm%pmax
        write(*,'(T2,"Ntotal = ",ES10.3)') fbm%n_tot
        write(*,*) ''
    endif

end subroutine read_f