Subroutine to advance GC variables
Comment this section further with evolution equations, numerical methods, and descriptions of both.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(PARTICLES), | intent(inout) | :: | vars | |||
| integer, | intent(in) | :: | pp | |||
| integer(kind=ip), | intent(in) | :: | tt | time iterator. |
||
| type(KORC_PARAMS), | intent(inout) | :: | params | Core KORC simulation parameters. |
||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | Y_R | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | Y_PHI | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | Y_Z | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | V_PLL | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | V_MU | ||
| real(kind=rp), | intent(in) | :: | q_cache | |||
| real(kind=rp), | intent(in) | :: | m_cache | |||
| integer(kind=is), | intent(inout), | DIMENSION(params%pchunk) | :: | flagCon | ||
| integer(kind=is), | intent(inout), | DIMENSION(params%pchunk) | :: | flagCol | ||
| type(FIELDS), | intent(in) | :: | F | |||
| type(PROFILES), | intent(in) | :: | P | |||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | B_R | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | B_PHI | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | B_Z | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | E_R | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | E_PHI | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | E_Z | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | PSIp | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | curlb_R | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | curlb_PHI | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | curlb_Z | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | gradB_R | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | gradB_PHI | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | gradB_Z | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | ne | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | ni | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | Te | ||
| real(kind=rp), | intent(out), | DIMENSION(params%pchunk) | :: | Zeff | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk,params%num_impurity_species) | :: | nimp | ||
| type(C_PTR), | intent(inout), | DIMENSION(params%pchunk) | :: | hint |
subroutine advance_GCinterp_fio_vars(vars,pp,tt,params,Y_R,Y_PHI,Y_Z, &
V_PLL,V_MU,q_cache,m_cache,flagCon,flagCol,F,P,B_R,B_PHI,B_Z, &
E_R,E_PHI,E_Z,PSIp,curlb_R,curlb_PHI,curlb_Z, &
gradB_R,gradB_PHI,gradB_Z,ne,ni,Te,Zeff,nimp,hint)
USE omp_lib
IMPLICIT NONE
!! @note Subroutine to advance GC variables \(({\bf X},p_\parallel)\)
!! @endnote
!! Comment this section further with evolution equations, numerical
!! methods, and descriptions of both.
TYPE(KORC_PARAMS), INTENT(INOUT) :: params
!! Core KORC simulation parameters.
TYPE(PARTICLES), INTENT(INOUT) :: vars
TYPE(PROFILES), INTENT(IN) :: P
TYPE(FIELDS), INTENT(IN) :: F
REAL(rp) :: dt
!! Time step used in the leapfrog step (\(\Delta t\)).
INTEGER :: cc,pchunk
!! Chunk iterator.
INTEGER(ip),intent(in) :: tt
!! time iterator.
INTEGER,intent(in) :: pp
REAL(rp),DIMENSION(params%pchunk) :: Bmag
REAL(rp) :: a1 = 1./5._rp
REAL(rp) :: a21 = 3./40._rp,a22=9./40._rp
REAL(rp) :: a31 = 3./10._rp,a32=-9./10._rp,a33=6./5._rp
REAL(rp) :: a41 = -11./54._rp,a42=5./2._rp,a43=-70./27._rp,a44=35./27._rp
REAL(rp) :: a51 = 1631./55296._rp,a52=175./512._rp,a53=575./13824._rp,a54=44275./110592._rp,a55=253./4096._rp
REAL(rp) :: b1=37./378._rp,b2=0._rp,b3=250./621._rp,b4=125./594._rp,b5=0._rp,b6=512./1771._rp
REAL(rp),DIMENSION(params%pchunk) :: k1_R,k1_PHI,k1_Z,k1_PLL,k1_MU
REAL(rp),DIMENSION(params%pchunk) :: k2_R,k2_PHI,k2_Z,k2_PLL,k2_MU
REAL(rp),DIMENSION(params%pchunk) :: k3_R,k3_PHI,k3_Z,k3_PLL,k3_MU
REAL(rp),DIMENSION(params%pchunk) :: k4_R,k4_PHI,k4_Z,k4_PLL,k4_MU
REAL(rp),DIMENSION(params%pchunk) :: k5_R,k5_PHI,k5_Z,k5_PLL,k5_MU
REAL(rp),DIMENSION(params%pchunk) :: k6_R,k6_PHI,k6_Z,k6_PLL,k6_MU
REAL(rp),DIMENSION(params%pchunk) :: Y0_R,Y0_PHI,Y0_Z
REAL(rp),DIMENSION(params%pchunk),INTENT(INOUT) :: Y_R,Y_PHI,Y_Z
REAL(rp),DIMENSION(params%pchunk),INTENT(INOUT) :: B_R,B_PHI,B_Z
REAL(rp),DIMENSION(params%pchunk),INTENT(INOUT) :: E_PHI,E_R,E_Z
REAL(rp),DIMENSION(params%pchunk),INTENT(INOUT) :: PSIp
REAL(rp),DIMENSION(params%pchunk),INTENT(INOUT) :: ne,Te,ni
REAL(rp),DIMENSION(params%pchunk),INTENT(OUT) :: Zeff
REAL(rp),DIMENSION(params%pchunk,params%num_impurity_species)&
&,INTENT(INOUT) :: nimp
REAL(rp),DIMENSION(params%pchunk),INTENT(INOUT) :: curlb_R&
&,curlb_PHI,curlb_Z
REAL(rp),DIMENSION(params%pchunk),INTENT(INOUT) :: gradB_R,gradB_PHI,gradB_Z
REAL(rp),DIMENSION(params%pchunk),INTENT(INOUT) :: V_PLL,V_MU
REAL(rp),DIMENSION(params%pchunk) :: RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU
REAL(rp),DIMENSION(params%pchunk) :: V0_PLL,V0_MU
INTEGER(is),DIMENSION(params%pchunk),intent(INOUT) :: flagCon,flagCol
REAL(rp),intent(IN) :: q_cache,m_cache
TYPE(C_PTR), DIMENSION(params%pchunk), INTENT(INOUT) :: hint
INTEGER :: thread_num
thread_num = OMP_GET_THREAD_NUM()
dt=params%dt
pchunk=params%pchunk
!$OMP SIMD
! !$OMP& aligned(Y0_R,Y0_PHI,Y0_Z,V0_PLL,V0_MU,Y_R,Y_PHI,Y_Z,V_PLL,V_MU)
do cc=1_idef,pchunk
Y0_R(cc)=Y_R(cc)
Y0_PHI(cc)=Y_PHI(cc)
Y0_Z(cc)=Y_Z(cc)
V0_PLL(cc)=V_PLL(cc)
V0_MU(cc)=V_MU(cc)
end do
!$OMP END SIMD
!write(output_unit_write,*) 'MPI',params%mpi_params%rank,'OMP',thread_num,'R0',Y_R(1)
!write(output_unit_write,*) 'MPI',params%mpi_params%rank,'OMP',thread_num,'PHI0',Y_PHI(1)
!write(output_unit_write,*) 'MPI',params%mpi_params%rank,'OMP',thread_num,'Z0',Y_Z(1)
!write(output_unit_write,*) 'MPI',params%mpi_params%rank,'OMP',thread_num,'PPLL0',V_PLL(1)
!write(output_unit_write,*) 'MPI',params%mpi_params%rank,'OMP',thread_num,'MU0',V_MU(1)
call get_fio_GCmagnetic_fields_p(params,F,Y_R,Y_PHI,Y_Z, &
B_R,B_PHI,B_Z,gradB_R,gradB_PHI,gradB_Z, &
curlb_R,curlb_PHI,curlb_Z,flagCon,hint)
if (F%FIO_E .ge. 0) then
call get_fio_GCelectric_fields_p(params,F, &
Y_R,Y_PHI,Y_Z,E_R,E_PHI,E_Z,flagCon,hint)
end if
call get_fio_vector_potential_p(params,F,Y_R,Y_PHI,Y_Z, &
PSIp,flagCon,hint)
!write(output_unit_write,*) 'MPI',params%mpi_params%rank,'OMP',thread_num,'B',B_R(1),B_PHI(1),B_Z(1)
!write(output_unit_write,*) 'MPI',params%mpi_params%rank,'OMP',thread_num,'gradB',gradB_R(1),gradB_PHI(1)
!write(output_unit_write,*) 'MPI',params%mpi_params%rank,'OMP',thread_num,'curlB',curlB_R(1),curlB_PHI(1),curlB_Z(1)
call GCEoM1_fio_p(tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, &
B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, &
gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,Y_PHI,Y_Z,q_cache,m_cache,PSIp, &
ne,ni,nimp,Te,Zeff,flagCon,flagCol,hint)
!write(output_unit_write,*) 'R',Y_R(1)
!write(output_unit_write,*) 'PHI',Y_PHI(1)
!write(output_unit_write,*) 'Z',Y_Z(1)
!write(output_unit_write,*) 'PPLL',V_PLL(1)
!write(output_unit_write,*) 'MU',V_MU(1)
!write(output_unit_write,*) 'BR',B_R(1)
!write(output_unit_write,*) 'BPHI',B_PHI(1)
!write(output_unit_write,*) 'BZ',B_Z(1)
! write(output_unit_write,*) 'gradBR',gradB_R(1)
! write(output_unit_write,*) 'gradBPHI',gradB_PHI(1)
! write(output_unit_write,*) 'gradBZ',gradB_Z(1)
! write(output_unit_write,*) 'curlBR',curlB_R(1)
! write(output_unit_write,*) 'curlBPHI',curlB_PHI(1)
! write(output_unit_write,*) 'curlBZ',curlB_Z(1)
! write(output_unit_write,*) 'dt',params%dt
!write(output_unit_write,*) 'RHS_R',RHS_R(1)
!write(output_unit_write,*) 'RHS_PHI',RHS_PHI(1)
!write(output_unit_write,*) 'RHS_Z',RHS_Z(1)
!write(output_unit_write,*) 'RHS_PLL',RHS_PLL(1)
!write(output_unit_write,*) 'RHS_MU',RHS_MU(1)
!$OMP SIMD
! !$OMP& aligned(Y0_R,Y0_PHI,Y0_Z,V0_PLL,V0_MU,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, &
! !$OMP& RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, &
! !$OMP& k1_R,k1_PHI,k1_Z,k1_PLL,k1_MU)
do cc=1_idef,pchunk
k1_R(cc)=dt*RHS_R(cc)
k1_PHI(cc)=dt*RHS_PHI(cc)
k1_Z(cc)=dt*RHS_Z(cc)
k1_PLL(cc)=dt*RHS_PLL(cc)
k1_MU(cc)=dt*RHS_MU(cc)
Y_R(cc)=Y0_R(cc)+a1*k1_R(cc)
Y_PHI(cc)=Y0_PHI(cc)+a1*k1_PHI(cc)
Y_Z(cc)=Y0_Z(cc)+a1*k1_Z(cc)
V_PLL(cc)=V0_PLL(cc) +a1*k1_PLL(cc)
V_MU(cc)=V0_MU(cc) +a1*k1_MU(cc)
end do
!$OMP END SIMD
!write(output_unit_write,*) 'R1',Y_R(1)
!write(output_unit_write,*) 'PHI1',Y_PHI(1)
!write(output_unit_write,*) 'Z1',Y_Z(1)
!write(output_unit_write,*) 'PPLL1',V_PLL(1)
!write(output_unit_write,*) 'MU1',V_MU(1)
! call interp_fields_p(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, &
call get_fio_GCmagnetic_fields_p(params,F,Y_R,Y_PHI,Y_Z, &
B_R,B_PHI,B_Z,gradB_R,gradB_PHI,gradB_Z, &
curlb_R,curlb_PHI,curlb_Z,flagCon,hint)
if (F%FIO_E .ge. 0) then
call get_fio_GCelectric_fields_p(params,F, &
Y_R,Y_PHI,Y_Z,E_R,E_PHI,E_Z,flagCon,hint)
end if
call get_fio_vector_potential_p(params,F,Y_R,Y_PHI,Y_Z, &
PSIp,flagCon,hint)
call GCEoM1_fio_p(tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, &
B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, &
gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,Y_PHI,Y_Z,q_cache,m_cache,PSIp, &
ne,ni,nimp,Te,Zeff,flagCon,flagCol,hint)
!$OMP SIMD
! !$OMP& aligned(Y0_R,Y0_PHI,Y0_Z,V0_PLL,V0_MU,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, &
! !$OMP& RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, &
! !$OMP& k2_R,k2_PHI,k2_Z,k2_PLL,k2_MU)
do cc=1_idef,pchunk
k2_R(cc)=dt*RHS_R(cc)
k2_PHI(cc)=dt*RHS_PHI (cc)
k2_Z(cc)=dt*RHS_Z(cc)
k2_PLL(cc)=dt*RHS_PLL(cc)
k2_MU(cc)=dt*RHS_MU(cc)
Y_R(cc)=Y0_R(cc)+a21*k1_R(cc)+a22*k2_R(cc)
Y_PHI(cc)=Y0_PHI(cc)+a21*k1_PHI(cc)+a22*k2_PHI(cc)
Y_Z(cc)=Y0_Z(cc)+a21*k1_Z(cc)+a22*k2_Z(cc)
V_PLL(cc)=V0_PLL(cc) +a21*k1_PLL(cc)+a22*k2_PLL(cc)
V_MU(cc)=V0_MU(cc) +a21*k1_MU(cc)+a22*k2_MU(cc)
end do
!$OMP END SIMD
! call interp_fields_p(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, &
call get_fio_GCmagnetic_fields_p(params,F,Y_R,Y_PHI,Y_Z, &
B_R,B_PHI,B_Z,gradB_R,gradB_PHI,gradB_Z, &
curlb_R,curlb_PHI,curlb_Z,flagCon,hint)
if (F%FIO_E .ge. 0) then
call get_fio_GCelectric_fields_p(params,F, &
Y_R,Y_PHI,Y_Z,E_R,E_PHI,E_Z,flagCon,hint)
end if
call get_fio_vector_potential_p(params,F,Y_R,Y_PHI,Y_Z, &
PSIp,flagCon,hint)
call GCEoM1_fio_p(tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, &
B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, &
gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,Y_PHI,Y_Z,q_cache,m_cache,PSIp, &
ne,ni,nimp,Te,Zeff,flagCon,flagCol,hint)
!$OMP SIMD
! !$OMP& aligned(Y0_R,Y0_PHI,Y0_Z,V0_PLL,V0_MU,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, &
! !$OMP& RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, &
! !$OMP& k3_R,k3_PHI,k3_Z,k3_PLL,k3_MU)
do cc=1_idef,pchunk
k3_R(cc)=dt*RHS_R(cc)
k3_PHI(cc)=dt*RHS_PHI(cc)
k3_Z(cc)=dt*RHS_Z(cc)
k3_PLL(cc)=dt*RHS_PLL(cc)
k3_MU(cc)=dt*RHS_MU(cc)
Y_R(cc)=Y0_R(cc)+a31*k1_R(cc)+a32*k2_R(cc)+a33*k3_R(cc)
Y_PHI(cc)=Y0_PHI(cc)+a31*k1_PHI(cc)+a32*k2_PHI(cc)+ &
a33*k3_PHI(cc)
Y_Z(cc)=Y0_Z(cc)+a31*k1_Z(cc)+a32*k2_Z(cc)+a33*k3_Z(cc)
V_PLL(cc)=V0_PLL(cc) +a31*k1_PLL(cc)+a32*k2_PLL(cc)+a33*k3_PLL(cc)
V_MU(cc)=V0_MU(cc) +a31*k1_MU(cc)+a32*k2_MU(cc)+a33*k3_MU(cc)
end do
!$OMP END SIMD
! call interp_fields_p(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, &
call get_fio_GCmagnetic_fields_p(params,F,Y_R,Y_PHI,Y_Z, &
B_R,B_PHI,B_Z,gradB_R,gradB_PHI,gradB_Z, &
curlb_R,curlb_PHI,curlb_Z,flagCon,hint)
if (F%FIO_E .ge. 0) then
call get_fio_GCelectric_fields_p(params,F, &
Y_R,Y_PHI,Y_Z,E_R,E_PHI,E_Z,flagCon,hint)
end if
call get_fio_vector_potential_p(params,F,Y_R,Y_PHI,Y_Z, &
PSIp,flagCon,hint)
call GCEoM1_fio_p(tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, &
B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, &
gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,Y_PHI,Y_Z,q_cache,m_cache,PSIp, &
ne,ni,nimp,Te,Zeff,flagCon,flagCol,hint)
!$OMP SIMD
! !$OMP& aligned(Y0_R,Y0_PHI,Y0_Z,V0_PLL,V0_MU,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, &
! !$OMP& RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, &
! !$OMP& k4_R,k4_PHI,k4_Z,k4_PLL,k4_MU)
do cc=1_idef,pchunk
k4_R(cc)=dt*RHS_R(cc)
k4_PHI(cc)=dt*RHS_PHI(cc)
k4_Z(cc)=dt*RHS_Z(cc)
k4_PLL(cc)=dt*RHS_PLL(cc)
k4_MU(cc)=dt*RHS_MU(cc)
Y_R(cc)=Y0_R(cc)+a41*k1_R(cc)+a42*k2_R(cc)+a43*k3_R(cc)+ &
a44*k4_R(cc)
Y_PHI(cc)=Y0_PHI(cc)+a41*k1_PHI(cc)+a42*k2_PHI(cc)+ &
a43*k3_PHI(cc)+a44*k4_PHI(cc)
Y_Z(cc)=Y0_Z(cc)+a41*k1_Z(cc)+a42*k2_Z(cc)+a43*k3_Z(cc)+ &
a44*k4_Z(cc)
V_PLL(cc)=V0_PLL(cc) +a41*k1_PLL(cc)+a42*k2_PLL(cc)+ &
a43*k3_PLL(cc)+a44*k4_PLL(cc)
V_MU(cc)=V0_MU(cc) +a41*k1_MU(cc)+a42*k2_MU(cc)+ &
a43*k3_MU(cc)+a44*k4_MU(cc)
end do
!$OMP END SIMD
! call interp_fields_p(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, &
call get_fio_GCmagnetic_fields_p(params,F,Y_R,Y_PHI,Y_Z, &
B_R,B_PHI,B_Z,gradB_R,gradB_PHI,gradB_Z, &
curlb_R,curlb_PHI,curlb_Z,flagCon,hint)
if (F%FIO_E .ge. 0) then
call get_fio_GCelectric_fields_p(params,F, &
Y_R,Y_PHI,Y_Z,E_R,E_PHI,E_Z,flagCon,hint)
end if
call get_fio_vector_potential_p(params,F,Y_R,Y_PHI,Y_Z, &
PSIp,flagCon,hint)
call GCEoM1_fio_p(tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, &
B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, &
gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,Y_PHI,Y_Z,q_cache,m_cache,PSIp, &
ne,ni,nimp,Te,Zeff,flagCon,flagCol,hint)
!$OMP SIMD
! !$OMP& aligned(Y0_R,Y0_PHI,Y0_Z,V0_PLL,V0_MU,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, &
! !$OMP& RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, &
! !$OMP& k5_R,k5_PHI,k5_Z,k5_PLL,k5_MU)
do cc=1_idef,pchunk
k5_R(cc)=dt*RHS_R(cc)
k5_PHI(cc)=dt*RHS_PHI(cc)
k5_Z(cc)=dt*RHS_Z(cc)
k5_PLL(cc)=dt*RHS_PLL(cc)
k5_MU(cc)=dt*RHS_MU(cc)
Y_R(cc)=Y0_R(cc)+a51*k1_R(cc)+a52*k2_R(cc)+a53*k3_R(cc)+ &
a54*k4_R(cc)+a55*k5_R(cc)
Y_PHI(cc)=Y0_PHI(cc)+a51*k1_PHI(cc)+a52*k2_PHI(cc)+ &
a53*k3_PHI(cc)+a54*k4_PHI(cc)+a55*k5_PHI(cc)
Y_Z(cc)=Y0_Z(cc)+a51*k1_Z(cc)+a52*k2_Z(cc)+a53*k3_Z(cc)+ &
a54*k4_Z(cc)+a55*k5_Z(cc)
V_PLL(cc)=V0_PLL(cc) +a51*k1_PLL(cc)+a52*k2_PLL(cc)+ &
a53*k3_PLL(cc)+a54*k4_PLL(cc)+a55*k5_PLL(cc)
V_MU(cc)=V0_MU(cc) +a51*k1_MU(cc)+a52*k2_MU(cc)+ &
a53*k3_MU(cc)+a54*k4_MU(cc)+a55*k5_MU(cc)
end do
!$OMP END SIMD
! call interp_fields_p(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, &
call get_fio_GCmagnetic_fields_p(params,F,Y_R,Y_PHI,Y_Z, &
B_R,B_PHI,B_Z,gradB_R,gradB_PHI,gradB_Z, &
curlb_R,curlb_PHI,curlb_Z,flagCon,hint)
if (F%FIO_E .ge. 0) then
call get_fio_GCelectric_fields_p(params,F, &
Y_R,Y_PHI,Y_Z,E_R,E_PHI,E_Z,flagCon,hint)
end if
call get_fio_vector_potential_p(params,F,Y_R,Y_PHI,Y_Z, &
PSIp,flagCon,hint)
call GCEoM1_fio_p(tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, &
B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, &
gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,Y_PHI,Y_Z,q_cache,m_cache,PSIp, &
ne,ni,nimp,Te,Zeff,flagCon,flagCol,hint)
!$OMP SIMD
! !$OMP& aligned(Y0_R,Y0_PHI,Y0_Z,V0_PLL,V0_MU,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, &
! !$OMP& RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, &
! !$OMP& k6_R,k6_PHI,k6_Z,k6_PLL,k6_MU)
do cc=1_idef,pchunk
k6_R(cc)=dt*RHS_R(cc)
k6_PHI(cc)=dt*RHS_PHI(cc)
k6_Z(cc)=dt*RHS_Z(cc)
k6_PLL(cc)=dt*RHS_PLL(cc)
k6_MU(cc)=dt*RHS_MU(cc)
Y_R(cc)=Y0_R(cc)+b1*k1_R(cc)+b2*k2_R(cc)+ &
b3*k3_R(cc)+b4*k4_R(cc)+b5*k5_R(cc)+b6*k6_R(cc)
Y_PHI(cc)=Y0_PHI(cc)+b1*k1_PHI(cc)+b2*k2_PHI(cc)+ &
b3*k3_PHI(cc)+b4*k4_PHI(cc)+b5*k5_PHI(cc)+b6*k6_PHI(cc)
Y_Z(cc)=Y0_Z(cc)+b1*k1_Z(cc)+b2*k2_Z(cc)+ &
b3*k3_Z(cc)+b4*k4_Z(cc)+b5*k5_Z(cc)+b6*k6_Z(cc)
V_PLL(cc)=V0_PLL(cc)+b1*k1_PLL(cc)+b2*k2_PLL(cc)+ &
b3*k3_PLL(cc)+b4*k4_PLL(cc)+b5*k5_PLL(cc)+b6*k6_PLL(cc)
V_MU(cc)=V0_MU(cc)+b1*k1_MU(cc)+b2*k2_MU(cc)+ &
b3*k3_MU(cc)+b4*k4_MU(cc)+b5*k5_MU(cc)+b6*k6_MU(cc)
end do
!$OMP END SIMD
!$OMP SIMD
! !$OMP& aligned(Y_R,Y_PHI,Y_Z,V_PLL,V_MU,Y0_R,Y0_PHI,Y0_Z,V0_PLL,V0_MU)
do cc=1_idef,pchunk
if ((flagCon(cc).eq.0_is).or.(flagCol(cc).eq.0_is)) then
Y_R(cc)=Y0_R(cc)
Y_PHI(cc)=Y0_PHI(cc)
Y_Z(cc)=Y0_Z(cc)
V_PLL(cc)=V0_PLL(cc)
V_MU(cc)=V0_MU(cc)
end if
end do
!$OMP END SIMD
call get_fio_GCmagnetic_fields_p(params,F,Y_R,Y_PHI,Y_Z, &
B_R,B_PHI,B_Z,gradB_R,gradB_PHI,gradB_Z, &
curlb_R,curlb_PHI,curlb_Z,flagCon,hint)
if (F%FIO_E .ge. 0) then
call get_fio_GCelectric_fields_p(params,F, &
Y_R,Y_PHI,Y_Z,E_R,E_PHI,E_Z,flagCon,hint)
end if
call get_fio_vector_potential_p(params,F,Y_R,Y_PHI,Y_Z, &
PSIp,flagCon,hint)
call GCEoM1_fio_p(tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, &
B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, &
gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,Y_PHI,Y_Z,q_cache,m_cache,PSIp, &
ne,ni,nimp,Te,Zeff,flagCon,flagCol,hint)
!$OMP SIMD
do cc=1_idef,pchunk
vars%RHS(pp-1+cc,1)=RHS_R(cc)
vars%RHS(pp-1+cc,2)=RHS_PHI(cc)
vars%RHS(pp-1+cc,3)=RHS_Z(cc)
vars%RHS(pp-1+cc,4)=RHS_PLL(cc)
vars%RHS(pp-1+cc,5)=RHS_MU(cc)
end do
!$OMP END SIMD
!write(6,*) E_PHI
if (params%collisions) then
call include_CoulombCollisions_GCfio_p(tt,params,Y_R,Y_PHI,Y_Z, &
V_PLL,V_MU,m_cache,flagCon,flagCol,F,P,E_PHI,ne,ni,Te,Zeff&
&,nimp,PSIp,hint)
end if
end subroutine advance_GCinterp_fio_vars