read_equilibrium Subroutine

public subroutine read_equilibrium()

Reads in the interpolation grid, plasma parameters, and fields and stores the quantities in inter_grid and equil

Arguments

None

Calls

proc~~read_equilibrium~~CallsGraph proc~read_equilibrium read_equilibrium h5ltread_dataset_int_f h5ltread_dataset_int_f proc~read_equilibrium->h5ltread_dataset_int_f proc~ind2sub ind2sub proc~read_equilibrium->proc~ind2sub h5gopen_f h5gopen_f proc~read_equilibrium->h5gopen_f h5ltpath_valid_f h5ltpath_valid_f proc~read_equilibrium->h5ltpath_valid_f h5gclose_f h5gclose_f proc~read_equilibrium->h5gclose_f proc~colrad colrad proc~read_equilibrium->proc~colrad proc~h5ltread_dataset_int_scalar_f h5ltread_dataset_int_scalar_f proc~read_equilibrium->proc~h5ltread_dataset_int_scalar_f h5open_f h5open_f proc~read_equilibrium->h5open_f h5ltread_dataset_double_f h5ltread_dataset_double_f proc~read_equilibrium->h5ltread_dataset_double_f interface~randn randn proc~read_equilibrium->interface~randn interface~deriv deriv proc~read_equilibrium->interface~deriv h5close_f h5close_f proc~read_equilibrium->h5close_f h5fclose_f h5fclose_f proc~read_equilibrium->h5fclose_f h5fopen_f h5fopen_f proc~read_equilibrium->h5fopen_f proc~eigen eigen proc~colrad->proc~eigen proc~get_rate_matrix get_rate_matrix proc~colrad->proc~get_rate_matrix proc~linsolve linsolve proc~colrad->proc~linsolve proc~h5ltread_dataset_int_scalar_f->h5ltread_dataset_int_f proc~hqr2 hqr2 proc~eigen->proc~hqr2 proc~balback balback proc~eigen->proc~balback proc~elmtrans elmtrans proc~eigen->proc~elmtrans proc~elmhes elmhes proc~eigen->proc~elmhes proc~balance balance proc~eigen->proc~balance interface~interpol_coeff interpol_coeff proc~get_rate_matrix->interface~interpol_coeff proc~matinv matinv proc~linsolve->proc~matinv dgetrf dgetrf proc~linsolve->dgetrf dgetrs dgetrs proc~linsolve->dgetrs proc~cyl_interpol3d_coeff cyl_interpol3D_coeff interface~interpol_coeff->proc~cyl_interpol3d_coeff proc~interpol1d_coeff interpol1D_coeff interface~interpol_coeff->proc~interpol1d_coeff proc~interpol2d_coeff_arr interpol2D_coeff_arr interface~interpol_coeff->proc~interpol2d_coeff_arr proc~interpol2d_coeff interpol2D_coeff interface~interpol_coeff->proc~interpol2d_coeff proc~cyl_interpol3d_coeff_arr cyl_interpol3D_coeff_arr interface~interpol_coeff->proc~cyl_interpol3d_coeff_arr proc~interpol1d_coeff_arr interpol1D_coeff_arr interface~interpol_coeff->proc~interpol1d_coeff_arr proc~hqrvec hqrvec proc~hqr2->proc~hqrvec proc~ludcmp ludcmp proc~matinv->proc~ludcmp proc~cyl_interpol3d_coeff->proc~interpol2d_coeff proc~outerprod outerprod proc~ludcmp->proc~outerprod proc~interpol2d_coeff_arr->proc~interpol2d_coeff proc~cyl_interpol3d_coeff_arr->proc~cyl_interpol3d_coeff proc~cyl_interpol3d_coeff_arr->proc~interpol2d_coeff proc~interpol1d_coeff_arr->proc~interpol1d_coeff

Called by

proc~~read_equilibrium~~CalledByGraph proc~read_equilibrium read_equilibrium program~fidasim fidasim program~fidasim->proc~read_equilibrium

Contents

Source Code


Source Code

