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(2) :: dims
integer :: impc
integer :: error
integer, dimension(:,:), allocatable :: p_mask, f_mask
!!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)
allocate(inter_grid%r(inter_grid%nr),inter_grid%z(inter_grid%nz))
allocate(inter_grid%r2d(inter_grid%nr,inter_grid%nz))
allocate(inter_grid%z2d(inter_grid%nr,inter_grid%nz))
allocate(p_mask(inter_grid%nr,inter_grid%nz))
allocate(f_mask(inter_grid%nr,inter_grid%nz))
dims = [inter_grid%nr, inter_grid%nz]
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)
call h5ltread_dataset_double_f(gid, "/plasma/r2d", inter_grid%r2d, dims, error)
call h5ltread_dataset_double_f(gid, "/plasma/z2d", inter_grid%z2d, dims, error)
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(inputs%verbose.ge.1) then
write(*,'(a)') '---- Interpolation grid settings ----'
write(*,'(T2,"Nr: ",i3)') inter_grid%nr
write(*,'(T2,"Nz: ",i3)') inter_grid%nz
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))
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
!!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))
!!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, equil%fields%br, equil%fields%dbr_dr, equil%fields%dbr_dz)
call deriv(inter_grid%r, inter_grid%z, equil%fields%bt, equil%fields%dbt_dr, equil%fields%dbt_dz)
call deriv(inter_grid%r, inter_grid%z, equil%fields%bz, equil%fields%dbz_dr, equil%fields%dbz_dz)
!!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))
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
write(*,'(a)') "READ_EQUILIBRIUM: Plasma and/or fields are not well defined anywhere"
stop
endif
end subroutine read_equilibrium