Reads in the interpolation grid, plasma parameters, and fields and stores the quantities in inter_grid and equil
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(HSIZE_T), dimension(4) :: dims4
integer :: impc, ic, ir, iz, iphi, it, ind(3), i
type(LocalProfiles) :: plasma
real(Float64) :: photons, smix(max_species)
real(Float64), dimension(nlevs) :: rates, denn, rates_avg
real(Float64), dimension(3) :: vi
integer :: error
integer :: n = 50
logical :: path_valid
integer, dimension(:,:,:), allocatable :: p_mask, f_mask
real(Float64), dimension(:,:,:), allocatable :: denn3d
real(Float64), dimension(:,:,:,:), allocatable :: deni
!!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/denimp", equil%plasma%denimp, 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)
! Read in ion densities
dims4(1) = n_thermal
dims4(2:4) = dims
allocate(deni(n_thermal,inter_grid%nr, inter_grid%nz, inter_grid%nphi))
call h5ltread_dataset_double_f(gid, "/plasma/deni", deni, dims4, error)
do i=1,n_thermal
equil%plasma(:,:,:)%deni(i) = deni(i,:,:,:)
enddo
deallocate(deni)
impc = 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%denimp.lt.0.0)
equil%plasma%denimp = 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
do i=1,n_thermal
where(equil%plasma%deni(i).lt.0.0)
equil%plasma%deni(i) = 0.0
endwhere
enddo
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
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.
smix = plasma%deni/sum(plasma%deni)
do i=1,n_thermal
rates_avg = 0.0
do it=1,n
rates = 0.0
rates(1) = 1.d19
call mc_halo(plasma, thermal_mass(i), vi)
call colrad(plasma, thermal_mass(i), 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(:,i) = smix(i)*denn3d(ir,iz,iphi)*(rates_avg)/sum(rates_avg)
enddo
enddo loop_over_cells
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
!!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