subroutine read_equilibrium
    !+ Reads in the interpolation grid, plasma parameters, and fields
    !+ and stores the quantities in [[libfida:inter_grid]] and [[libfida:equil]]
    integer(HID_T) :: fid, gid
    integer(HSIZE_T), dimension(3) :: dims

    integer :: impc, ic, ir, iz, iphi, it, ind(3), i
    type(LocalProfiles) :: plasma
    real(Float64) :: photons
    real(Float64), dimension(nlevs) :: rates, denn, rates_avg
    real(Float64), dimension(3) :: vi, random3
    integer :: error
    integer :: n = 50
    logical :: path_valid

    integer, dimension(:,:,:), allocatable :: p_mask, f_mask
    real(Float64), dimension(:,:,:), allocatable :: denn3d

    !!Initialize HDF5 interface
    call h5open_f(error)

    !!Open HDF5 file
    call h5fopen_f(inputs%equilibrium_file, H5F_ACC_RDONLY_F, fid, error)

    !!Open PLASMA group
    call h5gopen_f(fid, "/plasma", gid, error)

    !!Read in interpolation grid
    call h5ltread_dataset_int_scalar_f(gid, "/plasma/nr", inter_grid%nr, error)
    call h5ltread_dataset_int_scalar_f(gid, "/plasma/nz", inter_grid%nz, error)
    call h5ltpath_valid_f(gid, "/plasma/nphi", .True., path_valid, error)
    if(path_valid) then
        call h5ltread_dataset_int_scalar_f(gid, "/plasma/nphi", inter_grid%nphi, error)
    else
        inter_grid%nphi=1
    endif

    inter_grid%dims = [inter_grid%nr, inter_grid%nz, inter_grid%nphi]

    allocate(inter_grid%r(inter_grid%nr),inter_grid%z(inter_grid%nz),inter_grid%phi(inter_grid%nphi))
    allocate(p_mask(inter_grid%nr,inter_grid%nz,inter_grid%nphi))
    allocate(f_mask(inter_grid%nr,inter_grid%nz,inter_grid%nphi))
    allocate(denn3d(inter_grid%nr,inter_grid%nz,inter_grid%nphi))

    dims = [inter_grid%nr, inter_grid%nz, inter_grid%nphi]

    call h5ltread_dataset_double_f(gid, "/plasma/r", inter_grid%r, dims(1:1), error)
    call h5ltread_dataset_double_f(gid, "/plasma/z", inter_grid%z, dims(2:2), error)
    if(path_valid) then
        call h5ltread_dataset_double_f(gid, "/plasma/phi", inter_grid%phi, dims(3:3), error)
    else
        inter_grid%phi=0.d0
    endif

    inter_grid%dr = abs(inter_grid%r(2)-inter_grid%r(1))
    inter_grid%dz = abs(inter_grid%z(2)-inter_grid%z(1))
    inter_grid%da = inter_grid%dr*inter_grid%dz
    if (inter_grid%nphi .eq. 1) then
        inter_grid%dphi = 2*pi
    else
        inter_grid%dphi = abs(inter_grid%phi(2)-inter_grid%phi(1))
    endif
    inter_grid%dv = inter_grid%dr*inter_grid%dphi*inter_grid%dz

    inter_grid%ntrack = inter_grid%nr+inter_grid%nz+inter_grid%nphi
    inter_grid%ngrid  = inter_grid%nr*inter_grid%nz*inter_grid%nphi

    if(inputs%verbose.ge.1) then
        write(*,'(a)') '---- Interpolation grid settings ----'
        write(*,'(T2,"Nr: ",i3)') inter_grid%nr
        write(*,'(T2,"Nz: ",i3)') inter_grid%nz
        if(inter_grid%nphi.gt.1) then
            write(*,'(T2,"Nphi: ",i3)') inter_grid%nphi
        endif
        write(*,'(T2,"dA: ", f5.2," [cm^2]")') inter_grid%da
        write(*,*) ''
    endif

    !!Read in plasma parameters
    allocate(equil%plasma(inter_grid%nr,inter_grid%nz,inter_grid%nphi))

    call h5ltread_dataset_double_f(gid, "/plasma/dene", equil%plasma%dene, dims, error)
    call h5ltread_dataset_double_f(gid, "/plasma/te", equil%plasma%te, dims, error)
    call h5ltread_dataset_double_f(gid, "/plasma/ti", equil%plasma%ti, dims, error)
    call h5ltread_dataset_double_f(gid, "/plasma/zeff", equil%plasma%zeff, dims, error)
    call h5ltread_dataset_double_f(gid, "/plasma/vr", equil%plasma%vr, dims, error)
    call h5ltread_dataset_double_f(gid, "/plasma/vt", equil%plasma%vt, dims, error)
    call h5ltread_dataset_double_f(gid, "/plasma/vz", equil%plasma%vz, dims, error)
    call h5ltread_dataset_int_f(gid, "/plasma/mask", p_mask, dims,error)

    impc = inputs%impurity_charge
    where(equil%plasma%zeff.lt.1.0)
        equil%plasma%zeff = 1
    endwhere

    where(equil%plasma%zeff.gt.impc)
        equil%plasma%zeff = impc
    endwhere

    where(equil%plasma%dene.lt.0.0)
        equil%plasma%dene = 0.0
    endwhere

    where(equil%plasma%te.lt.0.0)
        equil%plasma%te = 0.0
    endwhere

    where(equil%plasma%ti.lt.0.0)
        equil%plasma%ti = 0.0
    endwhere

    equil%plasma%denimp = ((equil%plasma%zeff-1.d0)/(impc*(impc-1.d0)))*equil%plasma%dene
    equil%plasma%denp = equil%plasma%dene - impc*equil%plasma%denimp

    call h5ltpath_valid_f(fid, "/plasma/denn", .True., path_valid, error)
    if(path_valid) then
        call h5ltread_dataset_double_f(gid, "/plasma/denn", denn3d, dims, error)
        where(denn3d.lt.0.0)
            denn3d = 0.0
        endwhere
    else
        if((inputs%calc_pnpa + inputs%calc_pfida + inputs%calc_cold).gt.0) then
            if(inputs%verbose.ge.0) then
                write(*,'(a)') "READ_EQUILIBRIUM: Cold neutral density was not provided"
                write(*,'(a)') "Continuing without passive calculations"
            endif
        endif
        inputs%calc_pnpa = 0
        inputs%calc_pfida = 0
        inputs%calc_cold = 0
    endif

    loop_over_cells: do ic=1, inter_grid%nr*inter_grid%nz*inter_grid%nphi
        call ind2sub(inter_grid%dims,ic,ind)
        ir = ind(1) ; iz = ind(2) ; iphi = ind(3)
        if(p_mask(ir,iz,iphi).lt.0.5) cycle loop_over_cells
        if(denn3d(ir,iz,iphi).le.0.0) cycle loop_over_cells
        plasma = equil%plasma(ir,iz,iphi)
        plasma%vrot = [plasma%vr, plasma%vt, plasma%vz]
        plasma%in_plasma = .True.
        rates_avg = 0.0
        do it=1,n
            rates = 0.0
            rates(1) = 1.d19
            call randn(random3)
            vi = plasma%vrot + sqrt(plasma%ti*0.5/(v2_to_E_per_amu*inputs%ai))*random3
            call colrad(plasma, thermal_ion, vi, 1.0d-7, rates, denn, photons)
            rates_avg = rates_avg + rates/n
        enddo
        if(sum(rates_avg).le.0.0) cycle loop_over_cells
        equil%plasma(ir,iz,iphi)%denn = denn3d(ir,iz,iphi)*(rates_avg)/sum(rates_avg)
    enddo loop_over_cells

    !!Close PLASMA group
    call h5gclose_f(gid, error)

    !!Open FIELDS group
    call h5gopen_f(fid, "/fields", gid, error)

    allocate(equil%fields(inter_grid%nr,inter_grid%nz,inter_grid%nphi))

    !!Read in electromagnetic fields
    call h5ltread_dataset_double_f(gid, "/fields/br", equil%fields%br, dims, error)
    call h5ltread_dataset_double_f(gid, "/fields/bt", equil%fields%bt, dims, error)
    call h5ltread_dataset_double_f(gid, "/fields/bz", equil%fields%bz, dims, error)
    call h5ltread_dataset_double_f(gid, "/fields/er", equil%fields%er, dims, error)
    call h5ltread_dataset_double_f(gid, "/fields/et", equil%fields%et, dims, error)
    call h5ltread_dataset_double_f(gid, "/fields/ez", equil%fields%ez, dims, error)
    call h5ltread_dataset_int_f(gid, "/fields/mask", f_mask, dims,error)

    !!Calculate B field derivatives
    call deriv(inter_grid%r, inter_grid%z, inter_grid%phi, equil%fields%br, &
        equil%fields%dbr_dr, equil%fields%dbr_dz, equil%fields%dbr_dphi)
    call deriv(inter_grid%r, inter_grid%z, inter_grid%phi, equil%fields%bt, &
        equil%fields%dbt_dr, equil%fields%dbt_dz, equil%fields%dbt_dphi)
    call deriv(inter_grid%r, inter_grid%z, inter_grid%phi, equil%fields%bz, &
        equil%fields%dbz_dr, equil%fields%dbz_dz, equil%fields%dbz_dphi)

    !!Close FIELDS group
    call h5gclose_f(gid, error)

    !!Close file
    call h5fclose_f(fid, error)

    !!Close HDF5 interface
    call h5close_f(error)

    allocate(equil%mask(inter_grid%nr,inter_grid%nz,inter_grid%nphi))
    equil%mask = 0.d0
    where ((p_mask.eq.1).and.(f_mask.eq.1)) equil%mask = 1.d0
    if (sum(equil%mask).le.0.d0) then
        if(inputs%verbose.ge.0) then
            write(*,'(a)') "READ_EQUILIBRIUM: Plasma and/or fields are not well defined anywhere"
        endif
        stop
    endif

end subroutine read_equilibrium