This subroutine calculates f0 on the grid from an external input file. The grid on the input file can differ from that of gs2. A cubic spline is used to interpolate between the two. The user can either specify df0/dE or it can be estimated internally. The radial derivative of f0 as a function of energy should also be given.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | is | |||
| real, | intent(in), | dimension(:,:) | :: | egrid |
subroutine calculate_f0_arrays_tabulated(is,egrid)
!
! The egrid is already given by get_legendre_grids.
!
! The input file can take two forms, as controlled by the
! control parameter mode, the first integer read from the file
! mode = 1: E points in first columns, F0(E) in second. To calculate
! generalized temperature, a cubic spline is used to
! interpolate and the slope is taken from that.
! = 2: First column is E, second is F0(E), the third is (1/F0)*dF0/dE
! = 3: E, F0(E), -(1/F0)*dF0/drho
! = 4: E, F0(E),dF0/dE, -(1/F0)*dF0/drho
use mp, only: broadcast
use constants, only: pi
use file_utils, only: run_name
use splines, only: spline, new_spline, delete_spline, splint, dsplint
implicit none
integer, intent(in) :: is
real, dimension(:,:),intent(in):: egrid
type(spline) :: spl
integer:: f0in_unit, mode, ie, numdat, negrid
logical:: acceptable_spline
real:: moment0, sigma, sigma_fac, df0dE
real, dimension(:), allocatable:: f0_values_dat, df0dE_dat, egrid_dat, df0drho_dat
! Open file and read column option
open(newunit = f0in_unit,file=trim(run_name)//'.f0in',status='old',action='read')
read(f0in_unit,*) mode
read(f0in_unit,*) numdat
allocate(f0_values_dat(numdat))
allocate(df0dE_dat(numdat))
allocate(df0drho_dat(numdat))
allocate(egrid_dat(numdat))
negrid = size(egrid(:,1))
sigma = 1.0
sigma_fac = 2.0
acceptable_spline = .false.
if (mode == 1) then
! Read f0 values
do ie = 1,numdat
read(f0in_unit,*) egrid_dat(ie), f0_values_dat(ie)
end do
close(f0in_unit)
! Perform cubic spline to get F0 and its slope at grid points
! If F0 ends up being negative or zero anywhere on the grid, double the tension
! factor and try again
do while (.not. acceptable_spline)
! Generate spline parameters
spl = new_spline(egrid_dat, f0_values_dat, tension = sigma)
do ie = 1,negrid
! Interpolate to get f0 at grid points
f0_values(ie,is) = splint(egrid(ie, is), spl)
! Calculate d/dE lnF0 to get generalised temperature
df0dE = dsplint(egrid(ie,is), spl) / (f0_values(ie,is) * spec(is)%temp)
nonmaxw_corr(ie,is) = -spec(is)%temp*df0dE
end do
call delete_spline(spl)
if (minval(f0_values(:,is)) > epsilon(0.0)) then
acceptable_spline = .true.
else
sigma = sigma * sigma_fac
end if
end do
dlnf0drho(:,is) = - spec(is)%fprim
else if (mode == 2) then
! Read both f0 and df0/dE
do ie = 1,numdat
read(f0in_unit,*) egrid_dat(ie), f0_values_dat(ie), df0dE_dat(ie)
end do
close(f0in_unit)
do while (.not. acceptable_spline)
! Generate spline parameters for f0
spl = new_spline(egrid_dat, f0_values_dat, tension = sigma)
do ie = 1,negrid
! Interpolate to get F0 at grid points
f0_values(ie,is) = splint(egrid(ie, is), spl)
end do
call delete_spline(spl)
if (minval(f0_values(:,is)) > epsilon(0.0)) then
acceptable_spline = .true.
else
sigma = sigma * sigma_fac
cycle
end if
! Generate spline parameters for df0/dE
spl = new_spline(egrid_dat, df0dE_dat, tension = sigma)
do ie = 1,negrid
! Interpolate to get f0 at grid points
df0dE = splint(egrid(ie,is), spl) / spec(is)%temp
nonmaxw_corr(ie,is) = -spec(is)%temp*df0dE
end do
call delete_spline(spl)
end do
dlnf0drho(:,is) = - spec(is)%fprim
else if (mode == 3) then
! Read f0 and df0/drho
do ie = 1,numdat
read(f0in_unit,*) egrid_dat(ie), f0_values_dat(ie), df0drho_dat(ie)
end do
close(f0in_unit)
do while (.not. acceptable_spline)
! Generate spline parameters for f0
spl = new_spline(egrid_dat, f0_values_dat, tension = sigma)
do ie = 1,negrid
! Interpolate to get F0 at grid points
f0_values(ie,is) = splint(egrid(ie,is), spl)
! Calculate d/dE lnF0 to get generalised temperature
df0dE = dsplint(egrid(ie,is), spl) / (f0_values(ie,is) * spec(is)%temp)
nonmaxw_corr(ie,is) = -spec(is)%temp*df0dE
end do
call delete_spline(spl)
if (minval(f0_values(:,is)) > epsilon(0.0)) then
acceptable_spline = .true.
else
sigma = sigma * sigma_fac
cycle
end if
! Generate spline parameters for df0/drho
spl = new_spline(egrid_dat, df0drho_dat, tension = sigma)
do ie = 1,negrid
! Interpolate to get f0prim at grid points
dlnf0drho(ie,is) = splint(egrid(ie,is), spl)
end do
call delete_spline(spl)
end do
else if (mode == 4) then
! Read f0, df0/dE, df0/drho
do ie = 1,numdat
read(f0in_unit,*) egrid_dat(ie), f0_values_dat(ie), df0dE_dat(ie), df0drho_dat(ie)
end do
close(f0in_unit)
do while (.not. acceptable_spline)
! Generate spline parameters for f0
spl = new_spline(egrid_dat, f0_values_dat, tension = sigma)
do ie = 1,negrid
! Interpolate to get F0 at grid points
f0_values(ie,is) = splint(egrid(ie,is), spl)
end do
call delete_spline(spl)
if (minval(f0_values(:,is)) > epsilon(0.0)) then
acceptable_spline = .true.
else
sigma = sigma * sigma_fac
cycle
end if
! Generate spline parameters for df0/dE
spl = new_spline(egrid_dat, df0dE_dat, tension = sigma)
do ie = 1,negrid
! Interpolate to get f0 at grid points
df0dE = splint(egrid(ie,is), spl) / spec(is)%temp
nonmaxw_corr(ie,is) = -spec(is)%temp*df0dE
end do
call delete_spline(spl)
! Generate spline parameters for df0/drho
spl = new_spline(egrid_dat, df0drho_dat, tension = sigma)
do ie = 1,negrid
! Interpolate to get f0 at grid points
dlnf0drho(ie,is) = splint(egrid(ie,is), spl)
end do
call delete_spline(spl)
end do
dlnf0drho(:,is) = - spec(is)%fprim
else
write(*,*) "ERROR. First line in f0 input file should be mode"
stop 1
end if
! Now calculate moments of F0 from trapezoidal rule.
! Should probably do this by quadrature,
! but le_grids is not defined yet
moment0 = 0.0
do ie = 1,numdat-1
moment0 = moment0 + 0.5*(egrid_dat(ie+1)-egrid_dat(ie)) * &
( sqrt(egrid_dat(ie))*f0_values_dat(ie) + &
sqrt(egrid_dat(ie+1))*f0_values_dat(ie+1) )
end do
moment0 = moment0*4.0*pi
deallocate(egrid_dat,f0_values_dat,df0dE_dat,df0drho_dat)
end subroutine calculate_f0_arrays_tabulated