Subroutine to advance GC variables
Comment this section further with evolution equations, numerical methods, and descriptions of both.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | pchunk | |||
| 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(pchunk) | :: | Y_R | ||
| real(kind=rp), | intent(inout), | DIMENSION(pchunk) | :: | Y_PHI | ||
| real(kind=rp), | intent(inout), | DIMENSION(pchunk) | :: | Y_Z | ||
| real(kind=rp), | intent(inout), | DIMENSION(pchunk) | :: | V_PLL | ||
| real(kind=rp), | intent(inout), | DIMENSION(pchunk) | :: | V_MU | ||
| real(kind=rp), | intent(in) | :: | q_cache | |||
| real(kind=rp), | intent(in) | :: | m_cache | |||
| integer(kind=is), | intent(inout), | DIMENSION(pchunk) | :: | flagCon | ||
| integer(kind=is), | intent(inout), | DIMENSION(pchunk) | :: | flagCol | ||
| type(FIELDS), | intent(in) | :: | F | |||
| type(PROFILES), | intent(in) | :: | P | |||
| real(kind=rp), | intent(out), | DIMENSION(pchunk) | :: | B_R | ||
| real(kind=rp), | intent(out), | DIMENSION(pchunk) | :: | B_PHI | ||
| real(kind=rp), | intent(out), | DIMENSION(pchunk) | :: | B_Z | ||
| real(kind=rp), | intent(out), | DIMENSION(pchunk) | :: | E_PHI | ||
| real(kind=rp), | intent(out), | DIMENSION(pchunk) | :: | PSIp | ||
| real(kind=rp), | intent(out), | DIMENSION(pchunk) | :: | curlb_R | ||
| real(kind=rp), | intent(out), | DIMENSION(pchunk) | :: | curlb_PHI | ||
| real(kind=rp), | intent(out), | DIMENSION(pchunk) | :: | curlb_Z | ||
| real(kind=rp), | intent(out), | DIMENSION(pchunk) | :: | gradB_R | ||
| real(kind=rp), | intent(out), | DIMENSION(pchunk) | :: | gradB_PHI | ||
| real(kind=rp), | intent(out), | DIMENSION(pchunk) | :: | gradB_Z | ||
| real(kind=rp), | intent(out), | DIMENSION(pchunk) | :: | ne |
subroutine advance_GCinterp_psi_vars(pchunk,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_PHI,PSIp, &
curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z,ne)
!! @note Subroutine to advance GC variables \(({\bf X},p_\parallel)\)
!! @endnote
!! Comment this section further with evolution equations, numerical
!! methods, and descriptions of both.
INTEGER,intent(in) :: pchunk
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
!! Chunk iterator.
INTEGER(ip),intent(in) :: tt
!! time iterator.
INTEGER,intent(in) :: pp
REAL(rp),DIMENSION(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(pchunk) :: k1_R,k1_PHI,k1_Z,k1_PLL,k1_MU
REAL(rp),DIMENSION(pchunk) :: k2_R,k2_PHI,k2_Z,k2_PLL,k2_MU
REAL(rp),DIMENSION(pchunk) :: k3_R,k3_PHI,k3_Z,k3_PLL,k3_MU
REAL(rp),DIMENSION(pchunk) :: k4_R,k4_PHI,k4_Z,k4_PLL,k4_MU
REAL(rp),DIMENSION(pchunk) :: k5_R,k5_PHI,k5_Z,k5_PLL,k5_MU
REAL(rp),DIMENSION(pchunk) :: k6_R,k6_PHI,k6_Z,k6_PLL,k6_MU
REAL(rp),DIMENSION(pchunk) :: Y0_R,Y0_PHI,Y0_Z
REAL(rp),DIMENSION(pchunk),INTENT(INOUT) :: Y_R,Y_PHI,Y_Z
REAL(rp),DIMENSION(pchunk),INTENT(OUT) :: B_R,B_PHI,B_Z
REAL(rp),DIMENSION(pchunk) :: E_R,E_Z
REAL(rp),DIMENSION(pchunk),INTENT(OUT) :: E_PHI
REAL(rp),DIMENSION(pchunk),INTENT(OUT) :: PSIp
REAL(rp),DIMENSION(pchunk),INTENT(OUT) :: ne
REAL(rp),DIMENSION(pchunk),INTENT(OUT) :: curlb_R,curlb_PHI,curlb_Z
REAL(rp),DIMENSION(pchunk),INTENT(OUT) :: gradB_R,gradB_PHI,gradB_Z
REAL(rp),DIMENSION(pchunk),INTENT(INOUT) :: V_PLL,V_MU
REAL(rp),DIMENSION(pchunk) :: RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU
REAL(rp),DIMENSION(pchunk) :: V0_PLL,V0_MU
REAL(rp),DIMENSION(pchunk) :: Te,Zeff
INTEGER(is),DIMENSION(pchunk),intent(INOUT) :: flagCon,flagCol
REAL(rp),intent(IN) :: q_cache,m_cache
dt=params%dt
!$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,*) 'R0',Y_R(1)
!write(output_unit_write,*) 'PHI0',Y_PHI(1)
!write(output_unit_write,*) 'Z0',Y_Z(1)
!write(output_unit_write,*) 'PPLL0',V_PLL(1)
!write(output_unit_write,*) 'MU0',V_MU(1)
!write(output_unit_write,*) 'ER',E_R(1)
!write(output_unit_write,*) 'EPHI',E_PHI(1)
!write(output_unit_write,*) 'EZ',E_Z(1)
! call interp_fields_p(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, &
call calculate_GCfields_p(pchunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, &
E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z, &
flagCon,PSIp)
!write(output_unit_write,*) 'ER',E_R(1)
!write(output_unit_write,*) 'EPHI',E_PHI(1)
!write(output_unit_write,*) 'EZ',E_Z(1)
call add_analytical_E_p(params,tt,F,E_PHI,Y_R,Y_Z)
!write(output_unit_write,*) 'ER',E_R(1)
!write(output_unit_write,*) 'EPHI',E_PHI(1)
!write(output_unit_write,*) 'EZ',E_Z(1)
call GCEoM1_p(pchunk,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,flagCon)
!write(output_unit_write,*) 'R0',Y_R(1)
!write(output_unit_write,*) 'PHI0',Y_PHI(1)
!write(output_unit_write,*) 'Z0',Y_Z(1)
!write(output_unit_write,*) 'PPLL0',V_PLL(1)
!write(output_unit_write,*) 'MU0',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,*) 'ER',E_R(1)
! write(output_unit_write,*) 'EPHI',E_PHI(1)
! write(output_unit_write,*) 'EZ',E_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 calculate_GCfields_p(pchunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, &
E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z, &
flagCon,PSIp)
call add_analytical_E_p(params,tt,F,E_PHI,Y_R,Y_Z)
call GCEoM1_p(pchunk,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,flagCon)
!$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 calculate_GCfields_p(pchunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, &
E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z, &
flagCon,PSIp)
call add_analytical_E_p(params,tt,F,E_PHI,Y_R,Y_Z)
call GCEoM1_p(pchunk,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,flagCon)
!$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 calculate_GCfields_p(pchunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, &
E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z, &
flagCon,PSIp)
call add_analytical_E_p(params,tt,F,E_PHI,Y_R,Y_Z)
call GCEoM1_p(pchunk,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,flagCon)
!$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 calculate_GCfields_p(pchunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, &
E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z, &
flagCon,PSIp)
call add_analytical_E_p(params,tt,F,E_PHI,Y_R,Y_Z)
call GCEoM1_p(pchunk,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,flagCon)
!$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 calculate_GCfields_p(pchunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, &
E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z, &
flagCon,PSIp)
call add_analytical_E_p(params,tt,F,E_PHI,Y_R,Y_Z)
call GCEoM1_p(pchunk,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,flagCon)
!$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 calculate_GCfields_p(pchunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, &
E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z, &
flagCon,PSIp)
call add_analytical_E_p(params,tt,F,E_PHI,Y_R,Y_Z)
call GCEoM1_p(pchunk,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,flagCon)
!$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
end subroutine advance_GCinterp_psi_vars