read_f Subroutine

public subroutine read_f(fid, error)

Reads in the fast-ion distribution function and stores the quantities in fbm

Arguments

Type IntentOptional AttributesName
integer(kind=HID_T), intent(inout) :: fid

HDF5 file ID

integer, intent(out) :: error

Error code


Calls

proc~~read_f~~CallsGraph proc~read_f read_f proc~h5ltread_dataset_int_scalar_f h5ltread_dataset_int_scalar_f proc~read_f->proc~h5ltread_dataset_int_scalar_f h5ltpath_valid_f h5ltpath_valid_f proc~read_f->h5ltpath_valid_f h5ltread_dataset_int_f h5ltread_dataset_int_f proc~h5ltread_dataset_int_scalar_f->h5ltread_dataset_int_f

Called by

proc~~read_f~~CalledByGraph proc~read_f read_f proc~read_distribution read_distribution proc~read_distribution->proc~read_f program~fidasim fidasim program~fidasim->proc~read_distribution

Contents

Source Code


Source Code

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(5) :: dims
    real(Float64) :: denp_tot
    integer :: ir
    logical :: path_valid

    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)
    call h5ltpath_valid_f(fid, "/nphi", .True., path_valid, error)
    if(path_valid) then
        call h5ltread_dataset_int_scalar_f(fid,"/nphi", fbm%nphi, error)
    else
        fbm%nphi=1
    endif

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

    if (fbm%nphi .eq. 1) then
        allocate(fbm%energy(fbm%nenergy), fbm%pitch(fbm%npitch), fbm%r(fbm%nr), fbm%z(fbm%nz), fbm%phi(1))
        allocate(fbm%denf(fbm%nr, fbm%nz,1))
        allocate(fbm%f(fbm%nenergy, fbm%npitch, fbm%nr, fbm%nz,1))
    else
        allocate(fbm%energy(fbm%nenergy), fbm%pitch(fbm%npitch), fbm%r(fbm%nr), fbm%z(fbm%nz), fbm%phi(fbm%nphi))
        allocate(fbm%denf(fbm%nr, fbm%nz, fbm%nphi))
        allocate(fbm%f(fbm%nenergy, fbm%npitch, fbm%nr, fbm%nz, fbm%nphi))
    endif

    if (fbm%nphi .eq. 1) then
        dims = [fbm%nenergy, fbm%npitch, fbm%nr, fbm%nz, 1]
        call h5ltread_dataset_double_f(fid, "/denf",fbm%denf, dims(3:4), error)
        call h5ltread_dataset_double_f(fid, "/f", fbm%f, dims(1:4), error)
    else
        dims = [fbm%nenergy, fbm%npitch, fbm%nr, fbm%nz, fbm%nphi]
        call h5ltread_dataset_double_f(fid, "/denf",fbm%denf, dims(3:5), error)
        call h5ltread_dataset_double_f(fid, "/f", fbm%f, dims(1:5), error)
    endif
    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 h5ltpath_valid_f(fid, "/phi", .True., path_valid, error)
    if(path_valid) then
        call h5ltread_dataset_double_f(fid, "/phi", fbm%phi, dims(5:5), error)
    else
        fbm%phi=0.d0
    endif

    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))
    if (fbm%nphi .eq. 1) then
        fbm%dphi = 2*pi
    else
        fbm%dphi = abs(fbm%phi(2)-fbm%phi(1))
    endif

    fbm%emin = minval(fbm%energy,1)
    fbm%emax = maxval(fbm%energy,1)
    fbm%e_range = fbm%emax - fbm%emin
    fbm%pmin = minval(fbm%pitch,1)
    fbm%pmax = maxval(fbm%pitch,1)
    fbm%p_range = fbm%pmax - fbm%pmin
    fbm%rmin = minval(fbm%r,1)
    fbm%rmax = maxval(fbm%r,1)
    fbm%r_range = fbm%rmax - fbm%rmin
    fbm%zmin = minval(fbm%z,1)
    fbm%zmax = maxval(fbm%z,1)
    fbm%z_range = fbm%zmax - fbm%zmin
    fbm%phimin = minval(fbm%phi,1)
    fbm%phimax = maxval(fbm%phi,1)
    fbm%phi_range = fbm%phimax - fbm%phimin

    denp_tot = 0.0
    do ir=1,fbm%nr
        fbm%n_tot = fbm%n_tot + fbm%dphi*fbm%dr*fbm%dz*sum(fbm%denf(ir,:,:))*fbm%r(ir)
        denp_tot = denp_tot + fbm%dphi*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
        if(fbm%nphi.gt.1) then
            write(*,'(T2,"Distribution type: ",a)') "Non-axisymmetric Fast-ion Density Function F(energy,pitch,R,Z,Phi)"
        else
            write(*,'(T2,"Distribution type: ",a)') "Axisymmetric Fast-ion Density Function F(energy,pitch,R,Z)"
        endif
        write(*,'(T2,"Nenergy = ",i3)') fbm%nenergy
        write(*,'(T2,"Npitch  = ",i3)') fbm%npitch
        write(*,'(T2,"Nr  = ",i3)') fbm%nr
        write(*,'(T2,"Nz  = ",i3)') fbm%nz
        if(fbm%nphi.gt.1) then
            write(*,'(T2,"Nphi  = ",i3)') fbm%nphi
        endif
        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,"R  range = [",f6.2,",",f6.2,"]")') fbm%rmin,fbm%rmax
        write(*,'(T2,"Z  range = [",f7.2,",",f6.2,"]")') fbm%zmin,fbm%zmax
        if(fbm%nphi.gt.1) then
            write(*,'(T2,"Phi  range = [",f5.2,",",f5.2,"]")') fbm%phimin,fbm%phimax
        endif
        write(*,'(T2,"Ntotal = ",ES10.3)') fbm%n_tot
        write(*,*) ''
    endif

end subroutine read_f