FIXME : Add documentation
subroutine init_lorentz_error
use le_grids, only: jend, al, ng2, nlambda, &
il_is_wfb, il_is_passing, il_is_trapped
use theta_grid, only: ntgrid, bmag
use array_utils, only: zero_array
implicit none
integer :: je, ig, il, ip, ij, im
real :: slb0, slb1, slb2, slbr, slbl
real, dimension (:), allocatable :: slb
real, dimension (:,:), allocatable :: dprod
real, dimension (:,:,:), allocatable :: dlcoef, d2lcoef
allocate(slb(2*nlambda))
allocate (dprod(nlambda,5))
allocate (dlcoef(-ntgrid:ntgrid,nlambda,5))
allocate (d2lcoef(-ntgrid:ntgrid,nlambda,5))
allocate (dtot(-ntgrid:ntgrid,nlambda,5))
allocate (fdf(-ntgrid:ntgrid,nlambda), fdb(-ntgrid:ntgrid,nlambda))
dlcoef = 1.0; call zero_array(d2lcoef); call zero_array(dtot)
call zero_array(fdf); call zero_array(fdb); slb = 0.0
do ig=-ntgrid,ntgrid
je = jend(ig)
if (il_is_passing(je) .or. il_is_wfb(je)) then ! no trapped particles
! calculation of xi and finite difference coefficients for non-boundary points
do il=2,ng2-1
slb(il) = safe_sqrt(1.0-al(il)*bmag(ig)) ! xi_{j}
slb2 = safe_sqrt(1.0-al(il+1)*bmag(ig)) ! xi_{j+1}
slb1 = slb(il)
slb0 = safe_sqrt(1.0-al(il-1)*bmag(ig)) ! xi_{j-1}
slbr = (slb2+slb1)*0.5 ! xi_{j+1/2}
slbl = (slb1+slb0)*0.5 ! xi_{j-1/2}
! finite difference coefficients
fdf(ig,il) = (1.0 - slbr*slbr)/(slbr - slbl)/(slb2 - slb1)
fdb(ig,il) = (1.0 - slbl*slbl)/(slbr - slbl)/(slb1 - slb0)
end do
! boundary at xi = 1
slb(1) = safe_sqrt(1.0-al(1)*bmag(ig))
slb0 = 1.0
slb1 = slb(1)
slb2 = slb(2)
slbl = (slb1 + slb0)/2.0
slbr = (slb1 + slb2)/2.0
! derivative of [(1-xi**2)*df/dxi] at xi_{j=1} is centered, with upper xi=1 and
! lower xi = xi_{j+1/2}
fdf(ig,1) = (-1.0-slbr)/(slb2-slb1)
fdb(ig,1) = 0.0
! boundary at xi = 0
il = ng2
slb(il) = safe_sqrt(1.0 - al(il)*bmag(ig))
slb0 = safe_sqrt(1.0 - bmag(ig)*al(il-1))
slb1 = slb(il)
slb2 = -slb1
slbl = (slb1 + slb0)/2.0
slbr = (slb1 + slb2)/2.0
fdf(ig,il) = (1.0 - slbr*slbr)/(slbr-slbl)/(slb2-slb1)
fdb(ig,il) = (1.0 - slbl*slbl)/(slbr-slbl)/(slb1-slb0)
slb(ng2+1:) = -slb(ng2:1:-1)
else ! run with trapped particles
do il=2,je-1
slb(il) = safe_sqrt(1.0-al(il)*bmag(ig))
slb2 = safe_sqrt(1.0-al(il+1)*bmag(ig))
slb1 = slb(il)
slb0 = safe_sqrt(1.0-al(il-1)*bmag(ig))
slbr = (slb2+slb1)*0.5
slbl = (slb1+slb0)*0.5
fdf(ig,il) = (1.0 - slbr*slbr)/(slbr - slbl)/(slb2 - slb1)
fdb(ig,il) = (1.0 - slbl*slbl)/(slbr - slbl)/(slb1 - slb0)
end do
! boundary at xi = 1
slb(1) = safe_sqrt(1.0-bmag(ig)*al(1))
slb0 = 1.0
slb1 = slb(1)
slb2 = slb(2)
slbr = (slb1 + slb2)/2.0
fdf(ig,1) = (-1.0 - slbr)/(slb2-slb1)
fdb(ig,1) = 0.0
! boundary at xi = 0
il = je
slb(il) = safe_sqrt(1.0-bmag(ig)*al(il))
slb0 = slb(je-1)
slb1 = 0.
slb2 = -slb0
slbl = (slb1 + slb0)/2.0
fdf(ig,il) = (1.0 - slbl*slbl)/slb0/slb0
fdb(ig,il) = fdf(ig,il)
slb(je+1:2*je-1) = -slb(je-1:1:-1)
end if
! compute coefficients (dlcoef) multipyling first derivative of h
do il=3,ng2
do ip=il-2,il+2
if (il == ip) then
dlcoef(ig,il,ip-il+3) = 0.0
do ij=il-2,il+2
if (ij /= ip) dlcoef(ig,il,ip-il+3) = dlcoef(ig,il,ip-il+3) + 1/(slb(il)-slb(ij))
end do
else
do ij=il-2,il+2
if (ij /= ip .and. ij /= il) then
dlcoef(ig,il,ip-il+3) = dlcoef(ig,il,ip-il+3)*(slb(il)-slb(ij))/(slb(ip)-slb(ij))
end if
end do
dlcoef(ig,il,ip-il+3) = dlcoef(ig,il,ip-il+3)/(slb(ip)-slb(il))
end if
dlcoef(ig,il,ip-il+3) = -2.0*slb(il)*dlcoef(ig,il,ip-il+3)
end do
end do
il = 1
do ip=il,il+2
if (il == ip) then
dlcoef(ig,il,ip) = 0.0
do ij=il,il+2
if (ij /= ip) dlcoef(ig,il,ip) = dlcoef(ig,il,ip) + 1./(slb(il)-slb(ij))
end do
else
do ij=il,il+2
if (ij /= ip .and. ij /= il) then
dlcoef(ig,il,ip) = dlcoef(ig,il,ip)*(slb(il)-slb(ij))/(slb(ip)-slb(ij))
end if
end do
dlcoef(ig,il,ip) = dlcoef(ig,il,ip)/(slb(ip)-slb(il))
end if
dlcoef(ig,il,ip) = -2.0*slb(il)*dlcoef(ig,il,ip)
end do
il = 2
do ip=il-1,il+1
if (il == ip) then
dlcoef(ig,il,ip-il+2) = 0.0
do ij=il-1,il+1
if (ij /= ip) dlcoef(ig,il,ip-il+2) = dlcoef(ig,il,ip-il+2) + 1/(slb(il)-slb(ij))
end do
else
do ij=il-1,il+1
if (ij /= ip .and. ij /= il) then
dlcoef(ig,il,ip-il+2) = dlcoef(ig,il,ip-il+2)*(slb(il)-slb(ij))/(slb(ip)-slb(ij))
end if
end do
dlcoef(ig,il,ip-il+2) = dlcoef(ig,il,ip-il+2)/(slb(ip)-slb(il))
end if
dlcoef(ig,il,ip-il+2) = -2.0*slb(il)*dlcoef(ig,il,ip-il+2)
end do
dprod = 2.0
! compute coefficients (d2lcoef) multiplying second derivative of h
do il=3,ng2
do ip=il-2,il+2
if (il == ip) then
do ij=il-2,il+2
if (ij /= ip) then
do im=il-2,il+2
if (im /= ip .and. im /= ij) d2lcoef(ig,il,ip-il+3) = &
d2lcoef(ig,il,ip-il+3) + 1./((slb(il)-slb(im))*(slb(il)-slb(ij)))
end do
end if
end do
else
do ij=il-2,il+2
if (ij /= il .and. ij /= ip) then
dprod(il,ip-il+3) = dprod(il,ip-il+3)*(slb(il)-slb(ij))/(slb(ip)-slb(ij))
end if
end do
do ij=il-2,il+2
if (ij /= ip .and. ij /= il) then
d2lcoef(ig,il,ip-il+3) = d2lcoef(ig,il,ip-il+3) + 1./(slb(il)-slb(ij))
end if
end do
d2lcoef(ig,il,ip-il+3) = dprod(il,ip-il+3) &
*d2lcoef(ig,il,ip-il+3)/(slb(ip)-slb(il))
end if
d2lcoef(ig,il,ip-il+3) = (1.0-slb(il)**2)*d2lcoef(ig,il,ip-il+3)
end do
end do
il = 1
do ip=il,il+2
if (il == ip) then
do ij=il,il+2
if (ij /= ip) then
do im=il,il+2
if (im /= ip .and. im /= ij) d2lcoef(ig,il,ip) = d2lcoef(ig,il,ip) + 1./((slb(il)-slb(im))*(slb(il)-slb(ij)))
end do
end if
end do
else
do ij=il,il+2
if (ij /= il .and. ij /= ip) then
dprod(il,ip) = dprod(il,ip)*(slb(il)-slb(ij))/(slb(ip)-slb(ij))
end if
end do
do ij=il,il+2
if (ij /= ip .and. ij /= il) then
d2lcoef(ig,il,ip) = d2lcoef(ig,il,ip) + 1./(slb(il)-slb(ij))
end if
end do
d2lcoef(ig,il,ip) = dprod(il,ip)*d2lcoef(ig,il,ip)/(slb(ip)-slb(il))
end if
d2lcoef(ig,il,ip) = (1.0-slb(il)**2)*d2lcoef(ig,il,ip)
end do
il = 2
do ip=il-1,il+1
if (il == ip) then
do ij=il-1,il+1
if (ij /= ip) then
do im=il-1,il+1
if (im /= ip .and. im /= ij) d2lcoef(ig,il,ip-il+2) &
= d2lcoef(ig,il,ip-il+2) + 1./((slb(il)-slb(im))*(slb(il)-slb(ij)))
end do
end if
end do
else
do ij=il-1,il+1
if (ij /= il .and. ij /= ip) then
dprod(il,ip-il+2) = dprod(il,ip-il+2)*(slb(il)-slb(ij))/(slb(ip)-slb(ij))
end if
end do
do ij=il-1,il+1
if (ij /= ip .and. ij /= il) then
d2lcoef(ig,il,ip-il+2) = d2lcoef(ig,il,ip-il+2) + 1./(slb(il)-slb(ij))
end if
end do
d2lcoef(ig,il,ip-il+2) = dprod(il,ip-il+2)*d2lcoef(ig,il,ip-il+2)/(slb(ip)-slb(il))
end if
d2lcoef(ig,il,ip-il+2) = (1.0-slb(il)**2)*d2lcoef(ig,il,ip-il+2)
end do
if (il_is_trapped(je)) then ! have to handle trapped particles
do il=ng2+1,je
do ip=il-2,il+2
if (il == ip) then
dlcoef(ig,il,ip-il+3) = 0.0
do ij=il-2,il+2
if (ij /= ip) dlcoef(ig,il,ip-il+3) = dlcoef(ig,il,ip-il+3) + 1/(slb(il)-slb(ij))
end do
else
do ij=il-2,il+2
if (ij /= ip .and. ij /= il) then
dlcoef(ig,il,ip-il+3) = dlcoef(ig,il,ip-il+3)*(slb(il)-slb(ij))/(slb(ip)-slb(ij))
end if
end do
dlcoef(ig,il,ip-il+3) = dlcoef(ig,il,ip-il+3)/(slb(ip)-slb(il))
end if
dlcoef(ig,il,ip-il+3) = -2.0*slb(il)*dlcoef(ig,il,ip-il+3)
end do
end do
do il=ng2+1,je
do ip=il-2,il+2
if (il == ip) then
do ij=il-2,il+2
if (ij /= ip) then
do im=il-2,il+2
if (im /= ip .and. im /= ij) d2lcoef(ig,il,ip-il+3) = &
d2lcoef(ig,il,ip-il+3) + 1./((slb(il)-slb(im))*(slb(il)-slb(ij)))
end do
end if
end do
else
do ij=il-2,il+2
if (ij /= il .and. ij /= ip) then
dprod(il,ip-il+3) = dprod(il,ip-il+3)*(slb(il)-slb(ij))/(slb(ip)-slb(ij))
end if
end do
do ij=il-2,il+2
if (ij /= ip .and. ij /= il) then
d2lcoef(ig,il,ip-il+3) = d2lcoef(ig,il,ip-il+3) + 1./(slb(il)-slb(ij))
end if
end do
d2lcoef(ig,il,ip-il+3) = dprod(il,ip-il+3) &
*d2lcoef(ig,il,ip-il+3)/(slb(ip)-slb(il))
end if
d2lcoef(ig,il,ip-il+3) = (1.0-slb(il)**2)*d2lcoef(ig,il,ip-il+3)
end do
end do
end if
end do
dtot = dlcoef + d2lcoef
deallocate (slb, dprod, dlcoef, d2lcoef)
end subroutine init_lorentz_error